Wave propagation in a nonlinear acoustic metamaterial beam considering third harmonic generation

Nonlinear acoustic metamaterials (NAMs) provide new ways to control elastic waves. In this work, flexural wave propagation in an infinite NAM beam consisting of periodic Duffing resonators is reported by considering the third harmonic generation. Different analytical methods are proposed for the homogenized medium. By combining analytical and numerical approaches, we unveiled extensive physical properties of NAMs, including the nonlinear resonance, the effective density, the nonlinear locally resonant (NLR) bandgap, passbands, and the propagation and coupling of the fundamental and third harmonics. These characteristics are highly interrelated and they feature an identical near-field bifurcation frequency, which facilitates the prediction of functionalities. Moreover, we found that the NLR bandgap characterizes a distance-amplitude-dependent behavior that leads to a self-adaptive bandwidth in the far field. Our work will promote future studies, constructions and applications of NAMs with novel properties.

A considerable range of unusual nonlinear phenomena and physical mechanisms for nonlinear electromagnetic metamaterials (NEMs) have been demonstrated [20], including the nonlinear cloaking [21][22][23][24] and nonlinear metasurfaces [25]. Rich nonlinear dynamic effects of waves in nonlinear elastic media can also break new ground for the manipulation of phonons and elastic waves. A nonlinear AM (NAM) is defined as an AM featuring the nonlinear dynamics. Recently, Fang et al [26][27][28] found the chaotic band in finite NAMs and clarified its bifurcation mechanism. Ultra-low and ultra-broad-band wave suppression (the double-ultra effect) is achieved by NAMs based on the chaotic bands [29]. The double-ultra effect provides new vitality for the NAM field. Bridging coupling of bandgaps can manipulate the bandwidth and efficiency for wave reduction in chaotic bands [30]. In addition, multistate behavior of a wave in the bandgap is observed in NAMs [29]. Bistable periodic structures [31][32][33] also lead to interesting phenomena such as bandgap transmission.
High-order harmonic generation in nonlinear systems is useful for diverse applications. In NEMs, the second-and third harmonic generation (SHG and THG) are adopted to realize phase matching [25,34], which is of paramount importance for lasers; and enhanced nonlinear responses are desired [35,36]. High-order harmonics in nonlinear crystals/metamaterials can break the reciprocity of the wave propagation and phonon transport [37][38][39]. The broken reciprocity may leads to potential applications in realizing phononic computing and information processing [40,41]. Moreover, these nonlinear waves have been extensively applied to damage detection in materials [42], ultrasonic diagnosis and therapy in biomedicine [43][44][45].
Local resonances (LR) in AMs give rise to wave suppression in the LR bandgap. In the passbands, different from the dense modal resonances of finite AMs, resonances of the entire structure are absent in infinite AMs [46]. In the proposed infinite AMs made of side holes, Helmholtz resonators or membranes, evolution of highly intense acoustic waves is determined by weakly nonlinear (WN) processes [47][48][49], in which the bandgap shifting and SHG [50][51][52] appear. For a dispersive wave, Nayfeh and Mook [53] studied the flexural wave propagation in the infinite uniform beam placed on the weakly cubic nonlinear elastic foundation. The perturbation approach was adopted to study these WN problems.
NAMs can boost the development of various methods for achieving tunable, switchable, nonlinear, and sensing functionalities for metamaterials. However, the wave propagation in infinite NAMs and its bifurcation mechanisms are not fully understood. Especially in the strongly nonlinear regime, many physical regimes have not been clarified, including: the properties of the nonlinear effective parameters; the mechanism and distribution of high-order harmonic generation, and its interaction with the fundamental wave; the wave transmissions in bandgaps and passbands; the relationships between the nonlinear resonances, bandgaps, effective parameters and high-order harmonics; and the corresponding mathematical frameworks.
To reveal these mechanisms, this paper focuses on flexural wave propagation in a typical infinite NAM consisting of a uniform beam and periodic Duffing resonators. The finite element (FE) method and different analytical approaches are proposed to study the properties of the nonlinear effective density, ρ eff , nonlinear locally resonant (NLR) bandgap, propagations of the fundamental wave, third harmonic and their couplings. Their relationships and mechanisms arising from bifurcations are elucidated. Moreover, the distributions of harmonics, including the new behavior of the nonlinear bandgap, are reported. This work promotes further understandings of NAMs. Figure 1 shows a schematic of the NAM beam consisting of a uniform beam attached with periodic Duffing oscillators, m r , with lattice constant, a. The density, width and thickness of the primary beam are ρ, b and h, respectively. The linear density is ρ 0 =ρbh. m 0 denotes the concentrated mass fixed on the beam, and the mass m r is coupled to m 0 through a cubic nonlinear spring whose linear stiffness and coefficient of nonlinear stiffness are k r and k n , respectively. In theory, the length of the NAM bean is infinite in direction x, and a steady incident w 0 =W 0 sinωt is applied at the boundary x=0, where ω=2πf is the angular frequency, t denotes time, and W 0 is the constant amplitude. Thus, it is a half-infinite model in FE simulations. As shown in figure 1(b), a perfect match layer (PML) is built to simulate the nonreflecting boundary in the FE model. The PML consists of ten segments of uniform beams (i.e. containing no resonators) with gradually increasing structural damping. The length of each segment is 20a, and its structure parameters are identical to the primary beam in NAM. Rayleigh damping, ξ i =β i K, is adopted to attenuate the wave in the PML, where K denotes the stiffness of the beam. The stiffness damping factor for the ith segment is β i =c max (i/N p ) 3 , where c max =1 and N p =10. This PML accurately simulates the nonreflecting boundary condition. Moreover, the nonlinear strength is k W k 3 .

Theories of wave propagation
n 0 2 r s =

Effective mass density
The wave equation of a conventional medium is generally written as: where D 0 is relevant to the stiffness, u is the displacement at a point in the medium, and L represents the linear partial differential operator. For a flexural wave in a uniform beam, L=∂ 4 /∂x 4 and D 0 =Ebh 3 /12. In our AM beam, local resonators lead to a negative effective mass density, ρ eff . Based on the effective-medium theory [54,55], one obtains the wave equation of the LAM beam by a substitution of ρ into ρ eff In this case, u denotes the transverse displacement of the primary beam. For the LAM beam, substituting the representation u=Ue i(ωt−kx) gives the solution Herein, k denotes the wave number. Therefore, the eigen-function is: A homogenization approach is adopted to derive the effective parameters. In figure 1(a), w 0 and w r are transverse displacements of m 0 and m r , respectively, and u w w .
where F(t) is the concentrated force applied on m 0 . To derive the nonlinear effective density of the unit cell, let us suppose that u r =U r sinωt and F(t)=Fsinωt. Substituting these into equation (5) and adopting the first-order harmonic balance method lead to the algebraic equations: With the same approach, the approximate wave number of the NAM model is Therefore, the phase velocity of the flexural wave is expressed as Equation (11) shows that the velocity of a wave in the NAM is controllable by the amplitude, U r . Unless otherwise mentioned, the parameters specified in the simulations are: m 0 =0, m r =0.01 kg, ω r =80π rad s −1 , k n =1×10 9 N m −3 , b=0.02 m, h=0.002 m, E=70 GPa, ρ 0 =2700 kg m −3 and a=0.04 m. There are two opposite approaches to solve U r in equation (6): (i) specify force F, but the displacement W 0 is unknown; (ii) inversely, specify W 0 , but F is unknown. The results solved by specifying F=0.2 N and W 0 =2×10 -4 m, respectively, are shown in figures 2-4. As the nonlinear strength for W 0 =2×10 -4 m reaches σ=0. 19, the specified parameters can model moderate to strong nonlinearity. Nonlinear resonant curves for the two cases are as follows: there are three branches; branch-2 and branch-3 originate from the saddle-node bifurcation point whose angular frequency is ω J . A jumping occurs at ω J . As extensively demonstrated [53], branch-2 corresponds to unstable periodic solutions, i.e. the property depends mainly on branch-1 and branch-3. Although the magnitudes of U r in the two cases are approximate, ω J for F=0.2 N is much higher than ω J for W 0 =2×10 -4 m, which indicates that specifying the force has a stronger influence on the properties of the NAM model. Figure 3 illustrates the generalized effective density ρ eff /ρ 0 of the AM model. The corresponding wave number's real part, Re(k), and imaginary part, |Im(k)|, are shown in figure 4. For the linear meta-cell, the band for ρ eff /ρ 0 <0 keeps constant for the specified structural parameters; the negative ρ eff gives rise to linear LR   bandgaps whose interval is (ω r , ω cL ), in which Im(k)>0. ω cL is the cutoff frequency. Here, the linear LR bandgap is 40<f<73 Hz and ω cL =1.825ω r .
The hardening cubic nonlinearity in the NAM model generates three branches of ρ eff and k attributed to the bifurcation of U r (see figure 2). However, there are great differences between specifying F or w 0 in equation (6). If ω ∞ is defined as the frequency for ρ eff /ρ 0 →∞. ω ∞ =ω r for LAM. As indicated by figures 3 and 4, when F=0.2 N, the nonlinearity shifts ω ∞ upward to a higher frequency, and a stronger nonlinearity results in a larger deviation to ω r . Moreover, ρ eff /ρ 0 on branch-1 features unusual behaviors: ρ eff /ρ 0 originates from −∞ near ω ∞ , and ρ eff /ρ 0 →0 when ω→∞, indicating that the entire branch-1 is a negative density band, and |Im (k)|>0 in this range. If ρ eff /ρ 0 jumps to branch-3 at ω J , its |Im(k)| is still nonzero though ρ eff is positive. Therefore, the wave in NAM is attenuated when ω>ω ∞ .
In contrast, in the case W 0 =2×10 -4 m, branch-1 represents ρ eff >0 and |Im(k)|=0, and the negative ρ eff depends on branch-3. ω J initiates the band for |Im(k)|>0, and its cutoff frequency increases to ω cN in the nonlinear case with f J =58.6 Hz, ω J =2πf J . Thus, the interval (ω J , ω cN ) denotes the nonlinear LR (NLR) bandgap in the NAM model. As the deviation of ω J from ω r is much larger than the shifting of ω cL to ω cN , the width of the bandgap becomes narrower in NAMs when specifying the displacement.

THG and interaction
To analytically describe the mechanism for THG and its coupling with the fundamental wave, the meta-cell is equivalent to a homogeneous medium. In our work, the interaction concentrated force, F(t), between the oscillator m r and the primary beam is equivalent to a uniformly distributive force, f (x, t). Moreover, the linear damping effect, u , r z  of m r is also taken into consideration. Thus, the coupled wave equation of the NAM beam reads: Furthermore, u r is also a function of the space x in the homogenized media, Therefore, equation (12) can be simplified as:  (13) is the coupled flexural wave equation in the equivalent nonlinear homogeneous media. Upon considering the third harmonic, the solution for equation (13) is supposed as the superposition of the fundamental and third harmonics.
in which, H and B are the wave amplitudes, and superscript ' * ' denotes the conjugation. By substituting equation (14) into (13) and adopting the harmonic balance approach, we yield The coupling terms in equation (15) lead to complicated solutions. To describe the property of the NAM beam, we derive analytical solutions in different situations: unidirectional coupled (UC) solution, WN solution, and fully coupled (FC) solution.

UC solution
To obtain a succinct but accurate solution for equation (15), we first neglect the influence of the third harmonic on the fundamental wave but consider the energy transformation from the fundamental wave into its third harmonic, i.e. UC. Then, one obtains equations relevant to the fundamental wave only: Supposing the solution of equation (16) takes the form: Herein, k 1 denotes the wave number of the fundamental wave, and h 1 and b 1 are amplitudes. Substituting equation (17) into (16) leads to In fact, if ζ=0, equation (19) takes the same form as equation (10); furthermore, with a given h 1 in equation (18), the analytical result derived with equation (19) should be equal to equation (10) when specifying the displacement W 0 =h 1 in equation (6).
If damping is taken into consideration, by decomposing h 1 and b 1 into real and imaginary parts = + = + one can factorize the second formula in equation (18): When the third harmonic propagates along the NAM beam, the nonlinear effect arising from itself is negligible if the nonlinearity is not extremely strong. On this occasion, we obtain equations about the third harmonic by removing the nonlinear terms about B 3 in equation (15) (20), equation (21) becomes a system of linear differential equations with a forward-propagating driving term, k b e 4. k x n 1 3 i3 1 -This is also the process for THG. The solution is thereby a superposition of the particular and general solutions, i.e.
Amplitudes h 3j and b 3j are determined by the initial conditions. Because the third harmonic is absent at the boundary x=0, one yields h h . 32 31 = -Thus, there are at least two wave modes in the third harmonic whose wave numbers are 3k 1 and k 3 , respectively. For nondispersive media and nondispersive waves, 3k 1 =k 3 . However, AMs are dispersive media and flexural waves are also dispersive, we cannot directly specify 3k 1 =k 3 . Substituting equation (22) This formula indicates that k 3 depends on the amplitude of the fundamental wave b 1 . Furthermore, a conversion of equation (23) leads to an explicit solution: ( )

WN solution
When the nonlinearity is weak, the nonlinear terms relevant to the fundamental wave in equation (15) In the WN solution, linear effect is considered only so that k k 3 . perfect phase matching in this case is unrealizable. If the nonlinearity is not weak, the peak frequency of h 30 not only depends on the attached mass ratio but also on the amplitude. The UC approach neglects influences of the third harmonic on the fundamental wave, and the nonlinear effects arising from the third harmonic itself. To consider all these couplings, this work derives the FC solution.

FC solution
To obtain the FC solution, damping is not considered. Based on the analyses above, the solution for the differential equation (13) is written as: sin  sin 3  sin 3  3  sin  sin 3  sin 3  3  .  30   1  1  3 1  3  3 2  1   r  1  1  3 1  3  3 Similarly, with the boundary condition at x=0, one yields h h . 32 31 = -Expression u r 3 is reduced to In addition to the modes k 1 , k 3 and 3k 1 introduced in equation (30), the first-order modes k 3 -2k 1 and 4k 1 -k 3 , and the third-order modes 2k 3 -3k 1 and 6k 1 -k 3 appear in equation (31). If k 3 =3k 1 , these modes will be reduced to k 1 and k 3 . This manuscript considers the more general case, k 3 ≠3k 1 , in the strongly dispersive medium. With substitution of equation (30) into (13), and upon adopting the harmonic balances of the three modes, k 1 , k 3 and 3k 1 , one yields a system of algebraic equations for the FC solution: By combining the condition h 32 =−h 31 , six unknown parameters are solvable by specifying h 1 . We mention that the boundary condition h 32 =−h 31 at x=0 defines a half-infinite medium only for the propagation of third harmonic waves. It is still an infinite NAM for the propagation of fundamental waves in theories above.

Simulations and analyses
To demonstrate the analytical methods and investigate the properties of wave propagation in the NAM beam, a COMSOL FE model composed of 300 cells and a PML is established, as shown in figure 1(b). Moreover, a finite FE model consisting of 1500 cells is also established for necessary comparison. In the FE models, the incident transverse wave at the boundary x=0 is W 0 ·sin(2πf F ·t), where W 0 =2×10 -4 m, and f F is the frequency of the fundamental wave. This is a displacement boundary. Although the force is used in calculating the effective mass density, it cannot be used in time-domain simulations, on which occasion the response of the free structure approaches infinity. Unless otherwise specified, other parameters are same as those for figure 2.
Typical profiles for wave propagation in the NAM beam are shown in figure 5. In time-domain, although the incident wave is steady, the profiles of the propagating waves are not steady. The analyzed time interval at the nth cell is (t st,n , t ct,n ). We keep Δt=t ct,n − t st,n for different points equal, however, the starting time, t st,n , increases with the propagation distance, x=na.

Frequency components in harmonics
The frequency components in waves should be clarified in advance of analyzing the property of wave propagation. Figure 6 shows the typical spectra of waves in the passband ( f F =35 Hz) and LR bandgap ( f F =50 Hz). As indicated by the frequency spectrum |P( f )|, high-order harmonics become remarkable as the propagation distance na increases; moreover, the third and fifth harmonics are main components among highorder waves. However, they do not always localize at frequencies 3f F and 5f F but also at 3.5f F and 5.5f F . Therefore, the peak amplitude of |P( f )| in the interval, 2.5f F − 3.5f F , is regarded as the amplitude of the third harmonic A 3 . A t is the maximum amplitude of the wave in the analyzed time interval. A 1 is the amplitude of the fundamental wave. A p denotes the peak amplitude of |P( f )|, and f p is the corresponding frequency.
The distribution of the generalized peak frequency, Ω=f p /f F , against na and f F is shown in figure 7. For the analyzed nonlinear strength, the results indicate that: the fundamental wave is generally the main component for f F in the passbands; in the frequency range for highly efficient THG (described later), Ω=3-3.5; however, in the range of 60-70 Hz insides the LR bandgap, the second harmonic is also obvious.  Both f F =35 Hz and f F =50 Hz come from the band for efficient THG. As illustrated in figures 6(b), (d), A p A 1 for both LAMs and NAMs; A 3 remains almost steady for n>100, and A p =A 3 after a certain distance n c where the third harmonic becomes the main component of the wave. Because of the bandgap effects on the fundamental wave, n c =10 for f F =50 Hz. However, n c =290 for f F =35 Hz.

Fundamental wave, TGH and their interactions
This section describes the properties of the fundamental and third harmonics by combining the analytical solutions with the numerical method without the damping effect. In the analytical approaches, real solutions are solved and the input amplitude of the fundamental wave is h 1 =W 0 =0.2 mm, and h 3 =|h 31 |. Figure 8 shows the relationship of h 3 -f F , and the corresponding wave numbers k 1 and k 3 are shown in figure 9.
In the WN model, k 1 (3ω)=k 3 (ω) and its curves feature no bifurcations. This model anticipates that THG has two peaks at ω F =ω r /3 and ω F ≈ω r , respectively. Near ω F ≈ω r , the incident wave gives rise to the fundamental nonlinear resonances of the oscillators, m r , that efficiently generate the third harmonic. At ω F =ω r /3, the THG arises from the third superharmonic resonances of m r . However, the efficiency for THG near ω r /3 is much lower than the band near ω r in theory. As indicated by Im(k j ) in figure 9(a), the third harmonic 3ω F with wave number 3k 1 is stopped for ω r <ω F <ω cL (the linear LR bandgap), and 3ω F with wave number k 3 is stopped if ω r <3ω F <ω cL . Therefore, the range ω r /3<ω F <0.61ω r (i.e. 13.5<f F <24.3 Hz here) is a valley interval for h 3 (see figure 8(a)).
There are three branches of THG solved using the UC model, as shown in figure 8(a). Because W 0 =h 1 is specified in theory, k 1 derived with ρ eff (described by equation (10)) is equal to the one solved using the UC  In the legends, 'R300, Max' represents the maximum A 3 for cells 2<n<300 in the FE model consisting of 300 cells and PML; 'R300, Aver.' is the average A 3 in that space; 'R1500, Max' represents the maximum A 3 for 2<n<300 in the finite FE model consisting of 1500 cells. ω J denotes the saddle-node bifurcation point. model. The curves in figure 9(a) confirm this result. Therefore, there is a one-to-one correspondence between each branch of h 3 in the UC solutions and each branch of the nonlinear resonance as illustrated in figure 2(b), so that they obtain identical bifurcation frequencies, ω J (58.6 Hz here). Accordingly, the nonlinear resonance in a unit cell is helpful to explain the properties for THG in the entire NAM. The UC model still predicts the peak of h 3 near ω F =ω r /3, and its valley region ω r /3<ω F <0.61ω r agrees well with the WN result. However, there is not an obvious peak near ω r on branch-1 of the UC solution. Instead, branch-1 gradually increases to a saturated value after the valley. According to the nonlinear resonance, branch-2 is the unstable solution, thus, branch-3 may dominate the characteristic of THG for ω F >ω J . h 3 may jump from branch-1 to branch-3 at ω J leading to a reduction in the THG efficiency. Moreover, it is surprising that a peak for h 3 appears on branch-3. As indicated by Im(k 1 ), the band ω J <ω F <ω cN is the NLR gap of the third harmonic with wave number 3k 1 and the fundamental wave. Im(k 3 ) presented by the UC solution coincides with the WN solution here.
Before the statements of the FC solution, we describe the numerical results based on the FE models, as shown in figure 8. There are unavoidable numerical errors arising from the fast Fourier transform when calculating a spectrum to obtain A 1 and A 3 in the FE models. Although there are significant differences between the numerical results and the UC solutions, their tendencies, jumping frequency and peaks of curves are consistent. The FE results verify the peak near ω r /3 and the valley at 14-30 Hz. For f F >30 Hz, h 3 gradually increases to a saturate value and features a highly efficient THG band in 35<f F <53 Hz. However, there is no obvious peak as predicted by the WN solution, thus, the WN model fails in describing the features behind the valley. At 53 Hz, h 3 suddenly jumps to a low-amplitude trajectory, which accurately reproduces the bifurcation frequency, ω J . Therefore, ω J not only initiates the NLR gap of fundamental waves, but is also the cutoff frequency for highly efficient THG near ω r . Moreover, the numerical results on this trajectory follow the same regularity as branch-3 in the UC solutions, including the peak near 75 Hz, which demonstrates that branch-3 dominates the property after the jump. To further verify that the peak near 75 Hz is an essential characteristic of the NAM, we calculate the response of the FE model consisting of 1500 cells and analyze the signal in the space a<na<300a. These results confirm this peak again. However, its amplitude is much lower than that predicted in the UC solution. Numerical analyses demonstrate the validity of the UC solutions for predicting the features of wave propagation and THG.
The FC solution considers more interactions. As shown in figures 8(b) and 9(b), it characterizes complicate bifurcations. There are three backbone curves that follow the same laws as the UC solution. We also correspondingly labeled them as branch-1/2/3. However, the FC model presents a higher h 3 that is closer to the numerical results. Except for other multiple branches, the properties of k 1 indicated by branch-1/2/3 in UC and FC solutions are approximately equal. For branch-3, they predict that the propagating third harmonic's mode is k 3 . However, the behavior of k 3 on branch-1 is different. For ω F >ω r /3, the entire branch-1 in the FC solution corresponds to Im(k 3 )>0, leading to attenuation of the third harmonic with wave number k 3 . In band 35-53 Hz for highly efficient THG, both modes 3k 1 and k 3 for the third harmonics can propagate. Moreover, there are multiple branches in the FC solutions, in which two large-amplitude branches originate from the bifurcation near ω r /3, but the low-amplitude branch has k 3 >600, indicating slow waves or attenuated waves. Generally, the FC model fits the numerical result better, however, its complex bifurcations make the analyses difficult.  Figure 10 illustrates the distributions of the generalized amplitudes, log 10 (A i /W 0 ), produced by the halfinfinite FE model. The energy is efficiently converted to the third harmonic for 35<f F <53 Hz. In most cases, as the propagation distance n increases, A 3 remains steady, except for inevitable fluctuations, i.e. the vast majority of the third harmonic generates in the near-field. There are exceptional frequencies, f F <20 Hz, for which A 3 increases with the propagation distance.
In contrast, significant changes of A t and A 1 are observed as n increases. More details are presented in figure 11. For the LAM model, A t and A 1 in the bandgap remain almost steady for n>20. For the specified group of parameters, one obtains ω J <ω cL for a NAM beam. In the near-field when n<20, the NAM generates an attenuation effect for fundamental waves only in the band when ω J <ω F <ω cL , indicating numerically that ω J initiates the near-field NLR bandgap. In this tunable gap, both the fundamental and the third harmonics with the frequency ω F are stopped. However, the NLR bandgap does not remain constant as the distance increases, instead, A 1 in ω r <ω F <ω J reduces with the propagation distance, which results in the range (ω r , ω cL ) becoming the bandgap in the far field. The reason for this behavior is due to the amplitude-dependent property of the dispersion relation. For ω r <ω F <ω J , as the propagation distance increases, the attenuation of the fundamental wave weakens the nonlinear strength, then, a weaker nonlinearity actually leads to a lower bifurcation frequency, ω J ; conversely, the lower ω J reduces A 1 further. This process then leads the far field ω J to approach to the beginning frequency of the linear bandgap, ω r . In short, ω J depends on the fundamental wave's amplitude, and the amplitude depends on the propagation distance. Such a distance-amplitude-dependent property of ω J results in a highly self-adaptive behavior for the NLR gap's width to the propagation distance.  However, for ω r <ω F <ω J , A t does not decreases as much as A 1 due to the high-amplitude third harmonics. Therefore, the near-field ω J is still the initial frequency for the far field nonlinear bandgap, if we observe A t . A shallow valley of A t near ω r in the far field is also attributed to the self-adaptive bandgap. In the passbands of the NAM beam with k n =1×10 9 , the wave propagation approximates the linear case.

Influences of the nonlinear strength
The nonlinear strength greatly influences the properties of NAM [28] because h 1 or k n determines the bifurcations. As verified [27], laws arising from changing h 1 and k n are same.
There is an interesting peak for THG on branch-3. As shown in figure 12, k n influences the peak's frequency. At a given frequency, h 3 also features three branches when k n changes, and a saddle-node bifurcation, k nJ , connects branch-2 and branch-3. Branch-1 and branch-2 for different frequencies are approximate. Because f F =75 Hz is near the cutoff frequency of the linear LR bandgap, the peak will not appear at 75 Hz no matter how strong the nonlinearity. Upon slightly increasing f F to 76 Hz, a distinct peak germinates at a weak nonlinearity, k n =1.1 × 10 8 . If f F >76 Hz, the peak always appears on branch-3, and a higher frequency requires a larger k n to generate that peak. Moreover, if k n <5e7, this peak disappears, as predicted by the WN solution. Damping in the oscillators has a major influence on the amplitude of the peak but not on other parts of the branch; unexpectedly, the peak at a lower frequency obtains a larger reduction, so that the peak at 76 Hz for ζ=0.02 even disappears.
Apparently, a stronger nonlinearity indicates a higher bifurcation frequency, ω J , for the nonlinear fundamental resonance and amplitude-frequency curves of THG, as shown in figure 13. As predicted by the analyses above, a higher ω J leads to a broader band for highly efficient THG inside the bandgap. The numerical results illustrated in figure 13 demonstrate this prediction. For the case N2: k n =1 × 10 N m −3 , σ=1.9, ω J >ω cL ; therefore, the entire LR bandgap is full of third harmonics attributed to the strong nonlinearity. The jumping frequency presented in the numerical results fits the UC solutions, although there are errors. Moreover, k n barely influences branch-1 for ω F >ω cL . For ω F <ω J , a weaker nonlinearity results in a smaller h 3 on branch-1, including the peak at ω r /3; meanwhile, the curve near ω r becomes steeper. However, a peak like the WN solution never appears, although the nonlinearity is very weak. Figure 14 shows the distributions of generalized amplitudes A t , A 1 and A 3 for k n =1 × 10 10 N m −3 . This strongly nonlinear case features highly efficient THG in 35-80 Hz. The distance-amplitude-dependent bandgap appears again. Because ω J >ω cL , this NAM beam does not have a near-field (n<20) bandgap for the fundamental wave in 40-73 Hz. A longer propagation distance is required to open that bandgap than the case k n =1 × 10 9 .
However, case N2 obtains obvious reductions of the fundamental waves for 75-100 Hz (the linear passband), as detailed in figure 15. Im(k 1 ) shown in figure 16(a) explains its mechanism. As clarified above, the interval (ω J , ω cN ) dominated by branch-3 is the NLR bandgap of fundamental waves, in which Im(k 1 )>0. In theory, there is still a narrow NLR bandgap in the interval, 95.6<f F <102.4 Hz, in the case N2, permitting attenuation of the fundamental waves of 75-100 Hz in the far field.

Conclusions
This work reports the wave propagation in an infinite NAM beam (consisting of periodic Duffing oscillators) considering the THG. The characteristics and bifurcation mechanisms of the nonlinear resonance in a cell, the effective density, wave numbers, NLR bandgap, the propagation and coupling of the fundamental and third harmonics, including the interrelationships between them, were elucidated.
First, analytical mathematical frameworks were proposed to describe the fundamental wave propagation, THG and their interactions in the equivalent homogenous medium. FE models demonstrate that both the UC and FC solutions correctly describe the amplitude-frequency relationship of the third harmonic. Second, the distribution and bifurcations of the third harmonic were clarified. The THG efficiency has peaks near nature frequency ω r , ω r /3 and a frequency on branch-3. The jumping frequency, ω J , for THG inside the linear LR bandgap initiates the near-field NLR bandgap of fundamental waves, whose width is narrower for a stronger nonlinearity. Moreover, we found an interesting distance-amplitude-dependent property of the NLR bandgap, which leads to the self-adaptive bandwidth: the far field NLR bandgap approaches the original linear LR gap. In addition, when ω J <ω cL , the wave propagation in the passbands of LAMs and NAMs are similar. This law is different from the chaotic band of finite NAMs in which the broadband resonances are reduced greatly [29].

Outlook
Nonlinear media provide useful applications in many fields. A NAM is a novel nonlinear inhomogeneous medium that enable the manipulation of elastic waves with unconventional and subwavelength functionalities.
The third harmonic was used to image a smaller object and a faster reaction process than its fundamental wave. As shown in this work, the NAM is a high-quality third harmonic source in its NLR bandgap. This characteristic brings implications in extensive applications, such as filtering, the damage detection, the Figure 13. Third harmonic amplitude h 3 solved using the analytical and numerical approaches mainly for the strong nonlinearity, k n =1 × 10 10 N m −3 . The dashed (dotted) curves are the three branches of k n =1 × 10 9 (k n =1 × 10 8 ) solved with the UC method. ultrasonic diagnosis and therapy in biomedicine, and the detection and control of stress waves. The lowfrequency property of the NAM may expand their horizons, and the jumping bifurcation promotes flexible ways to control harmonics. Moreover, the phase matching between the fundamental and high-order harmonics can be realized by engineering the band structure of a NAM, so that the applications analogous with the nonlinear optics are anticipated. Besides, the wave interaction in the NAM enables the transmission and detection of waves with sum or difference frequencies. Furthermore, we report here the self-adaptive width of an NLR bandgap in NAMs. This interesting feature opens up new possibilities to realize the low-frequency and broadband vibration and noise control, or to design a nonlinear acoustic element, such as transistors and diodes.
In short, our work proposed effective approaches to study the NAMs and unveiled essential properties of wave propagation in the NAMs. The clarified phenomena enable future studies, constructions and applications of NAMs with novel properties.  . Imaginary parts of wave numbers k 1 and k 3 for the two cases. N1: k n =1 × 10 9 ; N2: k n =1 × 10 10 .