On negative effective mass and negative group velocity in anharmonic seismic metamaterials

. In this work, an anharmonic mass-in-mass system that can be employed as a nonlinear seismic metamaterial is represented as an equivalent anharmonic mass-spring 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 fourth-order 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 ﬁ rst 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.


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][3][4][5][6][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][9][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][16][17][18][19][20][21][22], phononic crystals with negative effective mass density [23] and acoustic cloaking [24][25][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][30][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 k 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.

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 fourth-order 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 k (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À(v/v 0 ) 2 ), 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]: with M t = M i + M e and v 0 = p k i /M i , the angular resonance frequency of the internal mass. In both cases [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 eff p:m:s: ¼ 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 0 the zero-order term, the second-order term the fourth-order term being Ds j,jÀ1 = s j À s jÀ1 , Ds j+1,j = s j+1 À s j . The thirdorder term V 3 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 ±ikL where t is the time, is linked to that in the jth cell via the Bloch factor e ±ikL in the regime of weak anharmonicity. The Bloch factor depends on the wave vector k having the role of a Bloch wave vector for real values of k.
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: with v = v(kL) the angular frequency (rad/s), f (2) (kL) = 2(1Àcos(kL)), 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 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 0 , and the guess solution proportional to A, it was found the relation between the two amplitudes A ¼ A 0 [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 À p K 0 /K is the nome appearing in the Lambert expansion of the exact solution, K m ð Þ ¼ ∫  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 third-order anharmonic force as a function of the wave vector k, 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: where M ± eff ¼ M ± eff kL ð Þ, the superscript + (À) refers to the sign + (À) in front of the square root and M eff According to equation (4), , 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 (kL). The frequencies of the A and O modes of the original p.m.m. system turns out to be: The values k ∊ R correspond to the A and O allowed propagating branches, while the values k ∊ 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 Dn BG = n O (k = 0) Àn A (k = p/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 k 0 ≠p/L at which the A branch attains a maximum so that, in this case, the BG would correspond to the frequency difference n O (k = 0) Àn A (k = k 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 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 k = p/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 k = p/L and does not correspond to the canonical BG definition. Here, 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. n A (k = p/L) is the frequency of the A mode at the border of the 1BZ obtained by equation (5) and expressed in explicit form as n A k ¼ ð that the frequency of the A mode at the border of the 1BZ is real for 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 fourth-order anharmonic contribution is not larger than the harmonic one. The difference between the frequency of the O mode at k = 0 and of the A mode at k = p/L turns out to be [41]: 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 (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 = ∂v/∂k. It can be derived from the dispersion of equation (3) rewriting it in the compact Within the effective medium approach, it is interesting to find the value of k such that v g (k = k 0 ) = 0 in the range of wave vectors 0 < k L < p (1BZ). First, substituting g 0 (kL) and g(kL) into v g one gets: According to equation (7), the group velocity of both the A and the O modes is indeterminate for k = 0 and assumes a finite value for 0 < kL p. After straightforward algebra, the derivative of the effective mass with respect to the wave vector can be casted in the compact form: with: with h + (kL) > 0 and finite, and h À (kL) > 0 and finite, for 0 < kL p. Therefore, accounting for equation (8) and equation (9), the group velocity expressed in equation (7) takes the final form: It can be shown that the quantity between round brackets is real and positive in the regime of weak anharmonicity and for 0 < kL p and depends on the branch considered. Hence, for both A and O modes, the condition for canceling the group velocity reduces to G 0 (kL) = 0 that can be set in the form: 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 g. The condition of weak anharmonicity realizes for À1 < g < 0 in the case of softening and for 0 < g < 1 in the case of hardening. For g = 0 the harmonic regime is recovered and equation (11) has the trivial solution k 0 = p/L. In the case of softening (g < 0) that implies k (4) < 0 N/m 3 , equation (11) admits a non-trivial solution in the interval of wave vectors 0 < k < p/L (as it is well known v g (k = p/L) = 0 and the value k = 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: with 0 < k 0 < p/L, k 0 = k 0 (g) 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 k 0 < k < p/L. The solution to equation (12) is real for À1 < g À0. 25. In particular, for g = À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 k 0 = p/L. For À0.25 < g < 0 the solution to equation (12) becomes complex. The more g increases negatively, the more k 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 k (see also the next section) and the A and O branches attain a maximum at a lower value of k 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 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.

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. v.
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 v v O (k = p/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 g = ± 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 v 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 (g = 0) to 10.04 Hz (g = À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  À < 0 kg and k (4) > 0 N/m 3 . Let us discuss case 1) and case 4) that occur for anharmonic hardening (k (4) . As in the harmonic regime, passing from values just below v 0 to values just above v 0 there is a sudden jump of the real part of the wave vector from k = p/L to k = 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 k = p/L (A mode) and upper limit for k = 0 rad/m (O mode). Let us now discuss case 2) and case 3) that occur for anharmonic softening (k (4) Passing from values just below v 0 to values just above v 0 there is a sudden jump of the real part of the wave vector from k = 0 to k = p/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 p/L in the real part of the complex wave vector at the resonance frequency.
In all cases the complex wave vectors k l with l = 1,2, …, n depicted in Figures 2b and 2c can be expressed in the complex plane using the exponential representation as k l = |k l | exp[iu l ] with u l = arctan(Im[k l ]/Re[k l ]) the angle formed by the wave vector with the real axis. In particular, in cases 1) and 2) u l ≠ p/2, while in cases 3) and 4) u l = p/2 ∀ l with the complex wave vector aligned along the imaginary axis (Re[k 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 v 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 k = Re[k] +i Im[k] vs. v/v 0 in the case of hardening anharmonicity looks very similar to that in the harmonic dynamics whose calculation was done in [22]. Re[k] shows some discrepancies in the softening anharmonic regime if compared to Re[k] 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 k which correspond to the same value of v yielding a breaking of the monotonic behavior for 0 v v A (k = p/L) for the A branch and for v O (k = 0) v v O (k = p/L) for the O branch. In these two regions the maximum value of frequency is at v = v A (k = k 0 ) and at v = v O (k = k 0 ), respectively with v A (k = k 0 ) > v A (k = p/L) and v O (k = k 0 ) > v O (k = p/L). At v A (k = p/L), Re[k] passes from p/L to 0 and remains 0 till v =v 0 ; at v =v 0 it passes from 0 to p/L and at v O (k = 0) passes again from p/L to 0. Hence, Re[k] 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 p/L of the real wave vector induced by the softening regime.
On the other hand, Im[k] 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 The attenuation factor diverges for v = v 0 , 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 v 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. In the latter case the two branches exhibit a frequency maximum at k 0 L ≈ 2.25 (n A (k 0 L) = 4.28 Hz and n O (k 0 L) = 15.1 Hz, respectively) and the range of dimensionless wave vectors at which v g < 0 m/s is 2.25 < k L < p. The "BG" amplitude at the 1BZ border is Dn 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 n O (k = 0) Àn A (k = p/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 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 n A (k = p/L) compared to the little reduction of n O (k = p/L) (Eq. (5)) indicating a maximum in the dispersion relation occurring at k < p/L. n A (k = p/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, n A (k = p/L) decreases but still the maximum of the dispersion relation and the vanishing of the group velocity is realized at k = p/L.
Both the widening of the BG at the 1BZ border and the decrease of n A (k = p/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 k 0 L as a function of |g | for the weak anharmonic regime where k 0 is the wave vector expressed by the transcendental function of equation (12) at which v g (k 0 ) = 0 m/s. It can be seen that k 0 L moves away from k 0 L = p obtained for |g | = 0.25, with the minimum value g = À0.25 corresponding to A = 0.014 m with the used parameters according to which k 0 is real (Eq. (12)), reaching k 0 L ≈ 2.0 for |g | = 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 (k) < 0 m/s is 2.0 < k 0 L < p. 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 k 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 k 0 < k < p/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.

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 (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 k 0 the A wave that propagates in the metamaterial could be stopped (v g (k 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 + g 1 x 3 with g 1 = k e / k e (2) , the restoring force attains a maximum at the displacement value x M ¼ 1 ffiffi 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 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 k 0 and maximum frequency n(k 0 ) followed by the change of direction of its motion as soon as k > k 0 with n (k) < n (k 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 n(k 0 ) at which the A wave v g (k 0 L) = 0 m/s depending on the strength of the softening anharmonic springs.

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 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. Support was received by National Group of Mathematical Physics (GNFM À INdAM).