Emergent subharmonic band gaps in nonlinear locally resonant metamaterials induced by autoparametric resonance

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


I. INTRODUCTION
The interest in metamaterials originates from the challenge for designing a new generation of materials with superior properties, not found in nature. In mechanics, the features promoted by metamaterials are, in general, associated to negative mass, elastic (bulk, longitudinal, and/or shear) moduli or Poisson's ratio effects [1][2][3]. These negative effective properties attained by dynamical metamaterials are in many cases triggered by localized resonances and the associated strongly dispersive behavior. The most extraordinary feature of locally resonant metamaterials is the subwavelength size of their unit cells, which allows to break through the traditional limits in wave focusing and imaging, and in the size of structures promoting sound insulation [4].
Recently, the effect of nonlinearities in metamaterials has received great attention. Nevertheless, in most of the works dealing with nonlinearity, the considered metamaterials do not include locally resonant units. Instead, phononic structures involving nonlinear interactions are widely explored, both from the fundamental and applied viewpoints. In general, nonlinear interactions involve either cubic and/or quartic potentials, such as the famous Fermi-Pasta-Ulam-Tsingou (FPUT) problem [5,6], or more complex models such as the Hertz interaction law in granular materials [7][8][9][10]. Recent studies on FPUT-type phononic lattices focused mainly on * pri.brands@gmail.com weakly nonlinear interactions and the effect of nonlinearity on the dispersion relations [11][12][13][14][15][16]. In these papers, it was shown that nonlinearity induces amplitude-dependent dispersion and group velocity, i.e., demonstrating the tunability potential of nonlinear phononic structures. By making use of this feature, new functionalities have been proposed, such as mechanical switches and logic gates [17], acoustic diodes [18], and acoustic rectifiers [19]. In granular materials, the effect of highly nonlinear interactions has been shown to induce solitary wave propagation [20,21] and discrete breathers [22].
The interaction between nonlinearity and the local resonance phenomenon, i.e., at the level of the local coupling between the host and the local resonator, has been much less explored. Within the nonlinear dynamics community, several works have addressed the irreversible energy transfer mechanisms induced by a single purely nonlinear attachment, the so-called nonlinear energy sinks [23][24][25]. However, so far, very few works have considered cases in which nonlinear resonant attachments are densely or periodically distributed in a host material in a similar fashion as done in linear locally resonant acoustic metamaterials. Within the framework of acoustic metamaterials, Lazarov and Jensen [26] have used the harmonic balance method to investigate the amplitude dependence of the locally resonant band gap in a discrete lattice system including cubic interaction between the local resonators and the main chain atoms. The effect of both softening and stiffening (also known as hardening) cubic nonlinearity on the amplitude-dependent dispersion behavior of the nonlinear locally resonant lattice metamaterial has been addressed by Manimala and Sun [27]. The potential functionality of these materials as selective filters has thereby been demonstrated. Other recent publications [28][29][30] have further investigated similar lattice structures and showed that bifurcations induced by the local resonance mechanism may induce chaotic bands and influence the attenuation bandwidth.
The present paper shows that multiple attenuation zones can emerge in a nonlinear locally resonant metamaterial, thus promoting novel opportunities for wave manipulation. In particular, the emergence of subharmonic attenuation zones in locally resonant metamaterials is shown to be induced by nonlinear energy conversion. Energy exchange between propagating wave modes is well known in nonlinear wave dynamics. It originates from the nonlinear terms driving the linear waves resonantly and giving rise to energy transfer between them [31]. Rushchitscky et al. [32,33] demonstrated nonlinear wave interactions in a nonlinear material described by the Murnagham potential. Recently, nonlinear wave interactions in quadratic, cubic, and quadratic-cubic phononic lattices were investigated with the aim of allowing advanced manipulation of wave propagation in metamaterials [13,34,35]. By making use of second-harmonic generation, Guo et al. [36] designed a nonlinear metasurface with unusual reflection effects. The novelty of the present work consists in the discovered effect that wave interactions between propagating and evanescent waves may also occur induced by nonlinear local interactions in metamaterials, which constitutes a new mechanism for wave mitigation and control.
The metamaterial design under consideration here is a discrete version of the locally resonant metamaterial proposed by Liu et al. [37]. The key feature of the present model is the incorporation of a neo-Hookean nonlinear behavior, exhibiting pronounced nonlinear tension-compression asymmetry, representative of elastomeric rubber-coated inclusions [38,39]. It should be pointed out that the present work is fundamentally different from other publications [40][41][42] on soft elastomeric metamaterials, where wave tailoring is achieved by induced geometrical modifications during wave propagation, e.g., as a result of buckling [42], or from the recent work by Konarski et al. [43] in which a frequency-dependent effective medium theory for locally resonant metamaterials with hyperelastic inclusions and embedded mechanical instabilities has been derived. Although the asymmetric nonlinear constitutive behavior of inclusions was considered, the interest there was in investigating second-harmonic generation and the metamaterial response in multiple frequency regimes, i.e., from low to high frequencies.
The paper is organized as follows. In the first part, Sec. II, the mechanical model system under consideration is presented together with the results obtained by direct numerical simulations. These results suggest the emergence of a subharmonic attenuation zone and its characteristics are underlined. In the second part, Sec. III, semianalytical solutions based on the method of multiple scales are derived using an approximate quadratic nonlinear local interaction model. Results for the elastomeric and approximate quadratic metamaterials are compared. In Sec. IV, the analytical expressions for the quadratic metamaterial model are used to unravel the underlying physical phenomenon explaining the emergent subharmonic attenuation zone.

II. MECHANICAL MODEL SYSTEM OF A HYPERELASTIC LOCALLY RESONANT METAMATERIAL
The mechanical model system considered throughout the paper is based on a discrete analog of the three-dimensional (3D) structure originally proposed by Liu et al. [37] to demonstrate the local resonance mechanism, shown in Fig. 1(a). The unit cell consists of a linear matrix (host medium) with density ρ, containing an inclusion with density ρ o covered by a thick compliant coating layer which can potentially behave nonlinearly. Restricting attention to longitudinal waves propagating along one of the principal directions, and assuming that the wavelength is much larger than the unit cell size, the simplified discrete analog model system is shown in Fig. 1(b). In this model system, a unit cell is represented by a discrete main chain mass m, given by ρV , where V is the volume of the matrix material, connected to a local mass m o , equal to ρ o V o , with V o being the volume of the heavy inclusion. The interaction between the unit cells is described by a linear spring of stiffness k, which is proportional to the matrix elastic modulus and, possibly, viscous damping c. Considering that the density and volume of the inclusion are much larger than those of the compliant rubber layer, the interaction between the main chain mass m and the local mass m o within the unit cell is described by a massless dash-pot with damping c o and an elastic interaction force f o , resulting from the parallel connection between a connected massless linear spring with stiffness k o , which is proportional to the equivalent Young's modulus of a neo-Hookean solid in uniaxial extension, and a nonlinear massless spring with force function f NL .
For an arbitrary unit cell j, with j = 1, . . . , n cell , n cell ∈ Z being the number of unit cells, the nondimensional equations of motion in terms of the nondimensional absolute displacementsū j = u j /u ref andv j = v j /u ref of the main chain and the local masses, respectively, at unit cell j and at neighboring unit cells j − 1 and j + 1 are given by: , relative to the main chain and oscillator, respectively;f NL is a nondimensional force function in terms of q j = v j − u j , i.e., the relative motion between the local and main chain masses. The nonlinear force results from the nonlinear constitutive material model, which will be further introduced in Sec. II A.

A. Nonlinear material model
The effect of nonlinearity on the local resonance phenomenon in the periodic chain will be investigated and compared with the classical linear local interaction. Since this work aims to incorporate the nonlinear stress-strain relation representative of rubberlike materials, commonly neglected in wave propagation analyses of metamaterials, one of the basic elastomer material models will be used, i.e., the hyperelastic incompressible neo-Hookean material model. The constitutive relation for a one-dimensional incompressible neo-Hookean material under uniaxial tension-compression is given by [44]: where σ is the nominal, or engineering, stress and λ the stretch ratio given by λ = L/L 0 = (L 0 + q)/L 0 , with L and L 0 being the deformed and undeformed lengths, respectively; C is a material constant, which for an incompressible neo-Hookean solid at small strains is related to Young's modulus force f NL = f NH , which can be written as a function of the relative displacement q of the local resonator, as follows: where k o = 6CA 0 /L 0 with A 0 being the reference cross section over which the interaction force is applied. The total local interaction force between the main chain mass and the attached local oscillator (consisting of linear and nonlinear springs in parallel) is given by: Inserting a Taylor's series expansion of the nonlinear forcing term (presented subsequently in Sec. III) in Eq. (4) shows that the equivalent linear stiffness of the locally resonant unit is 2k o and the equivalent linear resonance frequency of the oscillator In Fig. 2, the normalized neo-Hookean force-elongation relation with the corresponding linear model and the approximate quadratic model (to be used for the analytical solution in Sec. III). Note its asymmetric tension-compression behavior typical of elastomeric materials: Indeed, it exhibits a decrease in the tangent stiffness modulus in tension, hereafter referred to as "softening," and an increase in tangent stiffness modulus in compression, hereafter referred to as "stiffening."

B. Numerical model
The dynamic behavior of the metamaterial with nonlinear local resonators is first investigated via direct numerical simulations. A sufficiently long chain is considered for this purpose to eliminate finite-size effects. This choice is dictated by the fact that the highly nonlinear behavior of the local resonators prevents the direct application of Bloch's theorem [45] (which would reduce the analysis of an infinite periodic system to a single unit cell only).
The numerical model of the system under consideration consists of two 1D discrete chains connected through a central excited unit cell, as shown in Fig. 3. The chains on the left and right sides are each composed of n cell locally resonant unit cells of the type shown in Fig. 1. The total number of locally resonant unit cells in this system is 2n cell + 1. A prescribed harmonic displacement, in the form: is applied to the main chain mass of the central unit cell at position j = 0. In (5), ε 1 is a small nondimensional parameter,û e is the nondimensional amplitude scaled to be = t ref is the nondimensional input angular frequency, and φ the phase angle. Note that for each simulation the harmonic displacement excitation at a fixed frequency has been prescribed for a sufficiently long time interval to capture long timescale effects, and no frequency sweep response has been studied in this work. On the left and right ends, the chain is connected to fixed nodes by means of linear springs and dash-pots with a critical damping value c cr = 2 √ km. The default values for the nondimensional material and geometrical parameters used in the numerical simulations are as follows:ω m = 1, α = 2.5 × 10 −5 , β = 0.1, ζ m = 0, and ζ o = 0.05. Then, by considering appropriate reference values the link to (any) real material can be made. The number of unit cells in the left and right chains is n cell = 5000. Hence, the total number of degrees of freedom in the system amounts to 2 × (2n cell + 1) = 20 002. The ordinary differential equations (1) describing the momentum balance of this system are integrated numerically using an explicit time integration method, i.e., the fourth-order Runge-Kutta method (ODE45, MATLAB). The total integration time is 5000 dimensionless time units. A constant time step of 0.10 is used in order to ease the fast Fourier transform (FFT) during postprocessing for the investigation of the frequency spectrum of the system response.

C. Evidence of an emergent phenomenon
Periodic structures with linear locally resonant units are characterized by the appearance of an attenuation frequency zone in which waves do not propagate or are highly evanescent, i.e., a band gap [1,3,46]. The negative effective dynamic mass [47] is the mechanism that drives the formation of such an attenuation zone in structures exhibiting dipolar resonance. For a linear locally resonant metamaterial, the band gap occurs aroundω R , or¯ /ω R = 1, i.e., the resonance frequency of the oscillator, whereby the band gap width is independent of the amplitude of the applied excitation.
First, as a reference, numerical simulations for wave propagation in a linear locally resonant metamaterial are performed. Normalized and scaled input amplitudesû e = [1, 3,5] are applied, as given in Eq. (5). In Fig. 4(a), the results of the numerical simulations are presented in terms of the displacement transmissibility evaluated in a root mean square (RMS) sense, i.e., the ratio between the RMS amplitude values of the displacement averaged over a number of periods (at least 7, at the lowest frequency of excitation, and at most 84, at the highest frequency of excitation) evaluated at the main chain masses located at cell 1000 and cell 0 (reference value). Transmissibility is calculated for a number of normalized input frequencies¯ /ω R between 0 and 4. As expected, the results for the linear case are amplitude independent and consistent with the analytical predictions of the attenuation zone around¯ /ω R ≈ 1.
Next, the chain with the neo-Hookean local interaction, as expressed in Eq. (3), is evaluated. The results of the numerical simulations in terms of RMS displacement transmissibility are shown in Fig. 4(b). The locally resonant band gap around ≈ 1 is still present independently of the value of the input displacement amplitude. However, in contrast to the linear case, three marked differences emerge: (1) the level of displacement attenuation slightly decreases with the increase of the input displacement amplitude; (2) the band gap slightly shifts to the right (i.e., to higher frequencies) with increasing input amplitude; (3) a second attenuation zone emerges around¯ /ω R ≈ 2 for higher excitation amplitudes.
The first observation is an inherent feature of a nonlinear system, periodic or not, i.e., revealing an amplitude-dependent response. The second observation is also a consequence of the amplitude-dependent behavior, but it is due to the competition between "softening"-"stiffening" in the tension-compression behavior. In a chain with positive cubic nonlinear local interaction, the "stiffening" effect in both tension and compression has previously been shown to induce a shift of the band gap position toward higher frequencies [26]. In the case of asymmetric neo-Hookean local interaction, both "stiffening" and "softening" effects are present (see Fig. 2). In Figs. 5(a) and 5(b), the response of the local oscillators through the chain as a function of the excitation frequency is shown for both linear and nonlinear local interactions, respectively. In these figures, the response of the local oscillators is obtained by quantifying the oscillator amplification, i.e., the average peak-to-peak response of the local oscillators normalized by the peak-to-peak input amplitude. In the linear case, Fig. 5(a), the amplification peaks occur at the local resonance frequency, while, in the nonlinear case, Fig. 5(b), the amplification peaks vary along the chain and occur at frequencies above the local resonance one, following a stiffening backbone trend. This explains the shift in the band gap position toward higher frequencies with an increase of the input amplitude in the transmission plots in Fig. 4(b). Indeed, Fig. 5(b) shows that, FIG. 5. Averaged peak-to-peak response of local oscillators in the metamaterial normalized by the peak-to-peak input amplitude forû e = 5 and ε = 0.1 in (5) due to the nonlinearity, the response of each local resonator follows a backbone curve, whose maximum amplification, as well as the associated frequency, decrease further away from the excitation. This is because the wave attenuates while it propagates due to damping and local resonance effects. Notice that the width of the first and second attenuation zones, indicated by the dashed lines in Fig. 4(b), are associated with the width of the backbone behavior, indicated by the dashed lines in Fig. 5(b).
To the best of the authors' knowledge, the third, emergent, phenomenon has not yet been reported in the literature. Therefore, in the remainder of this paper, the focus will be on the second attenuation zone. Indeed, additional simulations (not presented here) revealed that the emergent phenomenon occurs regardless of the damping value. Since the focus of this paper is on the mechanism responsible for the emergence of the second attenuation zone, a parametric analysis of the phenomenon is outside the present scope. More insight in the characteristics of this attenuation zone can be obtained by analyzing the numerical results for a single excitation frequency in the time domain as well as in the frequency domain, i.e., after applying the FFT to the time signature.
Space-time maps of instantaneous total energy at several unit cells of the nonlinear locally resonant metamaterial excited at¯ /ω R = 2.47 are plotted in Fig. 6. For low excitation amplitudes, i.e., when the emergent second attenuation zone does not occur, as for instance forû e = 1, energy propagates with approximately constant velocity. In Fig. 6(b), the space-time map shows the evolution of the instantaneous total energy forû e = 5, for which the second transmission dip was observed around this frequency. A space-time attenuation zone can be clearly identified as denoted by the letter A. Notice that it occurs after a transient propagative regime and far away from the excitation. The propagative regime, identified as P, corresponds to an initial transient state. This regime is followed by a second regime of intense vibration and slower wave propagation, denoted as S. The envelope of this zone shows a decreasing velocity profile, as depicted in Fig. 6.
Further insights regarding the dynamics in zone S that will lead to attenuation far away from the excitation point are provided in Fig. 7, where the responses of the main chain mass and the local oscillator at unit cell 1, within zone S, are compared. The analysis is performed after sufficiently long time to avoid transient effects. In Fig. 7(a), notice that the local oscillator shows a higher amplitude of oscillation compared to the main chain mass. The Fourier transform of the computed time data as τ → ∞ [ Fig. 7(b)] shows that the highest harmonic amplitude for the local resonator occurs at the one-half subharmonic of the primary input frequency, i.e., at 1/2¯ . The response of the local resonator also shows an important dc (static) component and smaller contributions of other harmonics. The dominant frequency of the main chain oscillation, on the other hand, still corresponds to the primary input frequency¯ /ω R 2. A smaller amplitude contribution is observed in the main chain at 1/2¯ . In the following sections, these observations are used to develop an approximate semianalytical model, with the aim of unraveling the mechanism underlying the second attenuation zone formation and the role of the nonlinear local interaction of neo-Hookean type.

III. SEMIANALYTICAL ANALYSIS OF QUADRATIC LOCALLY RESONANT METAMATERIAL
In this section, a perturbation method, the method of multiple scales [48], is used to shed light on the dynamical behavior of the nonlinear locally resonant metamaterial considered in this paper with the main goal to clarify the mechanism underlying the emergent second attenuation zone.
Perturbation methods consist in obtaining the response of a nonlinear system by perturbing the response of the corresponding linear system. Following the general procedure of perturbation methods, the local interaction force f o given by Eq. Taylor's series expansion of the nonlinear interaction force gives: which approximates a linear function with stiffness k o (a property well known for the neo-Hookean model). In this case, the approximation up to order q 3 j of the local interaction (4), with the nonlinear neo-Hookean model and a linear spring parallel to it, is given by: with γ 1 = 2 and γ 2 = −1/L 0 , following from the Taylor's series expansion in Eq. (6). Then, in accordance with the equivalent linear stiffness shown in Sec. II A, the equivalent linear resonance frequency of the oscillator is ω R = √ 2k o /m o . The nondimensional form of the quadratic approximation (truncated after the second term) of the neo-Hookean interaction is plotted in Fig. 2.
Without loss of generality, the forced wave propagation problem to be solved via the method of multiple scales is given in nondimensional form by: subject to boundary condition:ū 0 (τ ) =ū e (τ ), initial conditions:ū j (0) = 0,u j (0) = 0, whereū e (τ ) is the harmonic displacement applied at j = 0, The method of multiple scales seeks for a solution of the nonlinear dynamical problem in the form of a spatially and temporally uniformly valid solution through an asymptotic expansion. In the case of forced wave propagation problems, multiple time and space scales should be considered [49]. For first-order approximated solutions, these scales are defined as: where ε 1 is a small scale separation parameter which makes each time (respectively, space) scale slower (respectively, longer) than its predecessor. Then, the displacement solution for an arbitrary unit cell j (ū j ) can be expressed by an asymptotic expansion of the form: whereû (n) j =û (n) j (X 0 , X 1 , T 0 , T 1 ) is the normalized displacement at the order n.
Thus, within the framework of the method of multiple scales, the approximate solution to the forced wave propagation problem requires the solution of an initial boundary value problem at each order of perturbation (see Ref. [49] for further details).
The numerical results shown in Fig. 7 suggest interactions between the excitation frequency (i.e., the primary frequency) and its half subharmonic (or secondary frequency), within the local resonance band gap. Indeed, the analytical analysis of the forced harmonic response of a local resonance metamaterial demonstrates that the local resonant mode is excited in the transient regime. Thus, at the first-order of perturbation, we assume that the localized excitation (boundary condition) is given by a time-harmonic displacement applied at j = 0 of the form: u e (τ ) = εû e (τ ) = ε(Û e e iω p τ +φ ep + εÛ e e iω s τ +φ es ). (11) where i is the imaginary unit, the primary frequency corresponds to the angular frequency of excitation, i.e.,ω p ≡¯ , ω s is the secondary frequency,Û e is the magnitude of the excitation, and φ e p and φ e s are, respectively, the phases of the primary and secondary harmonics. Note that the amplitude of the secondary harmonic is taken much smaller than the primary one. Accordingly, the zeroth-order solutionû (0) j in Eq. (10) is given by a superposition of the wave modes associated with the primary and secondary angular frequencies 063003-7 FIG. 8. Dispersion relation for the approximate linear undamped locally resonant metamaterial depicting the wave modes selected for the wave-wave interaction analysis. Parameters of the model system: as follows: where p and s are the vectors of wave mode shapes corresponding to the linear primary (ω p , μ p ) and secondary (or subharmonic) (ω s , μ s ) waves, respectively, with μ p and μ s the corresponding nondimensional wave numbers, and c.c. stands for the complex conjugate. Herein, we consider the primary wave mode (ω p , μ p ) induced by the excitation located on the high-frequency acoustic branch of the corresponding linear metamaterial and its angular frequency satisfies a 2:1 relation with respect to the linear locally resonance frequency, i.e., where σ p is a primary frequency detuning parameter. Consequently, the half subharmonic frequency might be related to the local resonance frequency bȳ with σ s > 0 being a subharmonic frequency detuning parameter. Due to the proximity to the local resonance frequency, a subharmonic wave mode (ω s , μ s ) of evanescent nature can be induced. For the purpose of illustration, the dispersion diagram of the corresponding linear metamaterial depicting the primary propagating wave mode and an associated evanescent subharmonic wave mode is shown in Fig. 8. The primary and subharmonic wave modes to be considered within the multiple scales framework correspond to those of the approximate linear metamaterial. Therefore, the relation between the corresponding frequency and wave number is given by: (15) From Eqs. (13), (14), and (15), a closed-form relation for the analogous wave-number detuning parameters can be expressed in terms of the angular frequencies, providing the corresponding dispersion relations for the approximate linear metamaterial.
Following the multiple scales analysis for forced wave propagation problems as described by Silva et al. [49], a set of 1D advection-reaction equations with nonlinear source (reaction) term describing the slow and long wave modulation of the zeroth-order solution is obtained as follows: where σ μ i = ε −1 Im(−2μ s ). In terms of real coefficients, this system of advection-reaction equations can be written in matrix form as: wherê with V and Q being the nondimensional advection and damping matrices, respectively, andq NL the nonhomogeneous nonlinear vector. Herein,ĉ r and α r are, respectively, the damping and nonlinear interaction coefficients relative to the primary (r = p) and secondary (or subharmonic) (r = s) wave modulations; ν p is the velocity of the envelope of the primary wave; ν s is the purely complex velocity of the envelope over the evanescent secondary wave, which means that for the evanescent wave there is no pure convection but a phase change under convection;Û (r) R andÛ (r) I are, respectively, the real and imaginary parts of the complex amplitudeÛ (r) relative to the primary (r = p) and secondary (r = s) wave modulations; and θ = (σ p − 2σ s )T 1 − ε −1 μ p X 1 .

Semianalytical results
In this section, we investigate the numerical solutions to the set of advection-reaction equations [Eq. (17)], which resulted from the consideration of multiple time and space scales in the dynamic resonance of the approximate quadratic locally resonant metamaterial. The 1D advection-reaction equations involve first-order derivatives in space and in time. Since the advection-reaction problem involves a positive real velocity ν p , the backward finite difference scheme is used for space discretization. On application of the boundary conditions, the resulting nonlinear initial value problem consists of a set of first-order differential equations solved using the implicit Euler integration scheme coupled with the Newton-Raphson method.
In what follows, we analyze results from the semianalytical model, assuming u ref = L 0 such thatγ 2 = −1, for primary and secondary waves with frequenciesω p = 2ω s = 2.05ω R , i.e.,ω s is located within the local resonance band gap and ω p is located on the high-frequency acoustic branch of the dispersion diagram of the linear metamaterial, as depicted in Fig. 8; the frequency detuning parameters are σ p = 0.5ω R , σ s = 0.25ω R . The system is considered initially at rest, and excited in the middle node, at j = 0, at the primary frequencȳ ω p with amplitudeÛ e . A residual excitation at the secondary frequency with amplitudeÛ (s) = εÛ e is considered to enable the nonlinear wave-wave interaction. Indeed, due to the localized nature of the excitation and the nonlinearity, a residual excitation at the secondary frequency might occur.
The slow-long wave modulations of the primary wave in space and time for two excitation amplitudesÛ e = {0.5, 1.5} are depicted in Fig. 9. For the lower amplitude of excitation, propagation of the primary wave is observed with the amplitude of the same order as the excitation amplitude close to the boundary and slowly decaying away from the excitation point for the entire simulation time. This is in accordance with the numerical results for the elastomeric locally resonant metamaterial, which exhibit high transmissibility (near to the unity) for low amplitudes of excitation and excitation frequencies outside the local resonance band gap (see Fig. 4). However, for the larger amplitude,Û e = 1.5, results from the semianalytical model show that the propagation at the level of the excitation amplitude is sustained only for a short interval of time, after which, the amplitude of the propagative primary wave decays considerably. On the contrary, the secondary wave response increases with the excitation amplitude [compare Figs. 9(b) and 9(d), and note the difference in the scales], but its effect is confined to the near-field only.
In Fig. 10, the evolution of the total energy for several unit cells is reconstituted from the semianalytical results forÛ e = 0.5 andÛ e = 1.5. As expected, for the lower amplitude of excitation, wave propagation with velocity ν ≈ 1 is observed for the entire simulation time, while for the higher excitation amplitude, the propagation with ν ≈ 1 occurs for limited time, after which only near-field excitation with very slow wave velocity is observed. Thus, for large-enough excitation and as τ → +∞, the response of the quadratic metamaterial is strongly attenuated, in accordance with the transmissibility results for the neo-Hookean locally resonant metamaterial. The difference between the numerical results for the elastomeric metamaterial (Fig. 6) and those obtained via the method of multiple scales for the simplified quadratic metamaterial model is in the spatial extent of zone S, where the response is dominated by the subharmonic wave mode. Despite that difference, the simplified semianalytical quadratic nonlinear model is able to capture the phenomenon observed previously for the more complex nonlinear interaction.
For four spatial positions along the chain (X 0 = {10, 100, 200, 300}), the time signatures of the primary and secondary wave modulations are shown in Fig. 11. For the lower amplitude of excitation,Û e = 0.5, the amplitude of the primary wave mode grows smoothly with time till the steadystate amplitude is reached. Similar behavior is observed for the secondary wave. However, the subharmonic modulation amplitudes grow slower with time and reach considerably lower steady-state values. For the higher amplitude of excita-tionÛ e = 1.5 [ Fig. 10(b)], a transient regime with high modulation amplitude occurs for the primary wave, after which the primary wave amplitude decreases until a steady-state value is reached. It is remarkable, especially in the near field, that while the primary wave amplitude decreases, the secondary wave amplitude increases. The cause for this is the energy exchange taking place between the primary and the secondary wave modes.
To investigate the energy exchange phenomenon over a range of excitation amplitudes, the steady-state amplitudes of the wave modes (ω s , μ s ) and (ω p , μ p ), i.e., lim τ →∞Û (p) and lim τ →∞Û (s) , for sufficiently long time τ → +∞, are evaluated for a range of excitation amplitudesÛ e and shown in Fig. 12 for several positions along the metamaterial. Close enough to the excitation, i.e., X 0 = 10 [ Fig. 12(a)], and for low excitation amplitudes, the system behaves as a linear one, i.e., the amplitudes of the primary and subharmonic wave modes increase linearly with the excitation amplitude, with the subharmonic growth rate much lower than the primary one, thus the amplitude of the subharmonic wave mode remains much smaller than the amplitude of the primary wave mode. Similar linear behavior for low amplitudes of excitation is observed for X 0 = {100, 200, 300}. Indeed, due to the near-field feature of the subharmonic wave mode, far enough from the excitation, i.e., for X 0 = {200, 300}, the subharmonic wave amplitude tends to zero. For all positions, the alteration of wave-wave interaction is evident when the excitation amplitude reaches a critical value,Û crit e ≈ 0.6. For larger excitation amplitudes and in the very near field, X 0 = 10, the primary wave mode amplitude reaches saturationÛ (p) =Û (p) sat , while 063003-9 the rate of the subharmonic wave amplitude growth with excitation becomes larger and nonlinear. Close enough to the excitation, the subharmonic wave amplitude may even exceed the primary wave amplitude [see Fig. 12(a)]. Then, as subharmonic wave motion becomes dominant over the primary one, and due to the nature of the subharmonic wave mode, energy is being transfered to the local oscillator.

IV. PHYSICAL INTERPRETATION: EMERGING SUBHARMONIC ATTENUATION DRIVEN BY ENERGY CONVERSION MECHANISM
Based on the above results, the second attenuation zone observed in the numerical simulations of the elastomeric metamaterial results from the energy exchange between the interacting primary (ω p , μ p ) and subharmonic (ω s , μ s ) wave modes due to the autoparametric resonance. Similarly to the second-harmonic resonance condition, a 2:1 ratio between the angular frequencies of the interacting wave modes is required. However, from the above results, for interacting waves of different natures, i.e., one being propagative and the other one, evanescent, phase matching is not required. Instead, the subharmonic frequency must be close enough to the local resonance frequency, i.e.,ω s ≈ω R , thus inducing evanescent behavior.
Rewriting Eq. (16) in terms of the wave componentsû (0,s) j andû (0,p) j , it becomes with ν p ,ĉ p ,ĉ s , α p , α s > 0, sgn(ν s ) < 0. The asymmetric nonlinear coupling termû (0,p) jû (0,s) * j appears in the evolution equation (18b) for the subharmonic wave mode amplitude. This corresponds to a parametric excitation of the subharmonic wave mode by the primary wave mode amplitude. If the amplitude of the primary wave is large enough such that the parametric interaction contribution can overcome the losses in the subharmonic wave and the internal resonant convective force, then the secondary wave amplitude may grow. On the other hand, the dynamics of the primary wave mode, Eq. (18a), might be affected by the subharmonic wave Since this is not a cross-coupling term, and the coefficient α p is positive, it can be regarded as an external sink term. Thus, as the amplitude of the secondary wave becomes sufficiently large, the nonlinear term in (18a) contributes to the saturation of the primary wave amplitude, stopping the linear amplification process. In the literature, the internal resonance phenomenon has often been used to promote energy exchange between propagating wave modes [48] and the autoparametric resonance phenomenon has typically been studied on finite systems involving few degrees of freedom or even a single nonlinear oscillator, such as mass-spring-pendulum systems and ship models [50]. Although the autoparametric resonance mechanism has been used to design simple nonlinear vibration absorbers [51], the role of this phenomenon in wave propagation and the possibility of inducing multiple attenuation zones in periodic systems has not been explored. In this paper, the focus has been restricted to the one half subharmonic; however, it is expected that similar attenuation zones may appear around frequencies corresponding to other tongues of instability of the Ince-Strutt diagram [52,53]. The conditions for stability of these multiple attenuation zones induced by autoparametric resonance are still unknown and will be the subject of future investigation.
The underlying phenomenon observed in the elastomeric metamaterial model and its approximate quadratic model shows similarities with quenching [52] in nonlinear dynamical oscillators, and the physics of optical parametric oscillator networks. In summary, forcing a lattice with nonlinear oscillators atω p ≈ 2ω R with sufficiently high amplitude, the trivial primary wave mode solution saturates and energy is transferred to the subharmonic wave mode. Due to the evanescent nature of the subharmonic wave mode (ω s , μ s ), the motion becomes localized in the oscillator, which is at resonance. Therefore, the overall energy in the main chain is reduced, giving rise to the subharmonic attenuation as demonstrated numerically in Sec. II C.

V. CONCLUSIONS
In this paper, discrete locally resonant metamaterials with nonlinear interaction between the local oscillator and the main chain mass have been investigated. The aim was to explore the effect of material nonlinearities of rubberlike materials, described by the neo-Hookean material model, on the dynamical behavior of metamaterials. Numerical simulations of a neo-Hookean locally resonant metamaterial revealed that not only the locally resonant band gap would be generated but also a half subharmonic attenuation region emerges.
In order to unravel the formation mechanism of the subharmonic attenuation zone, the dynamics of an approximate nonlinear model system involving only quadratic nonlinearity has been investigated using the method of multiple scales. Results of these analyses showed that the autoparametric resonance mechanism induced by the combination of the local resonance and quadratic nonlinearity drives the energy exchange between the primary wave mode (ω p , μ p ), excited by the applied harmonic displacement, and a subharmonic wave mode (ω s , μ s ) withω s ≈ (1/2)ω p . Indeed, for sufficiently high applied harmonic amplitudes, the solution with the dominant primary wave mode response saturates and energy exchange occurs. Since the subharmonic wave mode has evanescent behavior due to local resonance, the energy becomes localized in the local resonators and does not propagate further, thus giving rise to a second attenuation zone. In nonlinear dynamics, similar phenomenon is known as a 2:1 internal resonance or autoparametric resonance mechanism. In nonlinear wave dynamics, however, to the best of our knowledge, this phenomenon has not yet been demonstrated. Within the wave framework, it requires the consideration of the convective forces which does not play a role in the case of a single resonator. Moreover, since it occurs between waves of different nature, i.e., propagating and evanescent waves, subharmonic attenuation zones in metamaterials can emerge. This novel feature of nonlinear locally resonant metamaterials has the potential to inspire new metamaterial designs as a single local resonance may yield multiple attenuation zones.