Issue 
EPJ Appl. Metamat.
Volume 9, 2022
Metamaterials for Novel Wave Phenomena in Microwaves, Optics, and Mechanics



Article Number  10  
Number of page(s)  11  
DOI  https://doi.org/10.1051/epjam/2022008  
Published online  16 June 2022 
https://doi.org/10.1051/epjam/2022008
Research Article
On negative effective mass and negative group velocity in anharmonic seismic metamaterials
Istituto Nazionale di Alta Matematica (INdAM), Roma, Italy
^{*} email: roberto.zivieri@unife.it
Received:
14
October
2021
Accepted:
25
February
2022
Published online: 16 June 2022
In this work, an anharmonic massinmass system that can be employed as a nonlinear seismic metamaterial is represented as an equivalent anharmonic massspring system via an effective medium approach. The dispersion relation and the behavior of the effective mass as a function of the angular frequency obtained in the regime of weak anharmonicity deviate from those of the corresponding linear system because of the effect of the fourthorder potential anharmonicity. In the presence of anharmonic soft springs it is found a range of wave vectors close to the Brillouin border zone at which the group velocity of the acoustic and optical modes is negative, namely it is opposite to the phase velocity, and a wider band gap at the border of the first Brillouin zone with respect to that of the linear case. Both effects can be tuned by varying the anharmonicity strength. The huge band gap amplitude together with the strong reduction of the frequency of the acoustic mode could be exploited for the design of nonlinear seismic metamaterials at the basis of composite foundations operating in the stop band frequencies.
Key words: Acoustic metamaterials / seismic metamaterials / massinmass system / effective massspring system / effective mass / band gap / group velocity / anharmonic springs / softening regime / composite foundations
© R. Zivieri, Published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
After the seminal work by Veselago, who introduced the concept of negative electromagnetic media (sometimes called leftelectromagnetic media) as the first class of metamaterials and of negative refraction [1], the search for novel properties in artificial materials that are not available in conventional materials has been very fruitful during the last two decades. This line of research has led, for example, to the discovery of the phenomenon of cloaking and its applications, e.g., to antenna framework, superresolution and hyperbolic dispersion [2–7].
This topic has been recently extended to unconventional metamaterials such as magnetic media where it has been shown that “negative” group velocity and negative dynamic permeability is typical of volume spin excitations having a backwardlike nature in ferromagnetic nanostructures [8–10]. In this respect, several theoretical and experimental studies have been devoted to the study of nonlinearity in different classes of metamaterials such as, e.g., the investigation of nonlinear effects in electromagnetic metamaterials where the above mentioned special properties are even enhanced [11].
On the other hand, there has been a growing interest in the study of the metamaterial properties of acoustic media including sonic materials showing localized sonic resonances and effective negative elastic constants within sonic frequency ranges [12,13] and strong attenuation bands in the audible frequency range [14], materials with negative mass density [15–22], phononic crystals with negative effective mass density [23] and acoustic cloaking [24–26]. Since the introduction of the concept of effective mass M_{eff} made by modifying Newton’s second law starting, e.g., from a onedimensional (1D) model where N cylindrical cavities are created inside a bar of a rigid material [20], several works on this topic have appeared in the literature. For example, Yao et al. [21] have carried out an experimental study on a periodic massspring system having a negative effective mass investigating the transmission properties, while Huang et al. [22] have theoretically demonstrated that the introduction of a negative effective mass M_{eff} is needed if a linear and periodic massinmass (p.m.m.) acoustic metamaterial is treated, via an effective medium approach, as an equivalent periodic massspring (p.m.s) system within the classical elastic continuum model. As a result of the introduction of a negative M_{eff} there is the creation of a band gap (BG) effect.
During the last years, there has been also a growing interest in the challenging field of seismic metamaterials that are a special class of acoustic metamaterials used, e.g., in composite foundations of buildings to filter seismic waves [27]. They represent an exciting field of investigation because of their compelling physical properties and for allowing several engineering applications. Most of the recent works were devoted to the implementation of mechanical metamaterials [28] and of seismic metamaterials having the form of a metasurface or metawedge acting on surface seismic waves [29–31]. However, less attention has been paid to bulk seismic waves. Measurements have shown that bulk shear waves (Swaves) exhibit nonlinear dynamical properties [32]. Therefore, it is important to describe seismic waves designing nonlinear seismic metamaterials for a full understanding of seismic dynamics and wave filtering. Fermi, Pasta and Ulam developed, for the first time, a dynamical model for a nonlinear massspring system [33]. Later, an attempt to describe nonlinear dynamics was done by Tsurui who studied wave modulation in anharmonic lattices [34]. Recently, Fang et al. [35] have investigated the effect of nonlinearity on wave propagation [36] in two types of 1D acoustic metamaterials made by diatomic and tetraatomic unit cells, respectively with special regard to the attenuation of waves, band structures and bifurcations. In very recent works nonlinear effects to study the dynamical properties of seismic metamaterials have been introduced via a fourthorder anharmonicity [37,38]. This has been accomplished by suggesting a general guess solution to the nonlinear equations of motion consisting of a rectangular bipolar pulses distribution to study the dynamical properties of periodic massspring (p.m.s.) and periodic massinmass (p.m.m.) systems in the weak anharmonic regime (anharmonic interaction less than the harmonic interaction) [37]. This analysis has been extended to the strong anharmonic regime (anharmonic strength greater than the harmonic one) [38] subjected to large input displacements. It has been found, by means of simulations, an attenuation mechanism associated to anharmonicity characterized by an exponential amplitude decay that is explainable by means of a scattering mechanism between a phonon and a spectrum of other phonons characterizing a multiphonon process [38,39]. Finally, the amplitude and the period of the rectangular pulse distribution were found by matching, for a nonlinear p.m.s. system, the guess solution with the exact solution to the thirdorder Duffing differential equation of motion, the Jacobi elliptic sine function [40].
The trend of the effective mass as a function of the wave vector and the dispersion relation in an anharmonic p.m.m. subjected to the FloquetBloch periodic boundary conditions in the regime of weak anharmonicity (anharmonic interaction not larger than the harmonic interaction) have been recently investigated via an effective model [41]. This has been accomplished by mapping the nonlinear p.m.m. system onto an effective nonlinear p.m.s. system by introducing M_{eff} in the presence of fourthorder anharmonic potential energy. It has also been theoretically shown that the amplitude of the BG for a nonlinear p.m.m. generated by the negative M_{eff} depends on the strength of the anharmonicity.
In this work, the effective model introduced for a p.m.m. system in [41] is discussed more deeply by giving more details about the calculation of M_{eff} and of the dispersion relation. It is also shown that the widening of the BG amplitude in the softening anharmonic regime is mainly related to the not monotonic behavior of the A branch. This effect, in turn, results in a negative group velocity region for a range of wave vectors close to the border of the first Brillouin zone (1BZ) where the group velocity is opposite to the phase velocity. However, as the anharmonic interaction affects the dynamics of the external mass of the p.m.m. system, the resonance frequency within the stop band does not change with respect to that of the harmonic case. It is also found that the group velocity of both the A and O modes vanishes at a κ_{0} value in the 1BZ depending on the metamaterial parameters used in the numerical calculations. Basing on these results, it is finally proposed the starting point for realizing a negative stiffness nonlinear seismic metamaterial characterized by softening anharmonic springs by showing a realistic trend of the restoring force as a function of the displacement. A nonlinear seismic metamaterial based on the above analysis could be designed and realized in the future for civil and mechanical engineering applications in composite foundations of buildings.
2 Modeling of anharmonic seismic metamaterials
In this section, the concept of effective mass M_{eff} is introduced in the anharmonic regime by mapping a p.m.m. system onto a p.m.s. system. This mapping is similar to that proposed in the harmonic regime [22]. According to the effective medium approach, M_{eff} of the p.m.s. system must produce the same motion of the external mass M_{e} of the p.m.m. system. This condition occurs by means of the absorption of the effect of the internal mass M_{i} obtained treating the acoustic system in the limit of a classical elastic solid object and considering the motion of rigid spheres. A pictorial representation of this mapping is shown in Figure 1 where a 1D p.m.m. is represented with j, j−1 and j+1 the unit cells separated by a lattice constant or periodicity L. The p.m.m. system is characterized by the elastic constant k_{e}^{(2)} and by the fourthorder anharmonic constant k_{e}^{(4)} of the external spring assuming that the nonlinearity affects directly the dynamics of the external mass, and by the elastic constant k_{i} of the internal spring. In the effective p.m.s. it is k^{(2)} = k_{e}^{(2)} and k^{(4)} = k_{e}^{(4)}. The time profile of the main propagating wave at a given Bloch wave vector κ (rad/m) in the 1BZ having the form of an elliptic sine function (see the following discussion for more details) is also displayed.
M_{eff} is also called Pmass because it has been defined starting from the momentum P of a rigid bar. For a single unit within a 1D model of N cylindrical cavities placed inside a bar of a rigid material subjected to a sinusoidal force it takes the form M_{eff} = M_{e }+ M_{i} 1/(1−(ω/ω_{0})^{2}), namely it is frequency dependent [20,22]. For a periodic system subjected to the FloquetBloch periodic boundary condition, by imposing that the dispersion relation of the effective massspring system is identical to that of the massinmass system, the effective mass takes the similar form [22]: $${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}={M}_{\mathrm{t}}+{M}_{\mathrm{i}}\frac{{\left(\frac{\omega}{{\omega}_{0}}\right)}^{2}}{1{\left(\frac{\omega}{{\omega}_{0}}\right)}^{2}}$$(1)
with M_{t} = M_{i }+ M_{e} and ω_{0 }= √k_{i}/M_{i}, the angular resonance frequency of the internal mass. In both cases for ω = 0 rad/s it is M_{eff} = M_{t}. In particular, M_{eff} → +∞ (in Kg) as ω → ω_{0}^{−}, M_{eff} → −∞ (in Kg) as ω → ω_{0}^{+}, namely M_{eff} is huge near resonance, M_{eff} vanishes for ω = ω_{0} √M_{t}/M_{e} [20,21] and, in some cases, M_{eff} can be complex.
The fourthorder anharmonic Lagrangian referred to the jth unit cell of the effective p.m.s. system shown in Figure 1 with no damping reads: $${L}_{\mathrm{p}.\mathrm{m}.\mathrm{s}.}^{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{1}{2}\left{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right{\left({{\displaystyle \dot{s}}}^{j}\right)}^{2}\left({V}_{0}+{V}_{2}+{V}_{4}\right)\text{.}$$(2)
Here, the first term on the second member is the kinetic energy E_{k} of the jth cell of the effective p.m.s. system where “ · ” denotes the time derivative and the second term between round brackets is the potential energy modeled via a fourthorder anharmonic expansion (thirdorder anharmonic force) with V_{0} the zeroorder term, ${V}_{2}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{1ex}{$2$}\right.{k}^{(2)}\left({\left(\mathrm{\Delta}{s}^{j,j1}\right)}^{2}+{\left(\mathrm{\Delta}{s}^{j+1,j}\right)}^{2}\right)$ the secondorder term and ${V}_{4}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{1ex}{$4!$}\right.{k}^{(4)}\left({\left(\mathrm{\Delta}{s}^{j,j1}\right)}^{4}+{\left(\mathrm{\Delta}{s}^{j+1,j}\right)}^{4}\right)$ the fourthorder term being Δs^{j,j−1} = s^{j} − s^{j−1}, Δs^{j+1,j} = s^{j+1} − s^{j}. The thirdorder term V_{3} is set equal to zero because it would introduce a nonreciprocal propagation. The displacement in the j±1 cell, s^{j±1} (t) = u^{j} (t) e^{±iκL} where t is the time, is linked to that in the jth cell via the Bloch factor e^{±iκL} in the regime of weak anharmonicity. The Bloch factor depends on the wave vector κ having the role of a Bloch wave vector for real values of κ.
In principle, the displacement u^{j} is a complex quantity being proportional to a complex wave amplitude U^{j}. However, in the derivation of the general dispersion relation for the p.m.s. and p.m.m. systems [37] it was taken the real part of the displacement so that u^{j} is assumed real with A the real amplitude. Since, for M_{eff} < 0, the kinetic energy of the effective p.m.s. system is negative contradicting classical mechanics laws according to which E_{k} ≥ 0, a possible definition of E_{k} is taking the effective mass in modulus as appears in equation (2). However, other definitions leaving the effective mass with its sign and adding other terms that make the kinetic energy positive would be, in principle, possible.
The nonlinear dispersion relation for the effective p.m.s. can be derived from the effective Lagrangian of equation (2) following the same steps as done in [37] for the p.m.s. yielding: $${\omega}^{2}=\frac{{k}^{(2)}{f}^{(2)}\left(\kappa L\right)+2{k}^{(4)}{A}^{2}{f}^{(4)}\left(\kappa L\right)}{{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}}$$(3)
with ω = ω(κL) the angular frequency (rad/s), f ^{(2)} (κL) = 2(1−cos(κL)), A the amplitude of the propagating wave representing the guess solution in the jth unit cell written in the form of rectangular bipolar pulses distribution and f ^{(4)} (κL) =1/3![1−cos(3κL) + 3cos(2κL)−3 cos(κL)]. Recently, by matching, in the propagating regime for a p.m.s. system, the exact solution of the cubic Duffing differential equation of motion, derived from a fourthorder potential energy, that expresses the displacement of the mass M of the jth cell, namely the Jacobian elliptic sine function solution of the cubic Duffing equation proportional to A_{0}, and the guess solution proportional to A, it was found the relation between the two amplitudes $A={A}_{0}\left(\frac{{\pi}^{2}}{2\sqrt{m}K}\frac{\sqrt{q}}{1q}\right)$ [40]. The amplitude A_{0} depends on the metamaterial parameters M, k^{(2)}, k^{(4)}, m is the argument of the Jacobian elliptic function, q = e^{− π K ′/K} is the nome appearing in the Lambert expansion of the exact solution, $K\left(m\right)={\int}_{0}^{\pi /2}\frac{d\theta}{\sqrt{1{m}^{2}{\mathrm{sin}}^{2}\theta}}$ and ${K}^{\prime}\left({m}_{1}\right)={\int}_{0}^{\pi /2}\frac{d\theta}{\sqrt{1{m}_{1}^{2}{\mathrm{sin}}^{2}\theta}}$ are elliptic integrals with K and iK^{′} the quarterperiods of the elliptic function and m_{1 }= 1−m. In the effective p.m.s. system the mass M is replaced by the effective mass with M_{eff} ≥ 0 (propagating regime).
The effective mass, in the presence of a thirdorder anharmonic force as a function of the wave vector κ, has been recently obtained in explicit form by replacing equation (3) into M_{eff} of equation (1) [41]. The effective mass rewritten in a compact form reads: $${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}=\frac{1}{2{\omega}_{0}^{2}}\left[G\left(\kappa L\right)\pm \sqrt{{\left(G\left(\kappa L\right)\right)}^{2}4{M}_{\mathrm{e}}{\omega}_{0}^{2}\left(G\left(\kappa L\right){M}_{t}{\omega}_{0}^{2}\right)}\right]$$(4)
where ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}={M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)$, the superscript + (−) refers to the sign + (−) in front of the square root and M_{eff}^{+} (M_{eff}^{−}) to the effective mass related to the A (O) branch with M_{eff}^{+} >0, M_{eff}^{−} ≥0 (M_{eff}^{+} > M_{eff}^{−}) and G(κ L) = k^{(2)} f ^{(2)} (κ L) +2 k^{(4)}A^{2}f^{(4)}(κL) + M_{t}ω_{0}^{2}. According to equation (4), M_{eff} ^{±} = M_{eff} ^{±}(κL, A), namely the effective mass depends also on the amplitude A of the propagating wave in turn function of A_{0}. From the nonlinear dispersion of equation (3), that is referred to the effective p.m.s. system onto which the p.m.m. is mapped, the nonlinear dispersion of the original p.m.m. system can be easily derived using the definition of G (κL). The frequencies of the A and O modes of the original p.m.m. system turns out to be: $${\nu}^{\mathrm{A},\mathrm{O}}=\frac{1}{2\pi}\sqrt{\frac{G\left(\kappa L\right){M}_{t}{\omega}_{0}^{2}}{{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)}}\text{.}$$(5)
The values κ ∊ R correspond to the A and O allowed propagating branches, while the values κ ∊ C refer to the stop band where the propagating wave of the equivalent p.m.s. system is spatially attenuated. For k^{(4)} = 0 N/m^{3} (harmonic regime) the effective mass independent of the wave amplitude is obtained. The BG amplitude is defined as the difference between the bottom of the O branch and the top of the A branch. In the harmonic regime and for anharmonic hardening (k^{(4)} > 0 N/m^{3}) the above definition corresponds to the frequency difference Δν_{BG} = ν^{O} (κ = 0) −ν ^{A} (κ = π/L). Instead, for anharmonic softening (k^{(4)} < 0 N/m^{3}), the above definition is not anymore true because the top of the A branch is realized for κ_{0} ≠π/L at which the A branch attains a maximum so that, in this case, the BG would correspond to the frequency difference ν^{O} (κ = 0) −ν ^{A} (κ = κ_{0}). However, here one focuses on the widening of the gap in correspondence of the A frequencies close to the border of the 1BZ due to the fourthorder periodic potential perturbation that adds to the opening effect caused by the negative effective mass in the linear case. In a periodic system, such as the p.m.m. system under study, it can be realistically assumed that this effect occurs for wave vectors lying on Bragg planes determined by the reciprocal translation vector on which Bragg reflections can take place in a way analogous to what occurs to electrons in a weak periodic potential in crystals [42]. In the special case the wave vector κ = π/L lies on the Bragg plane of interest and corresponds to the A frequency at the border of the 1BZ. Therefore, without loss of generality, the BG definition given for the harmonic and the anharmonic hardening regimes can be extended also to the anharmonic softening regime remembering that, in this latter case, it is referred to the BG opening occurring at κ = π/L and does not correspond to the canonical BG definition.
Here, ν^{O} (κ = 0) = ω_{0}/2π $\sqrt{{M}_{\mathrm{\backslash !t}}/{M}_{\mathrm{e}}}$ is the frequency of the O mode in the infinite wavelength limit corresponding to M_{eff}̄= 0 Kg and to the upper border of the BG. ν^{A} (κ = π/L) is the frequency of the A mode at the border of the 1BZ obtained by equation (5) and expressed in explicit form as ${\nu}^{A}\left(\kappa =\pi /L\right)=1/\pi \left(\sqrt{\left({k}^{(2)}+2/3{k}^{(4)}{A}^{2}\right)/{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa =\pi /L\right)}\right)$. Note that the frequency of the A mode at the border of the 1BZ is real for ${k}^{\left(2\right)}\ge \frac{2}{3}{k}^{\left(4\right)}{A}^{2}$. For hardening (k^{(4)} > 0 rad/m) this condition is always fulfilled, while for softening (k^{(4)} <0 rad/m) the condition is fulfilled in the weak anharmonicity regime according to which k^{ (4)}A^{2} < k^{ (2)} and, at most, k^{ (4)}A^{2} ≈ k^{ (2)}, namely the fourthorder anharmonic contribution is not larger than the harmonic one. The difference between the frequency of the O mode at κ = 0 and of the A mode at κ = π/L turns out to be [41]: $${\nu}^{\mathrm{O}}\left(\kappa =0\right){\nu}^{\mathrm{A}}\left(\kappa =\pi /L\right)=\frac{1}{\pi}\left(\frac{{\omega}_{0}^{\mathrm{}}}{2}\sqrt{1+\frac{{M}_{\mathrm{i}}}{{M}_{\mathrm{e}}}}\sqrt{\frac{{k}^{\left(2\right)}+\frac{2}{3}{k}^{\left(4\right)}{A}^{2}}{{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa =\frac{\pi}{L}\right)}}\right)\text{.}$$(6)
Equation (6) expresses the BG according to the above consideration. The fourthorder anharmonic BG amplitude depends on the amplitude of the propagating mode with A in turn function of the metamaterial parameters of the nonlinear seismic metamaterial [40]. For k^{(4)} = 0 N/m^{3} the BG reduces to the one of the harmonic regime and has the usual explicit dependence on the metamaterial parameters k^{(2)}, M_{i} and M_{e}.
Another key quantity that gives important information about the metamaterial properties of seismic metamaterials is the group velocity v_{g} = ∂ω/∂κ. It can be derived from the dispersion of equation (3) rewriting it in the compact form $\omega \left(\kappa L\right)=\sqrt{g\left(\kappa L\right)}$ yielding ${\mathrm{v}}_{\mathrm{g}}=\frac{1}{2}\frac{{g}^{\prime}\left(\kappa L\right)}{\sqrt{g\left(\kappa L\right)}}$ where $g\left(\kappa L\right)=\frac{G\left(\kappa L\right){M}_{t}{\omega}_{0}^{2}}{{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)}$ and ${g}^{\prime}\left(\kappa L\right)=\frac{\partial g\left(\kappa L\right)}{\partial \kappa}$ that gives ${g}^{\prime}\left(\kappa L\right)=\frac{{G}^{\prime}\left(\kappa L\right){M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right){M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{{\pm}^{\prime}}\left(\kappa L\right)(G\left(\kappa L\right){M}_{t}{\omega}_{0}^{2})}{{\left[{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)\right]}^{2}}$ with ${G}^{\prime}\left(\kappa L\right)=\frac{\partial G\left(\kappa L\right)}{\partial \kappa}$ and ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{{\pm}^{\prime}}\left(\kappa L\right)=\frac{\partial {M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)}{\partial \kappa}$ with ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)$ expressed in equation (4). In particular, ${G}^{\prime}\left(\kappa L\right)={k}^{\left(2\right)}{f}^{{\left(2\right)}^{\prime}}\left(\kappa L\right)+2{k}^{\left(4\right)}{A}^{2}{f}^{{\left(4\right)}^{\prime}}\left(\kappa L\right)$ with ${f}^{\left(2\right)\text{'}}\left(\kappa L\right)=\frac{\partial {\mathit{f}}^{\left(2\right)}\left(\kappa L\right)}{\partial \kappa}$ and ${f}^{\left(4\right)\text{'}}\left(\kappa L\right)=\frac{\partial {\mathit{f}}^{\left(4\right)}\left(\kappa L\right)}{\partial \kappa}$ being f^{ (2)′} (κL) = 2Lsin(3κL) and ${f}^{\left(4\right)\text{'}}\left(\kappa L\right)=\frac{L}{2}\left(\mathrm{sin}\left(\mathrm{}\kappa L\right)2\mathrm{sin}\left(2\kappa L\right)+\mathrm{sin}\left(3\kappa L\right)\right)$.
Within the effective medium approach, it is interesting to find the value of κ such that v_{g} (κ = κ_{0}) = 0 in the range of wave vectors 0 < κ L < π (1BZ). First, substituting g^{′}(κL) and g(κL) into v_{g} one gets: $${\mathrm{v}}_{\mathrm{g}}^{\mathrm{A},\mathrm{O}}\left(\kappa L\right)=\frac{{G}^{\prime}\left(\kappa L\right){M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right){M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm \text{'}}\left(\kappa L\right)(G\left(\kappa L\right){M}_{\mathrm{t}}{\omega}_{0}^{2})}{2{\left({M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)\right)}^{\frac{3}{2}}\sqrt{G\left(\kappa L\right){M}_{\mathrm{t}}{\omega}_{0}^{2}}}\text{.}$$(7)
According to equation (7), the group velocity of both the A and the O modes is indeterminate for κ = 0 and assumes a finite value for 0 < κL ≤ π. After straightforward algebra, the derivative of the effective mass with respect to the wave vector can be casted in the compact form: $${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm \text{'}}\left(\kappa L\right)={h}^{\pm}\left(\kappa L\right){G}^{\prime}\left(\kappa L\right)$$(8)
with: $${h}^{\pm}\left(\kappa L\right)=\frac{\mathrm{\Delta}{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)}{\sqrt{{\left(G\left(\kappa L\right)\right)}^{2}4{M}_{\mathrm{e}}{\omega}_{0}^{2}\left(G\left(\kappa L\right){M}_{\mathrm{t}}{\omega}_{0}^{2}\right)}}$$(9)
where $\mathrm{\Delta}{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa L\right)={M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa L\right){M}_{\mathrm{e}}^{\mathrm{}}>0$ and $\mathrm{\Delta}{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{}\left(\kappa L\right)={M}_{\mathrm{e}}^{\mathrm{}}{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{}\left(\kappa L\right)>0$ being ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa L\right)>{M}_{\mathrm{e}}$ with a minimum value ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa L=0\right)={M}_{t}$, ${M}_{\mathrm{e}}>{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{}\left(\kappa L\right)$ ∀ κ (0 ≤ κL ≤ π) and h^{+} (κL) (h^{−} (κL)) referred to ${M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{+}\left(\kappa L\right)\text{}\left({M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{}\left(\kappa L\right)\right)$ with h^{+} (κL) > 0 and finite, and h^{−} (κL) > 0 and finite, for 0 < κL ≤ π. Therefore, accounting for equation (8) and equation (9), the group velocity expressed in equation (7) takes the final form: $${\mathrm{v}}_{\mathrm{g}}^{\mathrm{A},\mathrm{O}}\left(\kappa L\right)=\left(\frac{{M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right){h}^{\pm}\left(\kappa L\right)\left(G\left(\kappa L\right){M}_{\mathrm{t}}{\omega}_{0}^{2}\right)}{2{\left({M}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\pm}\left(\kappa L\right)\right)}^{\frac{3}{2}}\sqrt{G\left(\kappa L\right){M}_{\mathrm{t}}{\omega}_{0}^{2}}}\right){G}^{\prime}\left(\kappa L\right)\text{.}$$(10)
It can be shown that the quantity between round brackets is real and positive in the regime of weak anharmonicity and for 0 < κL ≤ π and depends on the branch considered. Hence, for both A and O modes, the condition for canceling the group velocity reduces to G^{′} (κL) = 0 that can be set in the form: $$\left(2+\gamma \right)\mathrm{sin}\left(\kappa L\right)+\gamma \left(\mathrm{sin}\left(3\kappa L\right)2\mathrm{sin}\left(2\kappa L\right)\right)=0$$(11)
for 0 < κL ≤ π, with $\gamma =\frac{{k}^{\left(4\right)}{A}^{2}}{{k}^{\left(2\right)}}$ a dimensionless parameter that expresses the strength of the anharmonic force if compared to the harmonic one. Equation (11) is a transcendental equation depending on the parameter γ. The condition of weak anharmonicity realizes for −1 < γ < 0 in the case of softening and for 0 < γ < 1 in the case of hardening. For γ = 0 the harmonic regime is recovered and equation (11) has the trivial solution κ_{0} = π/L. In the case of softening (γ < 0) that implies k^{(4)} < 0 N/m^{3}, equation (11) admits a nontrivial solution in the interval of wave vectors 0 < κ < π/L (as it is well known v_{g} (κ = π/L) = 0 and the value κ = 0 rad/m is excluded). Therefore, the group velocity of the A mode and of the O mode vanishes for the wave vector at which the two branches attain a maximum, viz. for: $${\kappa}_{0}=\frac{1}{L}\left[\pi {\mathrm{tan}}^{1}\left(\frac{\sqrt{2\gamma \left(\left(\gamma +1\right)\sqrt{\gamma \left(\gamma 2\right)}\right)}}{\gamma +\sqrt{\gamma \left(\gamma 2\right)}}\right)\right]$$(12)with 0 < κ_{0} < π/L, κ_{0} = κ_{0}(γ) and with the inverse function tan^{−1} positive. Straightforwardly, the group velocity of both the A and O modes is negative (v_{g} < 0 m/s) for κ_{0} < κ < π/L. The solution to equation (12) is real for −1 < γ ≤−0.25. In particular, for γ = −0.25 the argument x of tan^{−1}(x) vanishes and the solution to equation (12) merges with the trivial solution to equation (11) at κ_{0 }= π/L. For−0.25 < γ < 0 the solution to equation (12) becomes complex. The more γ increases negatively, the more κ_{0} moves away from the 1BZ border. In other words, at fixed elastic stiffness constant, the more the strength of anharmonicity increases the more the interval of wave vectors widens making the group velocity negative for a larger interval of κ (see also the next section) and the A and O branches attain a maximum at a lower value of κ_{0}. This effect accompanies the other effect due to anharmonicity, that is the widening of the BG evaluated at the border of the 1BZ in case of softening anharmonicity. Note that the BG arises from the negative effective mass as occurs in acoustic metamaterials and can be related to the kinetic energy but its widening at the border of the 1BZ is related to the weak fourthorder periodic potential analogously to the BG opening for electronic bands in crystals caused by a weak periodic potential [42].
These two important effects, due to the nonlinear dynamics caused by the anharmonicity, enhance the metamaterial properties of periodic seismic materials that could be used for composite foundations of buildings.
Fig. 1 Sketch of the anharmonic p.m.m. system and of the effective anharmonic p.m.s. system. The time profile of the main propagating wave in the effective anharmonic p.m.s. system is also shown. 
3 Results and discussion
3.1 Effective mass, wave vector, dispersion relation and group velocity in anharmonic seismic metamaterials
The parameters used in the numerical calculations are the ones already employed in previous studies of nonlinear seismic metamaterial [37,38,40,41]: M_{e }= 245 kg, k_{e}^{(2)} = 1.55× 10^{5} N/m, k_{e}^{(4)} = ± 2 × 10^{8} N/m^{3}, k_{i} =1.08 × 10^{6} N/m,M_{i} = 317 Kg. The values of k_{e}^{(2)} and k_{e}^{(4)} in the hardening anharmonicity regime were determined by fitting the experimental data of a compressive test of a low Young modulus natural rubber [27,37] with a small adjustment for the value of the elastic stiffness. Indeed, in the anharmonic model, it was reasonably assumed that the nonlinear effects are mainly experienced by the external mass M_{e}. The value of k_{e} ^{(4)} for the softening anharmonic regime was chosen as opposite to that fitted from the compressive test for the hardening anharmonic regime. The comparable values of M_{e} and M_{i} were chosen from the parameters of the composite foundations but interchanged with respect to the initial calculations [27]. Within the effective medium approach the following settings were made: k_{e}^{(2)} = k^{(2)} and k_{e}^{(4)} = k^{(4)}. In the present work the amplitude A of the main propagating mode has been chosen in the range 0 < A ≤ 0.027 m that, at fixed k^{(4)}, corresponds to the regime of weak anharmonicity.
In a very recent work [41] the effective mass was computed as a function of the wave vector, for real and complex values, according to equation (4) in both the harmonic and anharmonic regimes. Here, the effective mass is calculated as a function of the angular frequency by substituting equation (4) in the nonlinear dispersion of equation (3) that is, in turn, replaced in equation (1) expressing M_{eff} vs.ω.
In Figure 2a the trend of the effective mass as a function of the angular frequency is shown both in the propagating region corresponding to a real wave vector and in the stop band for 0 ≤ ω ≤ ω^{O} (κ = π/L). The qualitative behavior of the effective mass is similar making a comparison among the numerical results obtained in the harmonic regime (dotted black and red lines) and the corresponding ones in the hardening anharmonic regime (solid red and blue lines) and softening anharmonic regime (solid green and cyan lines), respectively for A = 0.02 m corresponding to γ = ± 0.52 where the plus sign refers to the hardening regime and the minus sign to the softening regime. However, the trends in the propagating regions and in the stop band region are quite different. In the propagating regions the curves of the effective mass are almost superimposable even though with slight numerical differences reflecting the similar behavior found for the effective mass as a function of the wave vector for real values [41]. Instead, some discrepancies among the effective mass curves are present in the stop band as a consequence of the different trends obtained calculating the effective mass versus the complex wave vector [41]. The main discrepancy is that, for both softening and hardening and for the range of masses considered, the effective mass becomes huge (either positively and negatively) more rapidly with respect to that in the harmonic regime. Specifically, while the curves of the negative mass are superimposable and shifted towards higher ω if compared to the harmonic curve, the profile of the positive effective mass for softening is shifted towards lower frequencies with respect to that for hardening and becomes huge more rapidly. The shift towards lower frequencies of the positive effective mass in the stop region occurring for softening can be related to the greater “BG” amplitude in this regime mainly due to the downshift of the A mode branch at the border of the 1BZ. As shown in [41], with the same parameters used for the calculation of M_{eff}, the BG amplitude evaluated at the border of the 1BZ passes from 9.26 Hz in the linear case (γ = 0) to 10.04 Hz (γ = −0.52) increasing of 9%.
A deeper understanding of the stop band formation in the anharmonic regime can be achieved by analyzing the trend of the wave vector which in the BG region becomes complex leading to the spatial attenuation in the wave amplitude. In this region the wave vector has no more the role of a Bloch wave vector. A representation in the complex plane of the complex wave vector is shown in Figure 2b and Figure 2c for different cases. In Figure 2b κ = π/L + i Im[κ] for: (1) M_{eff} ^{+} > 0 kg and k^{(4)} >0 N/m^{3} and (2) M_{eff} ^{−} < 0 kg and k^{(4)} < 0 N/m^{3}, while in Figure 2c κ =i Im[κ] (Re[κ] = 0 rad/m) for: (3) M_{eff}^{+} > 0 kg and k^{(4)} < 0 N/m^{3} and (4) M_{eff.}^{−} < 0 kg and k^{(4)} > 0 N/m^{3}. Let us discuss case 1) and case 4) that occur for anharmonic hardening (k^{(4)} > 0 N/m^{3}). In case 1), Im[κ] → ∞ as ω → ω_{0}^{−} for ω ^{A}(κ = π/L) < ω < ω_{0}, while, in case 4), Im[κ]→ ∞ as ω → ω_{0}^{+} for ω_{0} < ω < ω^{O}(κ = 0). As in the harmonic regime, passing from values just below ω_{0} to values just above ω_{0} there is a sudden jump of the real part of the wave vector from κ = π/L to κ = 0 rad/m. This is in accordance with the fact that the BG evaluated at the border of the 1BZ has its lower limit for κ = π/L (A mode) and upper limit for κ = 0 rad/m (O mode).
Let us now discuss case 2) and case 3) that occur for anharmonic softening (k^{(4)} < 0 N/m^{3}). In case 3), Im[κ] → ∞ as ω → ω_{0}^{−} for ω ^{A}(κ = π/L) < ω < ω_{0}, while, in case 2), Im[κ]→ ∞ as ω → ω_{0}^{+} for ω_{0} < ω < ω^{O}(κ = 0). Passing from values just below ω_{0} to values just above ω_{0} there is a sudden jump of the real part of the wave vector from κ = 0 to κ = π/L marking an interchange with respect to that occurring in the anharmonic hardening and to what should be expected looking at the lower and upper limits of the BG. In other words, the softening regime causes a shift of π/L in the real part of the complex wave vector at the resonance frequency.
In all cases the complex wave vectors κ_{l} with l = 1,2, …, n depicted in Figures 2b and 2c can be expressed in the complex plane using the exponential representation as κ_{l} = κ_{l} exp[iθ_{l}] with θ_{l} = arctan(Im[κ_{l}]/Re[κ_{l}]) the angle formed by the wave vector with the real axis. In particular, in cases 1) and 2) θ_{l} ≠ π/2, while in cases 3) and 4) θ_{l} = π/2 ∀ l with the complex wave vector aligned along the imaginary axis (Re[κ_{l}] = 0 rad/m).
As a general observation the range of angular frequencies where the effective mass is negative is the same both in the harmonic and anharmonic regimes because the resonance frequency does not change in the presence of anharmonic interactions. This can be understood by taking into account that the fourthorder anharmonic contribution of the original p.m.m. system was included to describe the nonlinear dynamics of the external mass M_{e} that, in a real seismic metamaterial, was supposed to be the mass most affected by the nonlinearity. To modify the range of ω at which the effective mass is negative, the anharmonicity should be included also in the dynamics of the internal mass M_{i} leading to a change of the resonance frequency. However, in the softening anharmonic regime, the widening of the stop band at the border of the 1BZ occurs in the lower stop band region where the effective mass is positive and is due to the fourthorder periodic potential.
In order to fully understand the trend of the wave vector shown in Figures 2b and 2c it is crucial to study its behavior as a function of the angular frequency. In all regimes studied, the attenuation factor is expressed by the imaginary part of the wave vector which gives the spatial attenuation of the wave in the stop band region. The trend of κ = Re[κ] +i Im[κ] vs. ω/ω_{0} in the case of hardening anharmonicity looks very similar to that in the harmonic dynamics whose calculation was done in [22]. Re[κ] shows some discrepancies in the softening anharmonic regime if compared to Re[κ] in the case of hardening anharmonicity. The most relevant one is that, for a given range of real wave vectors for the propagating regime of the A and O modes, there are two values of κ which correspond to the same value of ω yielding a breaking of the monotonic behavior for 0 ≤ ω ≤ ω^{A}(κ = π/L) for the A branch and for ω^{O}(κ = 0) ≤ ω ≤ ω^{O}(κ = π/L) for the O branch. In these two regions the maximum value of frequency is at ω = ω^{A} (κ = κ_{0}) and at ω = ω^{O} (κ = κ_{0}), respectively with ω^{A} (κ = κ_{0}) > ω^{A} (κ = π/L) and ω^{O}(κ = κ_{0}) > ω^{O} (κ = π/L). At ω^{A} (κ = π/L), Re[κ] passes from π/L to 0 and remains 0 till ω =ω_{0}; at ω =ω_{0} it passes from 0 to π/L and at ω^{O}(κ = 0) passes again from π/L to 0. Hence, Re[κ] exhibits three discontinuity points and, in common with the harmonic and the anharmonic hardening cases, it has only the discontinuity point at the angular resonance frequency. The additional discontinuities stem from a translation of π/L of the real wave vector induced by the softening regime.
On the other hand, Im[κ] has a very similar trend to that of the harmonic and anharmonic hardening cases and is not affected by the softening. It is 0 for 0 ≤ ω ≤ ω ^{A}(κ = π/L), it monotonically increases for ω ^{A}(κ = π/L) < ω < ω_{0} with Im[κ] →∞ as ω →ω_{0}^{−} and it monotonically decreases for ω_{0} < ω ≤ ω^{O} (κ = 0) with Im[κ] →∞ as ω →ω_{0} ^{+} and Im[κ] = 0 for ω^{O} (κ = 0) ≤ ω ≤ ω^{O}(κ = π/L).
The attenuation factor diverges for ω = ω_{0}, the angular resonance frequency of harmonic dynamics because, according to this model, the thirdorder anharmonic force acts on the external mass M_{e} displacement without affecting ω_{0}. The divergence of the attenuation factor within the stop band corresponds to a huge displacement amplitude of the internal mass M_{i} of the p.m.m. system if compared to that of M_{e} with the two masses moving in opposite directions resulting in its negative momentum and in the negative effective mass. This trend is very similar to the one of the harmonic case.
Figure 3a displays the dispersion relation of the A and O branches of the p.m.m. system calculated within the effective medium according to equation (5) via equation (4) in the harmonic regime (dotted lines) and in the anharmonic softening regime (solid lines) for k^{(4)} = −2 × 10^{8} N/m^{3} and A = 0.02 m. The A (O) mode corresponds to inphase (outofphase) displacements of M_{e} and M_{i} in the propagating regime. In the latter case the two branches exhibit a frequency maximum at κ_{0}L ≈ 2.25 (ν^{A} (κ_{0} L) = 4.28 Hz and ν^{O} (κ_{0} L) = 15.1 Hz, respectively) and the range of dimensionless wave vectors at which v_{g} < 0 m/s is 2.25 < κ L < π. The “BG” amplitude at the 1BZ border is Δν_{BG} = 10.04 Hz showing an increase of about 9% with respect to that of the harmonic regime [41]. Note that in [41] the same dispersion relation was shown and compared also to that in the anharmonic hardening regime.
The frequency difference ν^{O} (κ = 0) −ν^{A} (κ = π/L) as a function of the wave amplitude computed according to equation (6) in the weak anharmonic softening regime (k ^{(4)} = −2 × 10^{8} N/m^{3}) is shown in Figure 3b. Basing on the considerations done in Section 2, this frequency difference can be regarded as a BG at the border of the 1BZ that is enhanced by the perturbation induced by the fourthorder periodic potential. The trend is an increase with increasing A: for vanishing amplitude corresponding to the harmonic regime the BG amplitude is 9.26 Hz, while for A = 0.027 m (approximately the maximum limit for weak anharmonicity) the BG amplitude is 10.95 Hz undergoing an increase of about 18% with respect to the harmonic case. In other words, as already found in general terms in [41], the more the anharmonicity is stronger the more the BG is wider. The widening of the BG reflects the huge reduction of ν^{A} (κ = π/L) compared to the little reduction of ν^{O}(κ = π/L) (Eq. (5)) indicating a maximum in the dispersion relation occurring at κ < π/L. ν ^{A} (κ = π/L) passes from 4.81 Hz in the harmonic regime (A = 0 m) to 3.12 Hz for A = 0.027 m corresponding to a frequency decrease of about 35% (see Fig. 3c). Note that, for 0 ≤ A ≤ 0.014 m, ν^{A} (κ = π/L) decreases but still the maximum of the dispersion relation and the vanishing of the group velocity is realized at κ = π/L.
Both the widening of the BG at the 1BZ border and the decrease of ν^{A}(κ = π/L) are dynamical metamaterial properties that could be exploited for the fabrication of nonlinear seismic metamaterials characterized by soft springs. Figure 3d displays the trend of κ_{0}L as a function of γ  for the weak anharmonic regime where κ_{0} is the wave vector expressed by the transcendental function of equation (12) at which v_{g}(κ_{0}) = 0 m/s. It can be seen that κ_{0} L moves away from κ_{0} L = π obtained for γ  = 0.25, with the minimum value γ = −0.25 corresponding to A = 0.014 m with the used parameters according to which κ_{0} is real (Eq. (12)), reaching κ_{0} L ≈ 2.0 for γ  = 0.9, the chosen upper limit of the weak anharmonic regime corresponding to A ≈ 0.027 m. In this limit, with the parameters used, the range of wave vectors such that v_{g}(κ) < 0 m/s is 2.0 < κ_{0}L < π. In this range both the A and the O modes of the p.m.m. system exhibit a negative group velocity that were found for other ranges of wave vectors also in the hardening regime in the presence of strong anharmonicity [38]. Interestingly, at κ_{0} the envelopes of the A and O waves stop and invert their direction so that the A and O modes assume a backward character for κ_{0} < κ < π/L with the group velocity in the opposite direction with respect to the phase velocity.
Finally, note that the numerical calculations could change using different parameters but the qualitative trend of the above described metamaterial would be qualitatively the same by varying the masses of the rigid spheres and the elastic stiffness constants used in Figure 2 and in Figure 3.
Fig. 2 (a) Effective mass as a function of the angular frequency in the propagating A and O regions and in the BG region. Dotted black lines: harmonic regime in the propagating region. Solid red lines: hardening anharmonic regime in the propagating region. Solid green lines: softening anharmonic regime in the propagating region. Dotted red lines: harmonic regime in the BG region. Solid blue lines: hardening anharmonic regime in the BG region. Solid cyan lines: softening anharmonic regime in the “BG” region. (b) Representation of the complex wave vector (κ = π/L + i Im[κ]) in the complex plane for 1) M_{eff} ^{+} >0 Kg and k^{(4)} > 0 N/m^{3} and for 2) M_{eff} ^{−} < 0 Kg and k^{(4)} < 0 N/m^{3}. (c) As in panel (b) but for 3) M_{eff} ^{+} >0 Kg and k^{(4)} < 0 N/m^{3} and for 4) M_{eff} ^{−} < 0 Kg and k^{(4)} > 0 N/m^{3}. The double arrow marks the jump of the real part of the wave vector either from κ = 0 to κ = π/L (case 3)) or from κ = π/L to κ = 0 (case 4)). 
Fig. 3 (a) Frequency dispersion relation of the A mode (black lines) and of the O mode (red lines). The Bloch wave vector characterizing the propagating regime is labeled as Re[κ]. Dotted lines: harmonic dispersion. Solid lines: nonlinear dispersion in the softening anharmonic regime for k^{(4)} = −2 × 10^{8} N/m^{3}. The double arrow indicates the BG amplitude at the 1BZ border in the latter case. The dashdotted line labels the resonance frequency ν_{0} = 9.29 Hz. (b) Continuous black line with full blue circles: difference between the frequency of the O mode at κ = 0 rad/m and the frequency of the A mode at κ = π/L amplitude as a function of the wave amplitude A in the softening regime for k^{(4)} = −2 × 10^{8} N/m^{3} corresponding to the BG at the border of the 1BZ. (c) Continuous blue line with full green circles: frequency of the A mode at the 1BZ boundary, ν_{A} (κ = π/L), as a function of the wave amplitude A in the softening regime for k^{(4)} = −2 × 10^{8} N/m^{3}. (d) Continuous red line with full blue circles: dimensionless wave vector κ_{0}L vs. modulus γ  of the anharmonic parameter. 
3.2 Applications in mechanical and civil engineering
The results obtained in the previous paragraphs could be employed for applications in real mechanical systems acting as seismic metamaterials. It has been found, by fitting the data measured by means of a compression test on the sample of the natural rubber of low Young modulus used for the realization of a composite foundation [27], that the restoring force F as a function of the displacement amplitude x (taken real) of the external mass M_{e} in a p.m.m. system is fitted by a thirdorder polynomial equation F = k_{e}^{(2)} x + k_{e}^{(4)} x^{3} with k_{e}^{(2)} =1.14 ×10^{5} N/m and k_{e}^{(4)} = 2 ×10^{8} N/m^{3} and monotonically increases with x up to the maximum displacement amplitude x ≈ 0.02 m [37]. Therefore, this system exhibited hardening anharmonicity being k_{e}^{(4)} > 0 N/m^{3}.
However, as shown in Section 3.1, to practically widen the BG amplitude at the 1BZ boundary by lowering the A mode frequency and to make the group velocity negative, it can be realized a nonlinear system characterized by softening anharmonicity. This means that, at a given κ_{0} the A wave that propagates in the metamaterial could be stopped (v_{g}(κ_{0}) = 0 m/s). Recalling the definition of stiffness as the rate of change of the force F with respect to the displacement x, many materials in mechanical systems exhibit a negative stiffness for a given range of displacement amplitudes. For those materials the force decreases as the displacement amplitude increases. This trend occurs when these materials reach the point of degradation and the restoring force can be even negative for some values of the displacement.
The negative stiffness could be modeled introducing an anharmonic softening that gives the nonlinear character to the seismic metamaterial. In this way, it is proposed a k_{e}^{(4)} < 0 N/m^{3} anharmonic constant to describe the interaction of the rubber e.g. with concrete (external mass M_{e}), materials already employed to build up the linear p.m.m. seismic metamaterial characterized by a periodicity L = 0.2 m [27]. The softening behavior could be related to the formation of creeps or to a type of material degradation that would lead to a change of the spring character from positive to negative. In turn, the rubber would wrap a steel cylindrical resonator (internal mass M_{i}) as occurred in the linear seismic metamaterial [27]. By writing down the expression F/k_{e}^{(2)} = x + γ_{1} x^{3} with γ_{1} = k_{e}^{(4)}/ k_{e}^{(2)}, the restoring force attains a maximum at the displacement value ${x}_{M}=\frac{1}{\sqrt{3}}\frac{1}{\sqrt{{\gamma}_{1}}}$. For the case under study, γ_{1} = −1290.32/m^{2} using k_{e}^{(2)} = 1.55 ×10^{5} N/m and k_{e}^{(4)} = −2 ×10^{8} N/m^{3}. The force as a function of the displacement amplitude of the external mass M_{e} is shown in Figure 4 and attains a maximum at x_{M} = 0.016 m.
Hence, the graph is subdivided into two regions: (1) a positive stiffness region for x < x_{M} and (2) a negative stiffness region for x > x_{M}. It is crucial to look for a working point lying in this latter region by imposing displacements x > 0.016 m which still lie in the regime of weak anharmonicity where the restoring force diminishes.
This result together with the realization of a negative group velocity for wave vectors above a threshold value resulting in the lowering of the A branch towards frequency well below 5 Hz could be employed for the design and realization of novel seismic metamaterials able to filter bulk S seismic waves. In this respect, measurements of transmittance defined, in a finite effective p.m.s. system of N unit cells, as x_{N}/x_{0}, where subscript 0 denotes the first cell and subscript N the last cell, could confirm the strong attenuation of the propagation in the stop band that is enhanced with respect to that occurring in the linear case due to the wider “BG”. Bulk seismic waves of Stype having frequencies lying in the “BG” at the border of the 1BZ that can be tuned by softening anharmonicity leading BG lower limit well below 5 Hz could be filtered by strongly reducing their amplitude. These effects would avoid their massive propagation through the buildings and their destructive power. Another result that might be practically exploited is the vanishing group velocity of the propagating wave within the metamaterial at a given wave vector κ_{0} and maximum frequency ν(κ_{0}) followed by the change of direction of its motion as soon as κ > κ_{0} with ν (κ) < ν (κ_{0}). Therefore, a further working point could consist of selecting a small interval of frequencies within the A branch (lower frequencies) by means of a filter able to include the frequency ν(κ_{0}) at which the A wave v_{g}(κ_{0 }L) = 0 m/s depending on the strength of the softening anharmonic springs.
Fig. 4 Restoring force F vs. displacement amplitude of the external mass M_{e} in a p.m.m. system with softening anharmonic springs (k_{e} ^{(4)} = − 2 ×10^{8} N/m^{3}). 
4 Conclusions
In this work, it has been shown that, in a nonlinear p.m.m. seismic metamaterial modeled by means of a fourthorder anharmonicity that is mapped onto a nonlinear effective p.m.s. system in the regime of weak anharmonicity, the effective mass and the wave vector as a function of the angular frequency deviate from the trend exhibited in the harmonic regime, especially in the stop band region where propagation is forbidden. The dispersion relation in the presence of softening anharmonic springs is not anymore monotonic and, as a result, the group velocity of both the A and O waves becomes negative. Interestingly, in the softening anharmonic regime, the gap between the frequency of the O mode for vanishing wave vector and that of the A mode at the border of the 1BZ, arising not only from the negative effective mass but also from the fourthorder periodic potential, that can be regarded as a BG, is wider with respect to that of the harmonic system and strongly depends on the amplitude of the propagating wave of the effective p.m.s. system.
By studying the realistic trend of the restoring force as a function of the external mass displacement in a p.m.m. system characterized by soft springs, it has been shown that it is possible to work in the softening region above a threshold displacement where the force decreases its intensity as the displacement increases. This can be accomplished by tuning the amplitude of the stop band using springs having negative stiffness and playing on the metamaterial wave amplitude. This analysis could be the starting point for the design and realization of a nonlinear seismic metamaterial to be employed in the composite foundations of buildings for civil engineering purposes. This novel realization is obtained via the enhancement of the filtering effect on the bulk seismic wave of S type resulting in a more efficient attenuation of the destructive effects of earthquakes on buildings.
The results derived in terms of an effective medium approach applied to a nonlinear periodic system can be generalized in further works by including, in the first place, the damping term in the effective Lagrangian and, secondly, by extending the analysis to the regime of strong anharmonicity.
Acknowledgments
Support was received by National Group of Mathematical Physics (GNFM – INdAM).
References
 V.G. Veselago, The electrodynamics of substances with simultaneously negative values of ε and µ, Sov. Phys. Usp. 10, 509 (1968) [CrossRef] [Google Scholar]
 U. Leonhardt, D. Smith, Focus on cloaking and transformation optics, New J. Phys. 10, 115019 (2008) [CrossRef] [Google Scholar]
 J.B. Pendry, Negative refraction makes a perfect lens, Phys. Rev. Lett. 85, 3966 (2000) [CrossRef] [Google Scholar]
 J.B. Pendry, D. Schurig, D.R. Smith, Controlling electromagnetic fields, Science 312, 1780 (2006) [Google Scholar]
 J. Valentine, S. Zhang, T. Zentgraf, E. UlinAvila, D.A. Genov, G. Bartal, X. Zhang, Threedimensional optical metamaterial with a negative refractive index, Nature 455, 376 (2008) [CrossRef] [Google Scholar]
 V. Drachev, V.A. Podolskiy, A.V. Kildishev, Hyperbolic metamaterials: new physics behind a classical problem, Opt. Express 21, 15048 (2013) [CrossRef] [Google Scholar]
 S. Vellucci, A. Monti, M. Barbuto, A. Toscano, F. Bilotti, Progress and perspective on advanced cloacking metasurfaces, EPJ Appl. Metamat. 8, 7 (2021) [CrossRef] [EDP Sciences] [Google Scholar]
 R. Zivieri, Magnetic matter spin waves with “negative” group velocity, in 2015 9th International Congress on Advanced Electromagnetic Materials in Microwaves and Optics (Metamaterials) (2015), pp. 529–531 [CrossRef] [Google Scholar]
 R. Zivieri, Dynamic negative permeability in a lossless ferromagnetic medium, in 2015 9th International Congress on Advanced Electromagnetic Materials in Microwaves and Optics (Metamaterials) (2015), pp. 532–534 [CrossRef] [Google Scholar]
 R. Zivieri, Dynamic permeability in a dissipative ferromagnetic medium, in 2016 10th International Congress on Advanced Electromagnetic Materials in Microwaves and Optics (Metamaterials) (2016), pp. 427–429 [CrossRef] [Google Scholar]
 M. Lapine, I.V. Shadrivov, Yu S. Kishvar, Colloquium: Nonlinear metamaterials, Rev. Mod. Phys. 86, 1093 (2014) [CrossRef] [Google Scholar]
 Z. Liu, X. Zhang, Y. Mao, Y.Y. Zhu, Z. Yang, C.T. Chan, P. Sheng, Locally resonant sonic materials, Science 289, 1734 (2000) [CrossRef] [Google Scholar]
 P. Sheng, X.X. Zhang, Z. Liu, C.T. Chan, Locally resonant sonic materials, Physica B 338, 201 (2003) [CrossRef] [Google Scholar]
 M. Hirsekorn, Smallsize sonic crystals with strong attenuation bands in the audible frequency range, Appl. Phys. Lett. 84, 3364 (2004) [CrossRef] [Google Scholar]
 J. Li, C.T. Chan, Doublenegative acoustic metamaterial, Phys. Rev. E 70, 055602 (2004) [CrossRef] [Google Scholar]
 N. Fang, D. Xi, J. Xu, M. Ambati, W. Srituravanich, C. Sun, X. Zhang, Ultrasonic metamaterials with negative modulus, Nat. Mater. 5, 452 (2006) [CrossRef] [Google Scholar]
 C.T. Chan, J. Li, K.H. Fung, On extending the concept of double negativity to acoustic waves, J. Zhejiang Univ. Sci. A 7, 24 (2006) [CrossRef] [Google Scholar]
 Y. Ding, Z. Liu, C. Qiu, J. Shi, Metamaterial with simultaneously negative bulk modulus and mass density, Phys. Rev. Lett. 99, 093904 (2007) [CrossRef] [Google Scholar]
 Y. Cheng, J.Y. Xu, X.J. Liu, Onedimensional structured ultrasonic metamaterials with simultaneously negative dynamic density and modulus, Phys. Rev. B 77, 045134 (2008) [CrossRef] [Google Scholar]
 G.W. Milton, J.R. Willis, On modifications of Newton’s second law and linear continuum elastodynamics, Proc. R. Soc. A 463, 855 (2007) [CrossRef] [Google Scholar]
 S. Yao, X. Zhou, G. Hu, Experimental study on negative effective mass in a 1D massspring system, New. J. Phys. 10, 043020 (2008) [CrossRef] [Google Scholar]
 H.H. Huang, C.T. Sun, G.L. Huang, On the negative effective mass density in acoustic metamaterials, Int. J. Eng. Sci. 47, 610 (2009) [CrossRef] [Google Scholar]
 Z.Y. Liu, C.T. Chan, P. Sheng, Analytic model of phononic crystals with local resonances, Phys. Rev. B 71, 014103 (2005) [CrossRef] [Google Scholar]
 S.A. Cummer, D. Schurig, One path to acoustic cloaking, New J. Phys. 9, 45 (2007) [CrossRef] [Google Scholar]
 L.W. Cai, S. SánchezDehesa, Acoustic cloaking in two dimensions: a feasible approach, New J. Phys. 9, 450 (2007) [CrossRef] [Google Scholar]
 H.Y. Chen, C.T. Chan, Acoustic cloaking in three dimensions using acoustic metamaterials, Appl. Phys. Lett. 91, 183518 (2007) [CrossRef] [Google Scholar]
 O. Casablanca, G. Ventura, F. Garescì, B. Azzerboni, B. Chiaia, M. Chiappini, G. Finocchio, Seismic isolation of buildings using composite foundations based on metamaterials, J. Appl. Phys. 123, 174903 (2018) [CrossRef] [Google Scholar]
 M. Miniaci, A. Krushynska, F. Bosia, N.M. Pugno, Large scale mechanical metamaterials as seismic shields, New J. Phys. 18, 083041 (2016) [CrossRef] [Google Scholar]
 A. Palermo, S. Krödel, K.H. Matlack, R. Zaccherini, V.K. Dertimanis, E.N. Chatzi, A. Marzani, C. Daraio, Hybridization of guided surface acoustic modes in unconsolidated granular media by a resonant metasurface, Phys. Rev. Appl. 9, 54026 (2018) [CrossRef] [Google Scholar]
 A. Colombi, D. Colquitt, P. Roux, S. Guenneau, R.V. Craster, A seismic metamaterial: the resonant metawedge, Sci. Rep. 6, 27717 (2016) [CrossRef] [Google Scholar]
 S. Brûlé, E.H. Javelaud, S. Enoch, S. Guenneau, Experiments on seismic metamaterials: molding surface waves, Phys. Rev. Lett. 112, 133901 (2014) [CrossRef] [Google Scholar]
 P.A. Johnson, Nonlinear acoustic/seismic waves in earthquake processes, in AIP Conf. Proc. 1474 (AIP Press, Tokyo, Japan, 2012), pp. 39–46 [Google Scholar]
 E. Fermi, J. Pasta, S. Ulam, Studies of Nonlinear Problems, Collected Papers II, University of Chicago Press, Chicago, IL, pp. 977–88 Document LA1940 (1955) [Google Scholar]
 A. Tsurui, Wave modulations in anharmonic lattices, Progress Theor. Phys. 48, 1196 (1972) [CrossRef] [Google Scholar]
 X. Fang, J. Wen, B. Bonello, J. Yin, D. Yu, Wave propagation in onedimensional nonlinear acoustic metamaterials, New J. Phys. 19, 053007 (2017) [CrossRef] [Google Scholar]
 L. Brillouin, Wave Propagation in Periodic Structures (Dover, New York, NY), 1953 [Google Scholar]
 R. Zivieri, F. Garescì, B. Azzerboni, M. Chiappini, G. Finocchio, Nonlinear dispersion relation in anharmonic periodic massspring and massinmass systems, J. Sound Vib. 462, 114729 (2019) [Google Scholar]
 S. Fiore, G, Finocchio, R. Zivieri, M. Chiappini, F. Garescì, Wave amplitude decay driven by anharmonic potential in nonlinear massinmass systems, Appl. Phys. Lett. 117, 124101 (2020) [CrossRef] [Google Scholar]
 R. Zivieri, G. Santoro, V. Bortolani, Multiphonon effects in the onephonon cross section of Al, Phys. Rev. B 58, 5429 (1998) [CrossRef] [Google Scholar]
 R. Zivieri, Dynamical properties of a periodic massspring nonlinear seismic metamaterial, in 2020 Fourteenth International Congress on Artificial Materials for Novel Wave Phenomena (Metamaterials) (2020), pp. 012–014 [CrossRef] [Google Scholar]
 R. Zivieri, Negative effective mass and stop band of nonlinear periodic seismic metamaterials, in 2021 Fifteenth International Congress on Artificial Materials for Novel Wave Phenomena (Metamaterials) (2021), pp. 481–483 [CrossRef] [Google Scholar]
 N.W. Ashcroft, N.D. Mermin, Solid State Physics (HoltSaunders, 1976) [Google Scholar]
Cite this article as: Roberto Zivieri, On negative effective mass and negative group velocity in anharmonic seismic metamaterials, EPJ Appl. Metamat. 9, 10 (2022)
All Figures
Fig. 1 Sketch of the anharmonic p.m.m. system and of the effective anharmonic p.m.s. system. The time profile of the main propagating wave in the effective anharmonic p.m.s. system is also shown. 

In the text 
Fig. 2 (a) Effective mass as a function of the angular frequency in the propagating A and O regions and in the BG region. Dotted black lines: harmonic regime in the propagating region. Solid red lines: hardening anharmonic regime in the propagating region. Solid green lines: softening anharmonic regime in the propagating region. Dotted red lines: harmonic regime in the BG region. Solid blue lines: hardening anharmonic regime in the BG region. Solid cyan lines: softening anharmonic regime in the “BG” region. (b) Representation of the complex wave vector (κ = π/L + i Im[κ]) in the complex plane for 1) M_{eff} ^{+} >0 Kg and k^{(4)} > 0 N/m^{3} and for 2) M_{eff} ^{−} < 0 Kg and k^{(4)} < 0 N/m^{3}. (c) As in panel (b) but for 3) M_{eff} ^{+} >0 Kg and k^{(4)} < 0 N/m^{3} and for 4) M_{eff} ^{−} < 0 Kg and k^{(4)} > 0 N/m^{3}. The double arrow marks the jump of the real part of the wave vector either from κ = 0 to κ = π/L (case 3)) or from κ = π/L to κ = 0 (case 4)). 

In the text 
Fig. 3 (a) Frequency dispersion relation of the A mode (black lines) and of the O mode (red lines). The Bloch wave vector characterizing the propagating regime is labeled as Re[κ]. Dotted lines: harmonic dispersion. Solid lines: nonlinear dispersion in the softening anharmonic regime for k^{(4)} = −2 × 10^{8} N/m^{3}. The double arrow indicates the BG amplitude at the 1BZ border in the latter case. The dashdotted line labels the resonance frequency ν_{0} = 9.29 Hz. (b) Continuous black line with full blue circles: difference between the frequency of the O mode at κ = 0 rad/m and the frequency of the A mode at κ = π/L amplitude as a function of the wave amplitude A in the softening regime for k^{(4)} = −2 × 10^{8} N/m^{3} corresponding to the BG at the border of the 1BZ. (c) Continuous blue line with full green circles: frequency of the A mode at the 1BZ boundary, ν_{A} (κ = π/L), as a function of the wave amplitude A in the softening regime for k^{(4)} = −2 × 10^{8} N/m^{3}. (d) Continuous red line with full blue circles: dimensionless wave vector κ_{0}L vs. modulus γ  of the anharmonic parameter. 

In the text 
Fig. 4 Restoring force F vs. displacement amplitude of the external mass M_{e} in a p.m.m. system with softening anharmonic springs (k_{e} ^{(4)} = − 2 ×10^{8} N/m^{3}). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.