© R. Zivieri, Published by EDP Sciences, 2022

1 Introduction

After the seminal work by Veselago, who introduced the concept of negative electromagnetic media (sometimes called left-electromagnetic 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, super-resolution 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 backward-like 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 one-dimensional (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 mass-spring 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 mass-in-mass (p.m.m.) acoustic metamaterial is treated, via an effective medium approach, as an equivalent periodic mass-spring (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 (S-waves) 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 mass-spring 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 tetra-atomic 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 fourth-order 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 mass-spring (p.m.s.) and periodic mass-in-mass (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 third-order 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 Floquet-Bloch 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 fourth-order 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 κ₀ 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⁽²⁾ and by the fourth-order anharmonic constant k_e⁽⁴⁾ 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⁽²⁾ = k_e⁽²⁾ and k⁽⁴⁾ = k_e⁽⁴⁾. 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 P-mass 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−(ω/ω₀)²), namely it is frequency dependent [20,22]. For a periodic system subjected to the Floquet-Bloch periodic boundary condition, by imposing that the dispersion relation of the effective mass-spring system is identical to that of the mass-in-mass 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 ω₀= √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 ω → ω₀⁻, M_eff → −∞ (in Kg) as ω → ω₀⁺, namely M_eff is huge near resonance, M_eff vanishes for ω = ω₀ √M_t/M_e [20,21] and, in some cases, M_eff can be complex.

The fourth-order 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 fourth-order anharmonic expansion (third-order anharmonic force) with V₀ the zero-order term, $V_2=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.k^{(2)}\left({\left(\mathrm{\Delta }{s}^{j,j-1}\right)}^2+{\left(\mathrm{\Delta }{s}^{j+1,j}\right)}^2\right)$ the second-order term and $V_4=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4!$}\right.k^{(4)}\left({\left(\mathrm{\Delta }{s}^{j,j-1}\right)}^4+{\left(\mathrm{\Delta }{s}^{j+1,j}\right)}^4\right)$ the fourth-order term being Δs^j,j−1 = s^j − s^j−1, Δs^j+1,j = s^j+1 − s^j. The third-order term V₃ is set equal to zero because it would introduce a non-reciprocal 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)+2k^{(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 ⁽²⁾ (κ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 ⁽⁴⁾ (κ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 fourth-order 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₀, 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}}{1-q}\right)$ [40]. The amplitude A₀ depends on the metamaterial parameters M, k⁽²⁾, k⁽⁴⁾, 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{\sin }^2\theta }}$ and $K^{\prime }\left(m_1\right)={\int }_0^{\pi /2}\frac{d\theta }{\sqrt{1-m_1^2{\sin }^2\theta }}$ are elliptic integrals with K and iK^′ the quarter-periods of the elliptic function and m₁= 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 third-order 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-4M_{\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⁽²⁾ f ⁽²⁾ (κ L) +2 k⁽⁴⁾A²f⁽⁴⁾(κL) + M_tω₀². 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₀. 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⁽⁴⁾ = 0 N/m³ (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⁽⁴⁾ > 0 N/m³) the above definition corresponds to the frequency difference Δν_BG = ν^O (κ = 0) −ν ^A (κ = π/L). Instead, for anharmonic softening (k⁽⁴⁾ < 0 N/m³), the above definition is not anymore true because the top of the A branch is realized for κ₀ ≠π/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 (κ = κ₀). 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 fourth-order 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) = ω₀/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/3k^{(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⁽⁴⁾ > 0 rad/m) this condition is always fulfilled, while for softening (k⁽⁴⁾ <0 rad/m) the condition is fulfilled in the weak anharmonicity regime according to which |k⁽⁴⁾|A² < k⁽²⁾ and, at most, |k⁽⁴⁾|A² ≈ k⁽²⁾, namely the fourth-order 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 fourth-order 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⁽⁴⁾ = 0 N/m³ the BG reduces to the one of the harmonic regime and has the usual explicit dependence on the metamaterial parameters k⁽²⁾, 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)+2k^{\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(\sin \left(\mathrm{}\kappa L\right)-2\sin \left(2\kappa L\right)+\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 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-4M_{\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)\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)\sin \left(\kappa L\right)+\gamma \left(\sin \left(3\kappa L\right)-2\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 κ₀ = π/L. In the case of softening (γ < 0) that implies k⁽⁴⁾ < 0 N/m³, equation (11) admits a non-trivial 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 -{\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 < κ₀ < π/L, κ₀ = κ₀(γ) and with the inverse function tan⁻¹ positive. Straightforwardly, the group velocity of both the A and O modes is negative (v_g < 0 m/s) for κ₀ < κ < π/L. The solution to equation (12) is real for −1 < γ ≤−0.25. In particular, for γ = −0.25 the argument x of tan⁻¹(x) vanishes and the solution to equation (12) merges with the trivial solution to equation (11) at κ₀= π/L. For−0.25 < γ < 0 the solution to equation (12) becomes complex. The more γ increases negatively, the more κ₀ 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 κ₀. 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 fourth-order 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⁽²⁾ = 1.55× 10⁵ N/m, k_e⁽⁴⁾ = ± 2 × 10⁸ N/m³, k_i =1.08 × 10⁶ N/m,M_i = 317 Kg. The values of k_e⁽²⁾ and k_e⁽⁴⁾ 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 ⁽⁴⁾ 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⁽²⁾ = k⁽²⁾ and k_e⁽⁴⁾ = k⁽⁴⁾. 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⁽⁴⁾, 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⁽⁴⁾ >0 N/m³ and (2) M_eff ⁻ < 0 kg and k⁽⁴⁾ < 0 N/m³, while in Figure 2c κ =i Im[κ] (Re[κ] = 0 rad/m) for: (3) M_eff⁺ > 0 kg and k⁽⁴⁾ < 0 N/m³ and (4) M_eff.⁻ < 0 kg and k⁽⁴⁾ > 0 N/m³. Let us discuss case 1) and case 4) that occur for anharmonic hardening (k⁽⁴⁾ > 0 N/m³). In case 1), Im[κ] → ∞ as ω → ω₀⁻ for ω ^A(κ = π/L) < ω < ω₀, while, in case 4), Im[κ]→ ∞ as ω → ω₀⁺ for ω₀ < ω < ω^O(κ = 0). As in the harmonic regime, passing from values just below ω₀ to values just above ω₀ 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⁽⁴⁾ < 0 N/m³). In case 3), Im[κ] → ∞ as ω → ω₀⁻ for ω ^A(κ = π/L) < ω < ω₀, while, in case 2), Im[κ]→ ∞ as ω → ω₀⁺ for ω₀ < ω < ω^O(κ = 0). Passing from values just below ω₀ to values just above ω₀ 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 fourth-order 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 fourth-order 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. ω/ω₀ 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 (κ = κ₀) and at ω = ω^O (κ = κ₀), respectively with ω^A (κ = κ₀) > ω^A (κ = π/L) and ω^O(κ = κ₀) > ω^O (κ = π/L). At ω^A (κ = π/L), Re[κ] passes from π/L to 0 and remains 0 till ω =ω₀; at ω =ω₀ 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) < ω < ω₀ with Im[κ] →∞ as ω →ω₀⁻ and it monotonically decreases for ω₀ < ω ≤ ω^O (κ = 0) with Im[κ] →∞ as ω →ω₀ ⁺ and Im[κ] = 0 for ω^O (κ = 0) ≤ ω ≤ ω^O(κ = π/L).

The attenuation factor diverges for ω = ω₀, the angular resonance frequency of harmonic dynamics because, according to this model, the third-order anharmonic force acts on the external mass M_e displacement without affecting ω₀. 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⁽⁴⁾ = −2 × 10⁸ N/m³ and A = 0.02 m. The A (O) mode corresponds to in-phase (out-of-phase) displacements of M_e and M_i in the propagating regime. In the latter case the two branches exhibit a frequency maximum at κ₀L ≈ 2.25 (ν^A (κ₀ L) = 4.28 Hz and ν^O (κ₀ 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 ⁽⁴⁾ = −2 × 10⁸ N/m³) 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 fourth-order 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 κ₀L as a function of |γ | for the weak anharmonic regime where κ₀ is the wave vector expressed by the transcendental function of equation (12) at which v_g(κ₀) = 0 m/s. It can be seen that κ₀ L moves away from κ₀ L = π obtained for |γ | = 0.25, with the minimum value γ = −0.25 corresponding to A = 0.014 m with the used parameters according to which κ₀ is real (Eq. (12)), reaching κ₀ 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 < κ₀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 κ₀ 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 κ₀ < κ <  π/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) Meff + >0 Kg and k(4) > 0 N/m3 and for 2) Meff − < 0 Kg and k(4) < 0 N/m3. (c) As in panel (b) but for 3) Meff + >0 Kg and k(4) < 0 N/m3 and for 4) Meff − < 0 Kg and k(4) > 0 N/m3. 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 × 108 N/m3. The double arrow indicates the BG amplitude at the 1BZ border in the latter case. The dash-dotted 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 × 108 N/m3 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 × 108 N/m3. (d) Continuous red line with full blue circles: dimensionless wave vector κ0L 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 third-order polynomial equation F = k_e⁽²⁾ x + k_e⁽⁴⁾ x³ with k_e⁽²⁾ =1.14 ×10⁵ N/m and k_e⁽⁴⁾ = 2 ×10⁸ N/m³ 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⁽⁴⁾ > 0 N/m³.

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 κ₀ the A wave that propagates in the metamaterial could be stopped (v_g(κ₀) = 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⁽⁴⁾ < 0 N/m³ 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⁽²⁾ = x + γ₁ x³ with γ₁ = k_e⁽⁴⁾/ k_e⁽²⁾, 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, γ₁ = −1290.32/m² using k_e⁽²⁾ = 1.55 ×10⁵ N/m and k_e⁽⁴⁾ = −2 ×10⁸ N/m³. 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₀, 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 S-type 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 κ₀ and maximum frequency ν(κ₀) followed by the change of direction of its motion as soon as κ > κ₀ with ν (κ) < ν (κ₀). 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 ν(κ₀) at which the A wave v_g(κ₀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 Me in a p.m.m. system with softening anharmonic springs (ke (4) = − 2 ×108 N/m3).

4 Conclusions

In this work, it has been shown that, in a nonlinear p.m.m. seismic metamaterial modeled by means of a fourth-order 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 fourth-order 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

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) Meff + >0 Kg and k(4) > 0 N/m3 and for 2) Meff − < 0 Kg and k(4) < 0 N/m3. (c) As in panel (b) but for 3) Meff + >0 Kg and k(4) < 0 N/m3 and for 4) Meff − < 0 Kg and k(4) > 0 N/m3. 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 × 108 N/m3. The double arrow indicates the BG amplitude at the 1BZ border in the latter case. The dash-dotted 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 × 108 N/m3 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 × 108 N/m3. (d) Continuous red line with full blue circles: dimensionless wave vector κ0L vs. modulus |γ | of the anharmonic parameter.

In the text

	Fig. 4 Restoring force F vs. displacement amplitude of the external mass Me in a p.m.m. system with softening anharmonic springs (ke (4) = − 2 ×108 N/m3).
In the text