Nonlinear energy localisation in a model of plane metamaterial

Applying the concepts of nonlinear normal modes and limiting phase trajectories introduced by Manevitch in Manevitch (Arch Appl Mech 77:301–312, 2007) to a two-dimensional mass–spring system, the authors propose a generalised method to tune a plane metamaterial and get the desirable resonant behaviour at short wavelengths. Indeed, the account of nonlinear coupling between the oscillators enables the localisation of energy leading the origin of a bandgap at short wavelengths regardless the existence of external disturbances. Moreover, further restrictions on the modes amplitude allow the observation of Fermi–Pasta–Ulam–Tsingou recurrence and super-recurrence in the two-dimensional metamaterial. These findings can open the way to further research in order to improve efficiency and performance of resonant metamaterials.

not present in natural materials, but emerge as result of a suitable theoretical-experimental design [1]. These peculiarities contribute to create an extremely fascinating research field. The smart capabilities of metamaterials can be achieved by their designable microgeometry and/or by their composition [2] and their potential applications range from wave impact mitigation [3,4] to sound and vibrations control [5][6][7][8][9][10][11]. At a very general conceptual level, it can be said that the study of wave propagation in periodic lattices, with the inherent dispersive properties, constitutes the theoretical framework which allows to give birth to the idea of metamaterial. More precisely, the exotic features, such as negative refraction [12,13] and tunnelling effect [14][15][16], characterising these new materials are generally recognised as a consequence of the existence of frequency ranges in which wave propagation is inhibited, so-called bandgaps. The frequencies and width of these bandgaps are ruled by the mismatch between the mechanical properties of the constituent phases and the lattice parameters. Indeed, in phononic crystals the two main mechanisms responsible of the formation of bandgaps rely in the Bragg interference of waves scattered by periodic inhomogeneites and hybridization of local resonances due to the resonant modes of the constituents interacting with the embedding continuum [17][18][19][20][21][22][23]. Lately, locally resonant metamaterials get attention because they are able to generate bandgaps at wavelengths much longer than the lattice size [24][25][26][27], for this reason they are good candidates for efficient energy harvesters, that is materials where energy harvesting is obtained converting the localised vibrational kinetic energy into electrical energy [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].
The mass-spring system is a well-known classical physical model which can explain the most of resonance phenomena although it has been applied mainly to 1D metamaterials [44][45][46][47][48][49][50][51][52]. Some attempts have been made to represent a 2D metamaterial into a mass-spring system, but nonlinear coupling terms between the oscillators are generally not considered: actually, the most of such research works focus on the linear behaviour [53][54][55][56][57][58].
To gain quantitative descriptions of these facts, in order to plan rational design tools, the analytical studies mostly focus on the spectral properties of the (linear) operator governing the motion equations of the system. Bandgaps dictate frequency ranges where vibrations propagation is forbidden, so enabling the localisation of vibrational kinetic energy in the form of oscillatory motion. Nevertheless, a clear explanation of the mechanism by which resonances lead to the occurrence of energy localisations seems to be a challenging issue which, in our opinion, can be addressed by considering the nonlinear coupling of a mass-spring system.
In the present work, we show, in a conservative 2D mass-spring metamaterial, the occurrence of energy localisation at short wavelengths. We note that the localisation arises only because oscillators interact (by springs) in resonance with one another and not because there is an external disturbance or host structure with which they resonate at long wavelengths [59,60].
In the light of the above, understanding the mechanism leading to the energy transfer and its localisation, which is the main goal of the present paper, is a key problem in the study of metamaterials and it is connected with the search of elementary excitations preserving their characteristic under mutual interactions. This leads to the concept of nonlinear normal modes (NNM) introduced in [61]. In order to find NNMs, we take advantage from the fact that we deal with symmetrical system whose equations are of the kind where n is the number of degrees of freedom (in our case n = 4) and the potential V is defined as where a (m) i j and b (m) represent elastic coefficients. This system admits similar NNMs as solutions and it possesses always solutions of the type x 1 = ±x m with m = 2, ..., n which correspond to in-phase and antiphase similar NNMs [61,62]. These NNMs correspond to the stationary points of the system. By Lyapunov Theorem [63], we know that in a neighbourhood of stable equilibrium points of an n degrees of freedom (DOF) system there exist at least n NNMs. This theorem was extended by Weinsten [64] and Moser [65] to systems with internal resonances. However, NNMs technique is effective only for the analysis of stationary states and cannot be applied for highly non-stationary resonant processes in weakly nonlinear oscillators. For this reason, the concept of limiting phase trajectory (LPT) was first introduced by Manevitch [66] for 1D mass-spring system with two degrees of freedom (DOF) and it has been showed that it places a key role in understanding intense resonant energy exchange for highly non-stationary resonant processes [66].
More precisely, localised excitations originate from two transitions in the phase space: the instability of bifurcating NNMs and, then, the instability of LPTs. On the one hand, NNMs individuate the stationary points and the encircling curves close enough to them. On the other hand, the LPT individuates the outer curve of a set of curves encircling the stationary points which represents a highly non-stationary resonant process [67,68]. The instability of the NNMs corresponds to a partial energy exchange between the oscillators and leads to weak energy localisation, whereas the instability of the LPT corresponds to the maximum energy transfer between the weakly interacting nonlinear oscillators and leads to strong energy localisation. When this transfer is broken the maximum possible energy is localised in one of the oscillators; in other words, the vibrations are concentrated mainly on one of the masses. As a result, strongly modulated non-stationary oscillations appear.
It is worth to highlight that the investigation of the NNMs and LPT evolution in the phase plane is applied only to 1D mass-spring system with two DOF in the Manevitch's works, and as soon as the degrees of freedom become more than two, the continuous limit (strong nonlinearity) is considered instead of the discrete limit (weak nonlinearity). This is due to the mathematical difficulties in dealing with more than two DOF and, a fortiori, with a 2D or 3D mass-spring system [69,70]. Manevitch's approach in more than two DOF has some weak points: • The investigation is conducted at long wavelengths (wavelengths larger than the oscillators interspaces) leaving aside the possibility of probing the mechanism of energy exchange between the oscillators detected at smaller wavelengths (high frequencies). • The searching for stationary points by the equations of motion is not possible in case of non-canonical transformations or if dealing with many DOF.
Here, we propose a more general method in order to find energy localisation in a 2D mass-spring system with more than two DOF. We consider weak nonlinear equations of motion in the discrete regime for a unit cell with 4 DOF generating by periodicity a plane metamaterial, we name it the Z cell. To detect the occurrence of energy localisation, our strategy is to focus on the Hamiltonian by computing its Hessian and related eigenvalues in order to classify the stationary points in the whole 4D phase space. Once done that, we take into account every 2D subspace and we see if the nature of the stationary points changes in dependence of the physical parameters. In such subspaces, we follow the evolution of NNMs and LPTs following Manevitch's method. It is worth to highlight that almost all critical points are saddle points rather than maxima or minima as the dimension of the space increases, this observation was first made in the physics of the p-spin spherical spin glasses (see, for instance, [71,72] and literature therein). This condition suggests us to focus on a 2D subspace instead of the entire phase space. Furthermore, we note that energy localisation occurs in the 2D subspace and not in the whole phase space. The weak nonlinearity approximation, and so the discrete framework, implies the existence of the slow timescale and the emergence of slowly modulated fast oscillations. For this reason, we employ the multiple timescale approximation [73] and introduce complex variables (see [74] and literature therein).
Eventually, we perform a numerical simulation adapting the Fermi-Pasta-Ulam-Tsingou (FPUT) code [75,76] to the present case with periodic initial conditions. We add a fixed boundary surrounding the 2D metamaterial and we impose restriction on the modes amplitude. In this particular case, we observe a partial exchange and weak energy localisation between two modes (i.e. the FPUT recurrence) and the maximum energy exchange between the other two other modes (i.e. the FPUT super-recurrence). Ultimately, these last modes contribute to the kinetic energy and to the nonlinear coupling, but not to the linear coupling. Since the trend of the energy does not change as the linear coupling parameter changes, we can state that the energy localisation for the 2D metamaterial with fixed boundaries and with these particular conditions occurs for any generic value of the stiffnesses. On the contrary, in the general model without fixed boundary discussed above, energy localisation takes place only for a particular value of the stiffnesses.
The numerical experiments are implemented for a 2D metamaterial composed by 4 and 16 moving oscillators and we observe energy localisation in both the cases.
Moreover, at the best of our knowledge, it is an important open question to know if the FPUT recurrence occurs in a 2D system. We have found the conditions at which FPUT recurrence occurs in the 2D metamaterial and we shall prove that the energy of the 2D metamaterial of N × N moving oscillators (with generic N ) can be meant as the sum of the energies of N FPUT chains.
The paper is organised as follows. In Sect. 2, we present the model of the 2D metamaterial and write down the equations of motion considering small oscillations of the masses around their equilibrium positions. In Sect. 3, we introduce the Z cell as a model generating a plane metamaterial and we carry out the time scale approximation. In Sect. 4, we lead the phase analysis and show that energy localisations occur via NNMs and LPTs instabilities for a particular value of the stiffnesses ratio. In Sect. 5, we generalise the FPUT numerical experiment to the case of the 2D metamaterial and we show the occurrence of energy localisations related to the observation of FPUT recurrence and FPUT super-recurrence.

Dynamics of a plane metamaterial
Let N ∈ N be given, we consider the periodic plane metamaterial in Fig. 1 with N × N number of nodes. Here, we distinguish nodes of type 0 (green dots) from the nodes of type 1 (blue dots). Any 0-node is connected with the two nearest 1-nodes by elastic rods with stiffness k 0 (dashed lines) and, in addition, any 1-node interacts by rods with stiffness k 1 (solid lines) with other four 1-nodes. All the types of nodes carry unitary masses.
We denote by t → r k νμ (t) ∈ R 2 , with k = 0, 1, the positions at time t ∈ R + of the 0-nodes and 1-nodes, respectively, and the p k their associated momenta relative to the positions (x k νμ (t), y k νμ (t)), for ν, μ = 1, ...N and k = 0, 1. The ν-index discretizes the x-axis, whereas the μ-index discretizes the y-axis as explained in details later. Let us suppose that all the nodes oscillate along the y-axis only, such that where the variablesw νμ =w νμ (t) andv νμ =v νμ (t) are the displacements from the rest positions along the y-axis of the 0-nodes and 1-nodes, respectively, whereas the constants a and b are the horizontal and vertical spacings between the nearest nodes at the rest position. In particular, the displacementsw νμ andv νμ are zero when the nodes are at the rest positions, i.e. w νμ = 0 andv νμ = 0. In this case, the coordinates of every node reduce to (x k νμ , y k νμ ) = (νa, μb) ( Fig. 1). We highlight that the mass in r 1 11 = (x 1 11 , y 1 11 ) is at the rest position (a, b) and it does not lay on the origin (0, 0) ( Fig. 1).
Taking into account that within the metamaterial the 0-nodes alternate their positions with the 1-nodes, the Hamiltonian can be written in the following way where K and U are the kinetic energy and the potential energy, respectively, and they are given by where δ ,k denotes the Kronecker delta. Moreover, we put to zero the potential energies U 0 and U 1 in any case they would contain self-interaction terms, that is the outer nodes placed at the boundaries are not fixed and some contributions to the potential should be vanished with the following criteria: Fig. 1 The metamaterial when all the nodes are in their rest positions. The coordinates r k νμ = (x k νμ , y k νμ ) are labelled by two subscripted indices: the μ-index discretizes the y-coordinate and the ν-index discretizes the x-axis, such that (x k νμ , y k νμ ) = (νa, μb), ν, μ = 1, ..., N . The constants a and b are the horizontal and vertical spacings between the two nearest nodes in their rest positions • when ν = N , for every μ = 1, ..., N , • when ν = 1, for every μ = 1, ..., N • when μ = N , for every ν = 1, ..., N , • when μ = 1, for every ν = 1, .., N We write the interaction potential between a 0-node with coordinates r 0 νμ and the two nearest nodes 1 with coordinates r 1 ν−1μ and r 1 ν+1μ , is and the interaction potential between a node 1 with coordinates r 1 νμ and the other two forwards nearest 1-node with coordinates r 1 ν+1μ−1 , r 1 ν+1μ+1 is and using the Taylor expansion, we get the perturbed potentials We compute the equations of motion from the Hamiltonian (4) with the potentials (11) and (12) after and the characteristic scale L such that the scaled solutions w → w L and v → v L are dimensionless. For any fixed ν and μ, the perturbed equations of motion arë that is where = N m,n=0 and we neglect all the terms of the order O(δ 6 ) and higher. The indices ν and μ for the displacements w are different from those of the displacements v, but we use the same indices just to not overload the formalism. Note that one has to keep in mind that some terms in Eq.16 are zero when the equations describe the dynamics of nodes at the boundaries (the equations of motion 16 do not have self-interaction terms).

Z Cell and multiscale approximation
In the following, we focus our analysis to the elementary cell, generating by periodicity the plane metamaterial, which will be called the Z cell hereinafter. The reference frame is the same pictured in Fig. 1. More in details, when every node is in its rest position, one of the 1-node has position vector r 1 11 = (a, b), the other 1-node has position vector r 1 22 = (2a, 2b), whereas the 0-nodes are in the positions pointed out by the vectors r 0 12 = (a, 2b) and r 0 21 = (2a, b). This rest configuration is outlined in Fig. 2, where the coordinates (x k νμ , y k νμ ) = (νa, μb) are specified, for ν, μ = 1, 2 and k = 0, 1, and only the position vector r 1 11 is showed to not overload the figure, but to figure out the reference frame in an easy way. The two 0-nodes have displacements w 12 and w 21 and the two 1-nodes have displacements v 11 and v 22 . Hereinafter, we drop out the second index to simplify the notation and we define To the aim of detecting slow and fast dynamical regimes, we introduce the following general scalings The system of ODEs for the Z cell is where we neglect all the terms O( n+m ).
Following the strategy introduced by Manevitch in [74], we proceed by complexifying the problem introducing the following complex variables, for ν = 1, 2, where i is the imaginary unit. Then, the ODEs (19) read at the first harmonic (we take all the terms multiplied by e i τ ) as follows where c.c. stands for "complex conjugated" terms and h.h. stands for terms that contribute to the higher harmonics (until the third harmonic) and 0.h. are the terms that contribute to the 0-frequency harmonic. The terms of order O( n+ m 2 ) contribute to the higher harmonics and 0-frequency terms. In order to ensure that the nonlinear and the linear terms do not blow up as → 0, it is necessary to have m ≥ 3 and n ≥ 1.
As emphasised in [70], a formulation of mechanical problems is not fully adequate if it is not written in appropriate scales, where "appropriate scales" are those at which non-trivial nonlinear effects are clearly manifested. Therefore, in the following discussion we find those particular values of m and n satisfying the conditions (22) such that the system (21) describes the nonlinear dynamics regime for the Z -cell in which we can clearly capture the phenomena of energy transfer and localisations. To this aim, we have to begin by checking the right timescale, so we first perturb in time the system (21) following the methodology in [73]. Let us introduce the various timescales τ = τ , = 0, 1, ... and consider the complex functions ψ ν and φ ν as functions of the variables τ 0 , τ 1 ,... . Differentiation w.r.t. τ gives the following expansion The functions ψ ν (τ ) and φ ν (τ ) are expanded as asymptotic series via the small parameter Besides the initial conditions on the functions, for any fixed ν, ψ ν0 (τ ) and φ ν0 (τ ), the functions ψ νk (τ ) and φ νk (τ ) are supposed to satisfy zero initial conditions: We highlight that the first harmonic can be written as e i τ ψ ν =ẇ ν +i w ν and e i τ φ ν =v ν +i v ν , where ψ ν and φ ν are the slow parts and e i τ is the fast part of the solutions. For this reason, the functions ψ ν (τ ) = ψ ν (τ 0 , τ 1 , τ 2 , ...) and φ ν (τ ) = φ ν (τ 0 , τ 1 , τ 2 , ...) must not depend on the natural time τ 0 at the principal order in the perturbation O( 0 ), that is the equations must be satisfied. As a consequence, we have nontrivial PDEs at the order O( ) or higher. Proof See "Appendix A".
Proof See "Appendix B".
After Proposition 3.1 and Proposition 3.2, the system (21) becomes Note that we leave out the c.c. and the 0.h. terms from the system (21). After introducing the slow times τ 1 = τ 0 , τ 2 = 2 τ 0 ,... where τ = τ 0 , in order to rule out trivial solutions and keep solutions which depend on slow times only, the solutions have to oscillate with small frequencies. In particular, these frequencies must be of the order of and cannot be finite (see "Appendix B"). The system (26), at the order 0 , reads thus, ψ ν0 and φ ν0 do not depend on τ 0 and depend on the slow times only, that is ψ ν0 = ψ ν0 (τ 1 , τ 2 , ...) and At the order , we have In order to avoid secular terms, we study the system Eventually, we look for solutions of (29) in the form of single-frequency exponentials, namely hence, we get which is the ODEs system we will focus on in the sequel of this paper.
The system (31) possesses Hamiltonian structure (see "Appendix C"), and so there exists the Hamiltonian function H such that where By introducing the total mass M of the system as it is easy to check that H and M are two invariants of the motion. Then, the unknown f ν and g ν can be written in terms of the scaled mass M = M 4 and of the frequency variables δ μ , μ = 1, .., 4, as follows Introducing the frequency gaps the Hamiltonian (34), in terms of these variables, takes the form We claim that Hamiltonian (38) exhibits the features allowing the study of energy transfer and localisation whose existence we are going to detect. This purpose will be achieved in the following section through an analysis in the four-dimensional phase space (θ, 1 , 2 , 3 ).

Phase space analysis
We aim to find conditions of existence of stationary and non-stationary localised excitations when resonant processes take place, that is we want to find for which values of the physical parameters energy localisations occur in the metamaterial. In the following, we fix the interspaces, a = b = 10 √ 2 [79], the scale, L = 1, and the intensity of the excitation, M = 10. We vary the parameterβ = 2 −2 k 1 k 0 , that is we vary the ratio between the stiffness k 1 and k 0 . The parameterβ tunes the linear coupling between the oscillators and the characteristic scale L tunes the nonlinear coupling between the oscillators.

Definition 4.1
NNMs are the stationary points of the nonlinear Hamiltonian system (31).
Therefore, in order to find and investigate the NNMs, one has to search for the stationary points of the Hamiltonian system (31). Then, we study the 2D subspace of the phase space associated to a particular couple of oscillators. We observe that there are several energy levels which encircle the NNMs within this subspace. These curves correspond to non-stationary processes and the energy that the oscillators exchange increases from the closest to the furthest curve. Among these curves, one can individuate the LPT as the outer energy level, i.e. it is the maximum energy possible that two oscillators can exchange when they interact via weak nonlinear interaction.

Definition 4.2
The LPT corresponds to highly nonstationary resonant processes, equivalently the LPT is the outer curve encircling the NNMs in the phase plane.
The LPT can transform into a curve that separates the NNMs one from another and, when it breaks, the energy exchange is not possible anymore. This kind of curve is named a separatrix and it is defined below.

Definition 4.3
The separatrix is a boundary between domains with distinct dynamical behaviour (phase curves) in a dynamical system.
First of all, we find some stationary points, i.e. the NNMs, and we compute the value ofβ at which they change their nature (for instance, they are minima and then become saddle points).
For any fixed θ ∈ R, the relations (36) describe simple oscillators, so we shall say that two oscillators are in-phase if there exist values of δ ν such that they oscillate with the same phase.
Hereinafter, we fix θ = π 4 , so that when δ 1 = δ 2 = δ 3 = δ 4 or when the phases δ j differ one from another by a value 2π n, n ∈ N, by (36), all the oscillators f 1 , f 2 , g 1 , g 2 oscillate with the same phase and, as a consequence, the entire system is in-phase. Otherwise, one can choose values of the phases δ ν such that we have at least a couple of oscillators that oscillate in anti-phase. For instance, when δ 1 = δ 2 = 0 and δ 3 = δ 4 = π , we have f 1 and f 2 , respectively g 1 and g 2 , which oscillate in-phase, but every oscillator f ν is in anti-phase with every oscillator g ν . Therefore, we can have two situations: • the case in which δ ν are equal one with another or they differ one with another only for a phase 2nπ with n ∈ N, so the entire system oscillates in-phase in a 4D phase space; • the case in which in a 2D phase subspace two oscillators are in-phase and in another 2D phase subspace two other oscillators are in anti-phase.
As a consequence, in the 4D phase space one cannot classify the NNMs as purely anti-phase, but one can individuate anti-phase NNMs in a 2D phase subspace.
In order to find the stationary points of the Hamiltonian H = H(θ, 1 , 2 , 3 ) given in (38), we solve the equations then we compute the Hessian Hess(H) and its eigenvalues to classify them. Because of the periodic structure of the phase space, we do not need to take into account all the stationary points. For θ > 0, two of the stationary points (θ, 1 , 2 , 3 ) are π 4 , π, π, π and π 4 , π, π, 0 corresponding to NNMs, respectively. For the stationary point π 4 , π, π, π we have that, for δ 1 = 2π , δ 2 = δ 3 = π and δ 4 = 0, then 1 = 2 = 3 = π , so that the oscillators w 1 and v 2 oscillate in-phase, but they are in anti-phase with the oscillators w 2 and v 1 which, in turn, oscillate in-phase one with another.
For the stationary point π 4 , π, π, 0 a possible choice is δ 1 = δ 2 = π and δ 3 = δ 4 = 0, such that the oscillators w 1 and w 2 are in-phase, whereas they are in anti-phase with the oscillators v 1 and v 2 , which in turn oscillate in-phase one with another.
The analysis via the second-derivative test shows that the eigenvalues of the Hessian matrix are − 4 5 (3+ 5 2β ), − 3 10 , − 3 10 , −β 2 in the stationary point π 4 , π, π, π ), such that they are associated to a local maxima for anyβ > 0. Then, this point is a local maximum in any 2D phase subspace. Therefore, the NNMs cannot bifurcate neither in the entire 4D space nor in any 2D subspace.
The 2D phase space (θ, 3 ) of in-phase oscillations is obtained fixing 1 = 2 = π at the stationary point We have a phase transition and bifurcation of NNMs in the 2D phase space (θ, 3 ) for β < 6 5 . Looking at the Z cell in Fig. 2, this means that the ratio of the stiffness k 1 between the two oscillators v 1 and v 2 and the stiffnesses k 0 between the two oscillators w 1 and v 1 (or w 2 and v 2 ) is very small, k 1 k 0 < 3 5 2 . Subsequently the formation of NNMs and their bifurcations, the LPT modifies: for instance, its deformation is clear forβ = 0.92, (Fig. 3c). Atβ = 0.9002, it becomes unstable and transforms in a separatrix (Fig. 3d). This separatrix separates the phase plane in two domains around the NNMs: the first includes all the circular curves which represent the periodic oscillations around the right-hand NNM, the second includes all the circular curves which represent the periodic oscillations around the left-hand NNM. This instability corresponds to a complete energy transfer between the oscillators: all the energy of an oscillator is surrendered to the other oscillator [66][67][68][69][70]74]. Forβ = 0.9, the LPT breaks (see Fig. 3e and its zoom-in Fig. 3g) and we have energy localisation forβ ≤ 0.9, that is k 1 k 0 ≤ 2 0.45, at which NNMs in the 2D phase subspace (θ, 3 ) drift apart and oscillators do not exchange energy anymore. At this point, all the energy is localised in just one of the oscillators (see Fig. 3f and its zoom-in in Fig. 3h). One notes that, the characteristic scale L is O(1), such that the dynamics takes place when the nonlinear coupling dominates the linear coupling.
By a further discussion, we can corroborate the observations of the evolution of the LPT in the 2D subspace. First of all, we can compute analytically the value ofβ at which the LPT becomes a separatrix.
Keeping in mind the Definition 4.3, we observe that the LPT can be individuated in the 2D phase plane by some initial conditions at τ 1 = 0. Indeed, θ(0) = the LPT breaks to move away the NNMs, (e) and its zoom-in (g). Finally, the last two figures are both forβ = 0. 80. In (f) it is clear that the NNMs do not exchange energy anymore, the zoom-in (h) shows that the energy is completely localised 0 implies g 1 (0) = 0 (see definition (36)) and, going backwards to the formulas (30), (24b) and (20b), we obtain v 1 = 0 andv 1 = 0 when τ = 0. Therefore, the oscillator g 1 is initially at rest, while the oscillator g 2 is excited.
Let us consider the reduced Hamiltonian in the 2D phase space Since, when τ 1 = 0 and for the initial condition θ(0) = 0, the segment (0, 3 ) must belong to the LPT, then the expression of the Hamiltonian on this segment is and, for any generic θ(0) on the LPT, the expression of the Hamiltonian H(θ, 3 ) on the LPT (at initial time τ 1 = 0) satisfies the following condition Furthermore, the LPT becomes a separatrix when (see Fig. 3d) that is satisfied whenβ = 0.9002. Later, we prove that the dynamics of the system on the LPT is described only by the evolution in time of the variable θ(τ 1 ). Indeed, by the numerical analysis and by varying the values ofβ, we see that the profile of θ(τ 1 ) changes its shape at the same value ofβ at which the LPT transforms in the 2D subspace. Indeed, the reduction to the 2D phase space leads to the following system Let us focus on the first equation of the system (44). By the condition (42), we express sin 3 as function of θ . Then, we obtain the following differential equation that is satisfied on the LPT [67] which, by separation of variables, gives The integral in the left-hand side has as solution the Elliptic function, whereas the integration of the righthand side yields ∓βτ 1 . Nevertheless, we do not compute explicitly the solution of (46), but we provide a numerical analysis of the evolution of its shape as the values ofβ decreases, as shown in Fig. 4. Following Manevitch's works (see, for instance, [67]), we know that the LPT is described by non-smooth functions. In particular, we plot a saw-tooth periodic function for β = 1.2 (Fig. 4a). Moreover, the numerical computations confirm that the LPT allows the maximum exchange of energy between the oscillators whenβ = 0.9002 (Fig. 4b), and it breaks whenβ = 0.90 (Fig. 4c), i.e. when k 0 k 1 = 2 0.45 (see Fig. 4).
Finally, the previous analysis does not need to be repeated for other characteristic scales L, due to the invariance of the Hamiltonian.

Numerical investigations on energy exchange
In this section, we show that a 2D metamaterial with fixed boundaries can be meant as the sum of Fermi-Pasta-Ulam-Tsingou (FPUT) 1D chains [75]. We investigate the occurrence of energy localisation by performing the well-known FPUT numerical experiment [76] generalised to the 2D case by the Runge-Kutta-45 algorithm (tolerance 10 −4 ). Then, similarly to [76,78], we think our problem in terms of normal modes depending on the time variable t ∈ R + .
We add fixed boundaries at μ = 0 and μ = N + 1, such that v ν0 = w ν0 = 0 and v ν N +1 = w ν N +1 = 0, with ν = 0, ..., N +1, and at ν = 0 and ν = N +1 such that v 0μ = w 0μ = 0, for μ = 0, ..., N +1. Actually, the  (7) are encountered. Indeed, the boundary conditions are different from those provided in the previous section because the outer moving oscillators interact with the fixed ones. As a consequence, the ODEs we are going to integrate numerically are the N × N equations (16), but the equations of the moving outer oscillators have self-interaction terms.
The initial conditions for the displacements v νμ and w νμ , for 1 ≤ μ ≤ N and 1 ≤ ν ≤ N , are set as follows: the initial positions are zero except those of the nodes with ν = 1 which take the expression of the initial conditions provided in the FPUT experiment (see code in [76] and the formulas (47) and (48)), whereas the initial velocities of all the nodes are zero.
Moreover, since the oscillators move along the ydirection only, the Fourier transform is performed only for the discrete variable μ. Similarly to [75], we introduce the normal modes A νk = A νk (t) and B νk = B νk (t) for the N × N moving nodes, and with frequencies The boundary conditions are and the initial conditions for the positions are set different from zero only for the ν = 1, that is whereas, for every μ, ν = 1, .., N , The numerical integration does not involve the quadratic terms in (16) because they contribute to harmonics greater than the first harmonic as explained in Sect. 5. Note that in (47) and (48) the subscripts μ and ν for the displacements w are different from those of the displacements v, indeed 0-nodes and 1-nodes alternate their positions in the metamaterial. In this respect, we can use the matrix formulation as follows where The entries of the matrix U (t) are the solutions of the nonlinear ODEs over the integration time and they are computed via the Runge-Kutta method. However, in the following discussion, we focus only on the energy contributing to the linear part of the dynamics by analogy with Fermi's research work [75].

Proposition 5.1 If the masses of the metamaterial in study oscillate only along the y-direction and the moving masses satisfy
even if ν 1 = ν 2 , then the kinetic and the potential energy contributing to the linear part of the equations of motion read i.e. the energy E is the sum of the energy of N FPUT chains [75], except the 0-nodes do not have the potential contribution and the coupling parameter encloses the geometry of the metamaterial.
Proof See "Appendix E".
We highlight that the condition (55) holds only when t > 0, instead at t = 0 the initial conditions (51) and (52) are satisfied. Nevertheless, it easy to check that the energy E defined in (56) holds even at t = 0 with the initial conditions (51) and (52). Moreover, by considering only the energy contributing the linear part of the dynamics, the conditions (55) implies that, for every fixed μ = 1, ..., N , the displacements satisfy w ν 1 μ = w ν 2 μ and v ν 1 μ = v ν 2 μ , with ν 1 = ν 2 .
We plot the time evolution of the energy (56) of all the modes in both the two cases with N = 4 and N = 16 moving oscillators. In all the plots, we put L = 1, = 0.1, a = b = 10/ √ 2 and = 10. In Fig. 5, we plot the energy E of all the 4 modes in a metamaterial composed by 4 moving oscillators at fixedβ = 0.9 and the integration time is t = 190. We set initial conditions different from zero only for the displacements u 11 (0) and u 12 (0), whereas the velocities are set to zero (i.e.u 11 (0) = 0 andu 12 = 0) such that only the energy of the mode (1, 1) has initial energy different from zero. Indeed, the displacement u 12 contributes with only the kinetic term to the energy E. After some time (t > 0), all the modes start sharing their energies, and we observe that both the energies of the modes (1, 1) and (1, 2) decrease (Fig. 5). As time goes on, the energy of the mode (2, 2) oscillates with maxima higher and higher, and simultaneously, the energy of the mode (1, 1) oscillates with maxima smaller and smaller. For instance, the mode (2, 2) gains the maximum of its energy when the energy of the mode (1, 1) has transferred about the 75% of its initial energy to the other modes. The most of the lost energy (∼ 97%) has been surrendered to the mode (2, 2). Later, the energy of the mode (1, 1) starts oscillate with maxima higher and higher, whereas the the energy of the mode (2, 2) oscillates with maxima smaller and smaller. Finally, the energy of the mode (1, 1) gains back the most of its initial energy except the 11%. This behaviour resembles the FPUT recurrence. The modes (2, 1) and (1,2) which have only kinetic terms in the energy E have similar behaviour, but they exchange the maximum energy possible. The mode (1, 2) gives away to and then gains back from the mode (2, 1) all the possible energy and this process is periodic. This resembles the FPUT super-recurrence [77]. Thus, it seems that the exchange of energy occurs in pairs between the modes. Probably, this is linked to the fact that two NNMs are stable and the other two are unstable in the phase space. Moreover, we underline that varyingβ (at least within the range 0.9 ≤β ≤ 1000), the trend of the energy E does not change, since the modes exchange the energy in the same way described above. However, asβ increases the recurrences occur in a time smaller and smaller and the energy into play is larger and larger. For instance, for β = 1000, the mode (1, 1) gains back almost its initial energy except the 0.5% just after an integration time t = 5. For the sake of completeness, we plot some modes of a metamaterial composed by 16 moving oscillators withβ = 0.9 (see Fig. 6). At the beginning of the simulation (t = 0), all the nodes on the first column of the metamatrial (ν = 1) have initial conditions for the dis-placements different from zero. After some time, the excited mode (1, 1) gives way its energy, and after a time t = 190, it gains back its energy up to around 37%. It is worth to highlight that there are some modes which are not excited at every time t ≥ 0 (for example, the modes (4, 3) and (4, 2)) and they are not plotted. Instead, we observe that the modes (4, 1) and (4,4) with initial energy zero gain their energy when some other modes are loosing their energy (see, for instance, the modes (1, 1) and (3, 1)). Even for the metamaterial with 16 moving oscillators, we observe the same trend of the energy E asβ varies (0.9 ≤β ≤ 100) and the same considerations for 4 moving oscillators hold as well.
Similarly to the outcomes obtained for the 2D toymodel via the phase space analysis in which the energy exchange takes place between two modes only, even in the numerical experiments we observe that the energy exchange may take place in pair between the modes in the metamaterial with 4 moving oscillators. However, without fixed boundaries and without further restrictions on the modes amplitude, the energy localisation occurs only for particular values of the stiffnesses. This result provides conditions for which FPUT recurrence and super-recurrence can be observed in a 2D metamaterial.
Analogously to the Fermi's work on a chain of nonlinearly interacting oscillators, we would expect equipartition of energy, that is the energy introduced in some modes should drift into the other modes. Instead, just like the FPUT experiment, we observe that, after some period of time, the initially excited modes gain back a large percentage of their energies. This is due to the fact that the energy E, which contributes to the linear dynamics of the system, is affected by the terms of the energy which contributes to the nonlinear dynamics, that is the nonlinear terms cause localisation of energy instead of equipartition of energy. Moreover, we should keep in mind that the linear dynamics can be rediscovered in the linear approximation, i.e. when the nonlinear coupling between the oscillators approaches zero. Thus, we can say that the study of the energy E allows us to know qualitatively how much the nonlinear interactions "disturb" the linear dynamics of the system.

An insight on the energy exchange in the Z cell
In the numerical experiment on the Z cell performed in Sect. 5, we set the initial conditions (51) and (52) for which all the first column (ν = 1) of the material is excited, whereas the other two oscillators are at rest. With these initial conditions, the plot in Fig. 5 shows partial energy exchange between the two modes (1, 1) and (2,2). The initial excited mode (1, 1) surrenders a certain quantity of its initial energy to the mode (2,2). This phenomenon corresponds to the existence and instability of two NNMs in the phase plane. We have showed that the instability of these NNMs leads to the formation and subsequent instability of the LPT in the 2D subspace. On the LPT, the two unstable modes exchange the maximum energy possible. This means that, at the time t = 0, the LPT in the 2D subspace can be individuated by the initial conditions for which one of the two unstable mode is completely excited and the other one is at rest (see discussion in Sect. 4). Therefore, we repeat the same experiment in Sect. 5, with the same boundary conditions (50), but now we set the initial conditions for which only the displacement u 22 (0) is different from zero, whereas all the other displacements and velocities are zero, as a consequence, at t = 0, only the mode (2, 2) is excited.
In Fig. 7, we observe that the energy of the mode (2, 2) oscillates with maxima smaller and smaller, whereas the energy of the mode (1, 1) reaches its maxima values around the half of the integration time. When the energy of the mode (1, 1) reaches its maximum, the energy of the mode (2, 2) looses around the 99.5% of its initial energy. At the end of the period of Note that the other two modes (1, 2) and (2, 1) have the same initial conditions in both the two numerical experiments and their behaviour is the same in both Figs. 5 and 7. We suggest that these two modes are associated to the two stable modes in the phase space.
The numerical experiments on the Z cell with different initial conditions suggest that its nonlinear behaviour depends on which oscillator is excited and which one is at rest at the initial time.
Finally, we aim to analyse deeper the connection between the occurrence of FPUT recurrence and superrecurrence with the instabilities of NNMs and LPTs in future research works, even for a system with number of oscillators bigger than 4.

Conclusions
A simplified model of plane metamaterial represented by a 2D mass-spring system with weakly nonlinear coupled oscillators seems to be suitable for understanding the formation of energy transfer and localisation. We have considered a toymodel composed by 4 oscillators in order to prove that weakly nonlinear coupling between them is a fundamental characteristic for internal resonant bandgap existence at short wavelengths (high frequencies). In this way, the working resonant bandgap of a metamaterial can be augmented by adding weakly nonlinear coupling between the oscillators besides the interaction with incident waves or excited host structures. The numerical investigation via a tailored FPUT code suggests that this result may be extended to a system with a finite number of oscillators (bigger than 4) and that the energy localisation may be observed for generic values of the stiffnesses by adding fixed boundaries, particular restrictions on the modes amplitude and proper initial conditions. These outcomes pave the road to other future research on resonant metamaterials. Data availability Enquiries about data availability should be directed to the authors.

Declarations
Competing interests The authors have not disclosed any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A: Proof of Proposition 3.1
Let us perturb in time the left-hand side of the ODEs (21) taking into account (23) and (24). Considering only the first harmonics, until the order O( 2 ), the system (21) reads Let k 1 = m − 3 = k 2 = n − 1, for m, n ∈ N, then conditions (22) require Let us focus on different cases.
(1) If k 1 = k 2 = 0, that is m = 3 and n = 1, we have that the nonlinear terms and the linear terms are multiplied by O( 0 ) and, after perturbing in time (57), at the principal order O( 0 ), we get that entails that ψ ν0 and φ ν0 depend on the natural time τ 0 for ν = 1, 2. Thus, the equations (25) are not satisfied. (2) If k 1 = 0 or k 2 = 0, that is m = 3, n > 1 or m > 3, n = 1, we have at least two solutions at the order O( 0 ) depending on the natural time τ 0 . Thus, the equations (25) are not satisfied. (3) If k 1 , k 2 ≥ 2, that is m ≥ 5 and n ≥ 3, the equations (25) are satisfied, but we have a system which does not describe a nonlinear dynamics at the order O( ). Indeed, in such a case we get and in order to avoid the secular terms we must have which does not describe a nonlinear dynamics.
We conclude that, in order to get a system describing a nonlinear dynamics at the O( ) order perturbation, the powers of at the right-hand sides must match those ones at the left-hand sides, so that k 1 = k 2 = 1, that is m = 4, n = 2.
Appendix B: Proof of Proposition 3.2 The system (21) can be written in a general form as Differentiating w.r.t. τ in (20a), we have then and multiplying by e i τ , it reads This last formula holds for φ ν and their conjugated formulas are and Summing up together (65) and (67) e i τ dψ ν dτ and summing up together (66) and (68) The expressions (69) and (70) are the left-hand side of the ODEs of the system 62 up to the factor e i τ . Let us focus on the right-hand side of the first ODE of the system 62. All the terms can be collected so that and using (20a), (20b) we obtain Equating the right-hand side of (69) with the right-hand side of (72), we get the first and the second ODEs of the system (19).
Similarly, but with a longer computation, the terms in the right-hand side of the second equation of the system 62 can be collected to give the right-hand side of the third and fourth ODEs of the system (19).
and the Hamiltonian is expressed in (38).