Flat bands and compactons in mechanical lattices

TL;DR

This paper explores flat bands and compact localized solutions in mechanical lattices, analyzing their stability and the effects of nonlinearity, with a focus on impact constraints and anharmonic oscillators.

Contribution

It introduces new mechanical lattice models with flat bands, enabling comprehensive stability analysis of compactons under impact and smooth nonlinearities.

Findings

Existence of flat bands and compactons in mechanical lattices.

Identification of two localization types: complete and exponential.

Two mechanisms for stability loss: internal and global resonance.

Abstract

Local configurational symmetry in lattice structures may give rise to stationary, compact solutions, even in the absence of disorder and nonlinearity. These compact solutions are related to the existence of flat dispersion curves (bands). Nonlinearity can destabilize such compactons. One common flat-band-generating system is the 1D cross-stitch model, in which compactons were shown to exist for the photonic lattice with Kerr nonlinearity. The compactons exist there already in the linear regime and are not generally destructed by that nonlinearity. Smooth nonlinearity of this kind does not allow doing a complete stability analysis for this chain. We consider a discrete mechanical system with flat dispersion bands, in which the nonlinearity exists due to impact constraints. There, one can use the concept of the saltation matrix for the analytic construction of the monodromy matrix. Besides, we consider a smooth nonlinear lattice with linearly connected massless boxes, each containing two symmetric anharmonic oscillators. In this model, the flat bands and discrete compactons also readily emerge. This system also permits performing comprehensive stability analysis, at least in the anti-continuum limit, due to the reduced number of degrees-of-freedom. In both systems, there exist two types of localization. The first one is the complete localization, and the second one is the more common exponential localization. The latter type is associated with discrete breathers (DBs). Two principal mechanisms for the loss of stability are revealed. The first one is the possible internal instability of the symmetric and/or antisymmetric solution in the individual unit cell of the chain. One can interpret this instability pattern as internal resonance between the compacton and the DB. The other mechanism is global instability related to resonance of the stationary solution with the propagation frequencies.

Equations

$$\begin{array}{r}\forall n\in \mathbb{N}\le N_s:m{\ddot{x}}_n+k_{\text{shear}}(x_n-y_n)\\ +k(4x_n-x_{n-1}-x_{n+1}-y_{n-1}-y_{n+1})=0,\\ m{\ddot{y}}_n+k_{\text{shear}}(y_n-x_n)+k(4y_n-x_{n-1}-\\ x_{n+1}-y_{n-1}-y_{n+1})=0,\forall |x_n|,|y_n|<\Delta ;\\ \frac{{\dot{y}}_n{|}_{t{|}_{|y_n|\to \Delta }+\delta }}{{\dot{y}}_n{|}_{t{|}_{|y_n|\to \Delta }-\delta }}\underset{\delta \to 0}{\to }\frac{{\dot{x}}_n{|}_{t{|}_{|x_n|\to \Delta }+\delta }}{{\dot{x}}_n{|}_{t{|}_{|x_n|\to \Delta }-\delta }}\underset{\delta \to 0}{\to }-1\end{array}$$

$$\begin{array}{r}m{\ddot{v}}_n+(4k+2k_{\text{shear}})v_n=0\end{array}$$

$$v_n\triangleq x_n-y_n,\forall n\in \mathbb{N}\le N_s$$

$${\omega }_{\text{FB}}\triangleq \sqrt{\frac{4k+2k_s}{m}}$$

$$\begin{array}{r}m{\ddot{u}}_n+2k(2u_n-u_{n-1}-u_{n+1})=0\end{array}$$

$$\begin{array}{r}{\omega }_1(q)=2\sqrt{\frac{k}{m}}(1-\cos q{)}^{1/2}\end{array}$$

$${\text{x}}^{\prime }=\text{A}\text{x}$$

$$\text{A}\triangleq \left[\begin{array}{cc}\text{0}& \text{I}\\ -\hat{\text{A}}& \text{0}\end{array}\right]$$

$$\begin{array}{r}{\hat{A}}_{ij}=(4+\gamma ){\delta }_{ij}-\gamma {\delta }_{i,j+1}⟨N-1-j{⟩}^0\\ -\gamma {\delta }_{i+1,j}⟨N-1-i{⟩}^0-{\delta }_{i,j+2}⟨N-2-j{⟩}^0\\ -{\delta }_{i+2,j}⟨N-2-i{⟩}^0-{\delta }_{i,j+3}⟨N-3-j{⟩}^0\\ -{\delta }_{i+3,j}⟨N-3-i{⟩}^0\end{array}$$

$$⟨x{⟩}^0\triangleq \{\begin{array}{l}1,\forall x\ge 0\\ 0,\forall x<0\end{array}$$

$$\text{x}(t_i^{+})={\text{C}}_j\text{x}(t_i^{-})$$

$${\text{C}}_j\triangleq \left[\begin{array}{cc}\text{I}& \text{0}\\ \text{0}& {\hat{\text{C}}}^{(j)}\end{array}\right],{\hat{C}}_{ik}^{(j)}\triangleq {\delta }_{ik}-2{\delta }_{ik}{\delta }_{ij}$$

$$\text{M}=(\text{L}\text{S}\text{L}{)}^2,\text{L}=\exp \left(\frac{\pi }{2\hat{\omega }}\text{A}\right),\text{S}={\text{S}}_2{\text{S}}_1$$

$${\text{S}}_j={\text{C}}_j+(\text{A}{\text{C}}_j-{\text{C}}_j\text{A}){\text{D}}_j,j=1,2$$

$${\text{D}}_1\triangleq \frac{\text{L}{\text{x}}^0[{\text{n}}^{(1)}{]}^{\top }}{[{\text{n}}^{(1)}{]}^{\top }\text{A}\text{L}{\text{x}}^0},{\text{D}}_2\triangleq \frac{{\text{C}}_1\text{L}{\text{x}}^0[{\text{n}}^{(2)}{]}^{\top }}{[{\text{n}}^{(2)}{]}^{\top }\text{A}{\text{C}}_1\text{L}{\text{x}}^0}$$

$$x_i^0=v_i^0⟨i-N_s-1/2{⟩}^0$$

$$v_i(t)=2{\bar{v}}_1(t){\delta }_{1i}$$

$$\begin{array}{r}{\bar{v}}_1(t)=\frac{\sqrt{\kappa }}{\tan \left(\frac{\pi \sqrt{\kappa }}{2}\right)}\{\left|\pi -2\pi {\left\{\frac{\omega t}{2\pi }\right\}}_{\text{frac}}\right|-\frac{\pi }{2}\\ +\frac{4\kappa }{\pi }\left[\Xi \left(\chi ,\sqrt{\kappa }\right)-\frac{1}{16}\Xi \left(2\chi ,\frac{\sqrt{\kappa }}{2}\right)\right]\}\end{array}$$

$$\begin{array}{r}\Xi (\zeta ,a)\triangleq \frac{6{\pi }^2\{\frac{\zeta }{2\pi }{\}}_{\text{frac}}\left(1-\{\frac{\zeta }{2\pi }{\}}_{\text{frac}}\right)-{\pi }^2}{6a^2}\\ -\frac{\pi \cos \left[\pi a\left(1-2\{\frac{\zeta }{2\pi }{\}}_{\text{frac}}\right)\right]}{2a^3\sin (\pi a)}\end{array}$$

$$|\lambda {|}_{\text{max}}\underset{\gamma \to 0}{\to }1+4\sqrt{\frac{2}{3}}{q}_0\sqrt{3-4q_0^2}\gamma +\mathcal{O}\left({\gamma }^2\right)$$

$$\hat{q}(\gamma )\triangleq \frac{\pi }{2\hat{\omega }}\underset{\gamma \ll 1}{\to }\frac{\pi {\omega }_{\text{FB}}}{4\omega }\to q_0{\gamma }^{1/2}+\mathcal{O}(\gamma )$$

$$q_0^{\text{cr}}=\frac{\sqrt{3}}{2}$$

$$\frac{{\omega }_{\text{cr}}}{{\omega }_{\text{FB}}}\underset{\gamma \to 0}{\to }\frac{\pi }{2\sqrt{3}}{\gamma }^{-1/2}$$

$$u_n(t)=2\sum _{l=0}^{\infty }{\hat{V}}_{nl}\cos [(2l+1)\omega t]$$

$$\begin{array}{r}{\hat{V}}_{1,l}\triangleq ({\bar{\alpha }}_l{\hat{\alpha }}_l^2{\hat{W}}_{1l}+{\hat{\alpha }}_l{\bar{\alpha }}_l^2{\hat{W}}_{2l})p\end{array}$$

$$\begin{array}{r}{\hat{V}}_{n\ne 1,l}\triangleq ({\hat{\alpha }}_l^n{\hat{W}}_{1l}+{\bar{\alpha }}_l^n{\hat{W}}_{2l})p\end{array}$$

$$\begin{array}{r}p\triangleq \Delta {\left[\sum _{l=0}^{\infty }({\bar{\alpha }}_l{\hat{\alpha }}_l^2{\hat{W}}_{1l}+{\hat{\alpha }}_l{\bar{\alpha }}_l^2{\hat{W}}_{2l})\right]}^{-1}\end{array}$$

$$\begin{array}{r}{A}_l\triangleq \frac{{\hat{\alpha }}_l^3}{{\bar{\alpha }}_l^3}\frac{2{\alpha }_l{\hat{\alpha }}_l^{N_s-3}-{\hat{\alpha }}_l^{N_s-4}-{\bar{\alpha }}_l/{\hat{\alpha }}_l}{2{\alpha }_l{\bar{\alpha }}_l^{N_s-3}-{\bar{\alpha }}_l^{N_s-4}-{\hat{\alpha }}_l/{\bar{\alpha }}_l}\end{array}$$

$$\begin{array}{r}{B}_l\triangleq A_l\left[\left(4{\alpha }_l^2-2\frac{{\alpha }_l}{{\bar{\alpha }}_l}-1\right){\bar{\alpha }}_l^{N_s}-{\bar{\alpha }}_l^2\right]\\ -\left[\left(4{\alpha }_l^2-2\frac{{\alpha }_l}{{\hat{\alpha }}_l}-1\right){\hat{\alpha }}_l^{N_s}-{\hat{\alpha }}_l^2\right]\end{array}$$

$$\begin{array}{r}{\hat{W}}_{1l}\triangleq \frac{2}{\pi }{B}_l^{-1},{\hat{W}}_{2l}\triangleq -A_l{\hat{W}}_{1l}\end{array}$$

Full text

††journal: Physical Review E

1 Introduction

The aim of the present research is to investigate the nonlinear dynamical phenomena related to the presence of flat bands (FBs) in the oscillatory spectra of mechanical and physical systems with special structural symmetry. Such structures are known to possess a number of unique and appealing physical and dynamical properties – in particular, completely localized, or compact, solutions as well as “hidden modes”, which cannot be excited externally in the linear setting. Flat bands in physical systems have been observed and broadly discussed in the recent literature. The effect of interactions and disorder on wave transport in periodic potentials, such as electrons in crystals, is strongly amplified if the bandwidth (kinetic energy) is small. A particularly interesting situation arises when some of the dispersion bands become strictly flat with macroscopically degenerate eigenstates. In this limit, any relevant perturbation will lift the degeneracy and determine the emerging highly correlated and nontrivial eigenstates. In these systems it may be possible to control the interactions, the potential energy associated with which successfully competes with kinetic energy, such that interesting wave transport phenomena are produced. Engineering FB lattice models have been extended to three-dimensional (3D), two-dimensional (2D), and even one dimensional (1D) settings. A number of FB construction pathways using graph theory were suggested in the literature. One example is the case where compact states that are fully localized on several lattice sites are used.

The origin of the compact flat band states (FBSs) is the destructive interference effectively decoupling the FBSs from the rest of the lattice, similar to the one in the known case where antisymmetric bound states embedded in and decoupled from the continuum are observed, as in [1], where the context is Fano lattices, the mechanism of resonance in which the notion of flat bands helps to explain. Fano lattices appear, for example, in the context of plasmonics, which is relevant for laser technology, and also in the broader context of meta-materials [2]. An interesting feature of the flat bands is that they are (still hypothetically) characterized by fractional power-law dependence between correlation length and level of disorder [3]. Furthermore, flat bands have extremely singular densities of states. Therefore, such systems are prone to form correlated states with broken symmetries [4]. To elaborate on the aforementioned feature of destructive interference, it may be argued that flat (dispersion) bands associated with the propagation of waves in latticed materials occur when perfect destructive interference allows for compact localized eigenstates (CLSs), or, modes with nonzero amplitude only at a finite number of lattice sites [5]. This characterization of the flat bands turns out to be closely related to the notion of discrete breathers.

Discrete breathers (DBs) have long been a subject of both theoretical analysis and experimental studies [6, 7, 8]. Recently, it was demonstrated that one can derive exact solutions for DBs in vibro-impact chain models. Such lattices have been investigated analytically both for the Hamiltonian case [9, 10] and for the forced-damped case [11, 12, 10]. Noteworthy in the perspective of the present research are the recent achievements in the analysis of FB structures, especially such as the ones in [3] and [5], and in what concerns the effects produced by the nonlinearity, also the recent success in the analysis of the dynamics and stability of discrete breathers, in specific, the analytic results presented in [11, 12].

In a spatially-extended system comprised of complex elements, the internal symmetry that produces the flat dispersion bands, leads to local dynamics that can exist with no long-range interactions. This means that a single element can be excited locally, which corresponds to the ideal compacton in the discrete lattice.

In continuous systems, compacton solutions with non-analytical profiles were first given in explicit form by Rosenau and Hyman in [13]. They emerge as solutions of a special class of nonlinear partial differential equations, the $K(m,n)$ equations, which possess essential gradient nonlinearity. Although physically motivated, these continuous equations were not derived as a rigorous limit of a specific discrete system [13]. In fact, in an earlier work by Nesterenko, in [14], a partial differential equation with gradient nonlinearity related to nonlinear dispersion, was derived consistently as a continuous limit of a granular mechanical lattice. That work presented a specific solution suggesting the interpretation of being arbitrarily close to having a compact support [15]. An important feature of the solution was that it was characterized by amplitude-independent width. This characteristic, as well as the high level of spatial localization of the obtained soliton solution were validated by experimental observation in a one-dimensional granular medium [16]. Shortly thereafter, in [17], a solitary wave with exactly compact support was shown to exist as a one-parameter class of functions for a specific scenario in the case of nonlinear waves in an elastic rod. In [18], a one-dimensional lattice of atoms with purely anharmonic interactions was considered, a compacton solution was found analytically, and as in [14] was found to be characterized by amplitude-independent width. The compacton solution was also considered in the context of dynamic friction following the Frenkel-Kontorova model [19]. In [20] the mechanical analog of the double-well nonlinear interaction potential that produces compacton-like solutions, was constructed in the form of a lattice of pendula on an elastic foundation with strong coupling. It was found there that static compactons are stable, whereas their mobile counterparts are unstable. In [21], exact results were obtained for linearly stable compactons on discrete 1D lattices with interaction of the Klein-Gordon type. In [22], a discrete chain with linear elastic inter-site interactions with the addition of non-smooth on-site potentials was examined in the context of emergence of compactons.

In the recent work in [23], discrete systems were examined, both in one dimension (a chain) and in two dimensions (a rectangular array), using the Klein-Gordon model with cubic nonlinearities on-site and between the sites. Compact breather solutions were obtained, which are stable and non-traveling, in the continuum limit. In [24] the discrete compact breathers were discussed in the context of Bose-Einstein condensates and nonlinear optics using various interaction models, such as the discrete nonlinear Shrödinger equation, as well as the Klein-Gordon and the Fermi-Pasta-Ulam models. The question of mobility of a discrete soliton in the context of a practical application is addressed in [25] for a chain with a Peierls-Nabarro potential barrier typically applied for a description of the dislocation movement in imperfect crystals. In [26] a discrete mechanical chain with Hertzian interaction is analyzed and traveling compactons with amplitude-independent wave-speed are revealed. In [27], DNA chains are studied using the extended discrete nonlinear Shrödinger equation and solitary waves with compact support are obtained. Some of the obtained compactons are robust and some quickly decompose, depending on the DNA molecule model stacking parameters. A recent work, Ref. [28], deals with the problem of compactification by tuning parameters, such as the measure of interaction nonlinearity (the exponent in a power-law dependence), on-site energies and anisotropic coupling. Compact solitary waves are then obtained for the three lattice-sites support, in the context of the discrete nonlinear Shrödinger equation. Stability analysis of the obtained compact solutions is performed numerically, by approximate application of Floquet theory.

It should be clarified that the compactons in the FBSs differ strongly from the compactons in the systems with the nonlinear gradient coupling, mentioned above. In the latter models, the solutions with the compact support appear only in the continuum limit; the discrete counterparts are characterized by hyper-exponential localization, but are, strictly speaking, not compact. The FBS compactons are strictly compact even in the discrete models, due to the symmetry-related destructive interference.

In [29], the flat bands and the compact states were observed in 2D photonic kagome lattices, with possible application to distortion-free image transmission. In [30], 1D and 2D diamond lattices were studied in the context of kinetic-energy-driven ferromagnetic transition in the Hubbard model, using perturbation theory. The transition was found to be related to perturbed flat bands acquiring small dispersion allowing originally localized electrons to correlate. In [31], a photonic lattice is analyzed, where a flat band emerges due to PT (parity time) symmetry for specific gain and loss values. In [32], a similar system is considered. It is shown there that both randomly positioned states and pinned states at the symmetry plane in the flat band can undergo thresholdless PT symmetry breaking. Furthermore, the flat band is shown to be unaffected by weak disorder, contrary to the modes in the dispersive bands.

In [33], a rhombic 1D optical waveguide lattice was analyzed taking into account cubic nonlinearity, and FB solutions were found to exist yet stay unstable due to mode-coupling until a power-threshold is reached. In [34], a dimerized Lieb lattice with next-to-nearest-neighbor couplings was studied. A flat band was found, corresponding to a point-defect zero-mode, spectrally isolated and spatially localized, producing sub-lattice polarization, which provides a route to mode selection via loss-imbalance. The revealed mechanism may be useful for photonic mode-shaping and guiding, and can, possibly, be extended to laser settings. In artificial 2D materials, the symmetry-breaking-originated flat-band structure provides a mechanism for the design of robust defect-states. In [35], scattering theory was applied to collisions in highly degenerate (flat) bands. It was shown how to construct localized states in flat bands that are the scattering states of impurity-potentials in one dimension, particularized to the sawtooth lattice. In [36] it was found that a planar composition of 2D square photonic lattices with Kerr nonlinearity (nonlinear dielectrics), with rotations corresponding to Pythagorean triples or, alternatively, co-rotated hexagonal lattices, produces localization-delocalization-transition with stable thresholdless solitons.

In [37] a stable flat band in a rhombic quasi-1D lattice optical waveguide is observed in the presence of external driving. Robust localization is found to originate from destructive interference associated with the symmetry in the lattice geometry. Excitation of dispersive bands then leads to Bloch oscillations and coherent destruction of tunneling. In [38], a method is presented for the construction of different single and multiple flat bands using 1D optical lattices with configurational symmetry, such as the sawtooth and rhombic chains. Cubic Kerr-effect-related nonlinearities and the continuum-limit are addressed. Next-to-nearest-neighbor and other extra-interactions are used to create localization for specific power and coupling values. This is an example of existence of compactons in discrete nonlinear quantum lattices, established numerically.

One additional result is relevant for the perspective of the present work. In [21], perfectly compact solutions are obtained analytically in discrete lattices for artificial on-site and inter-site potentials. That is, the potentials are not inspired by physical reasoning but rather are designed to produce compact solutions, the so-called compactlets. This result is noteworthy since it is a unique example, disregarding the present work, where compact solutions are obtained analytically for discrete nonlinear lattices, albeit artificial ones. Moreover, linear stability analysis is performed numerically, in an approximate fashion. A survey of different results related to local-support solutions in the aforementioned physical systems is given in [39].

The objective of the present work, in the discussed context, is to obtain analytic results for the existence and stability of perfectly compact solutions in physically reasonable discrete mechanical lattices, possible asymptotically.

In [40], FB generators were constructed for one-dimensional settings. It is shown there that there exist certain classes of lattices having configurational symmetry that lead to FBs. A typical one is the so-called cross-stitch lattice.

As mentioned, the compacton modes exist in FB systems already in completely linear settings, and persist if the added nonlinearities respect the symmetry. In the nonlinear setting, however, the compactons can become unstable. The stability analysis of such periodic solutions is rather complicated, as it involves direct computation of the monodromy matrix. This paper addresses certain special types of mechanical lattices with FBs. Nonlinear interactions in the considered models lead to possible instability of the compactons, and the exploration of this instability is one of our main research goals. We circumvent, to some extent, the difficulties related to the stability analysis, by taking, in the first instance, the nonlinearity to be introduced into the considered system in the form of impact constraints, which can easily be realized in mechanical lattices. We use an augmented version of the 1D cross-stitch lattice, and analyze the stability of the compacton solutions. The principal idea is to use the methodology of Floquet analysis employing analytically constructed saltation matrices, the way it was done for the case of discrete breathers, as in [11, 12]. Both numerical and asymptotic methods of analysis, available for models with vibro-impact potential, are employed herein for the compactons and the counterpart DB modes, to examine stability. The two basic instability mechanisms revealed in the vibro-impact cross-stitch lattice, namely the local and global instabilities, appear also in a chain of massless boxes with internal nonlinear oscillators, which was examined in the second part of this paper, following the insight gained in [41].

It was recently shown in [41], that in a discrete mechanical system comprised of a massless box in harmonic on-site potential, having two identical oscillators inside, characterized by two-term (linear-cubic) interaction forces with the box, there exists a flat band corresponding to antisymmetric dynamics of the internal masses. The compact dynamics remains stable in a finite parameter range in the nonlinear regime. A (more or less) rigorous stability analysis of the compact solution is enabled, in part, due to the choice of a massless box and the type of coupling. The single-element result can be naturally interpreted to apply to the anti-continuum limit (weak normalized inter-site potential) in a chain of such elements. The present work shows also, by numerical integration, that the compacton solution remains stable also for a chain of such elements, for finite inter-site stiffness values.

The structure of the paper is as follows: Section 2 introduces the augmented version of the cross-stitch model with non-smooth nonlinearities, Secs. 3-4 present the complete stability analysis and numerical verification of the results for the non-smooth model; Sec. 5 addresses the chain of massless boxes with internal nonlinear oscillators, related to the system analyzed in [41]. Section 6 establishes the existence of stable compacton solutions for the chain by direct numerical integration, and is followed by Sec. 7 with concluding remarks.

2 Description of the model– (augmented) cross-stitch lattice with vibro-impact on-site potential

As a model system, we employ the one-dimensional cross-stitch lattice discussed in [38] and [40], with the addition of identical shear springs of generically different stiffness coefficient, connecting the two chains site-to-site (rather than diagonally, as the cross-stitching). Furthermore, the on-site vibro-impact potential with unit restitution coefficient (elastic impacts) is added to each mass in the lattice, similarly to the case in [12]. Geometric nonlinearities are neglected, and only the motion in the axis direction is considered. The system is sketched in Fig. 1 below.

The equations of motion for the cross-stitch lattice given in Fig. 1 are as follows (where $x_{n}(t),y_{n}(t)$ are the horizontal displacements of the lower and upper chain elements, respectively):

[TABLE]

2.1 Linear dynamics

Taking the difference of the first two equations in Eq. (1), one obtains the equation of motion for the detached mode corresponding to the flat band:

[TABLE]

[TABLE]

The flat band frequency is given by ($k_{s}\triangleq k_{\text{shear}}$):

[TABLE]

Taking the sum of the two first equations in Eq. (1), one obtains the known equation of motion for a linear chain (with $2k$ for link stiffness):

[TABLE]

where $u_{n}\triangleq x_{n}+y_{n},\ \forall\ n\in\mathbb{N}\leq N_{s}$. The following dispersion relation stems from Eq. (5):

[TABLE]

Dispersion curves are presented in Fig. 2, where one can observe intersection between the flat and the dispersive bands (no self-avoidance, as there is no linear coupling, unlike in [42]).

To specify linear dynamics beyond dispersion analysis, the system (in the linear regime) can be represented in the following dimensionless dynamical systems’ vector form:

[TABLE]

where $()^{\prime}$ denotes differentiation with respect to time that is scaled by $\sqrt{m/k}$, and:

[TABLE]

in which 0 is a matrix of zeros of size $N\times N$, I is the identity matrix of rank $N$, and:

[TABLE]

As aforementioned, $\gamma\triangleq k_{s}/k$, and:

[TABLE]

Also, $N=2N_{s}$, $N_{s}$ being the number of sites.

2.2 Nonlinear dynamics (regime of impacts)

The nonlinearity enters into the play when some of the lattice elements are engaged in conservative impacts, described by the following relation:

[TABLE]

with $t_{i}$ being the temporal series of impact instances (measured in units of $\sqrt{m/k}$), and where:

[TABLE]

with $j=1,2$ denoting the index of the particle at the localization site.

3 Linear stability analysis

3.1 General framework

To analyze the stability of periodic (modal) solutions, the monodromy matrix can be represented as follows [43] (see A for details):

[TABLE]

where $\hat{\omega}=\omega/\sqrt{k/m}$ is the normalized frequency of the periodic solution ($\omega$ being the frequency) and:

[TABLE]

[TABLE]

where $n_{i}^{(j)}=\delta_{ij}$ for $j=1,2$ and $i=1,2,..,N$.

With no loss of generality, impulsive initial conditions are assumed:

[TABLE]

where Eq. (10) is used and where $x^{0}_{i}$ is a phase-space state vector, in which displacements are normalized by $\Delta$ and velocities – by $\Delta\sqrt{k/m}$.

3.2 Analytical integration of the compacton mode

For the compacton mode, using [12], one gets:

[TABLE]

where $v_{i}(t)$ is as defined in Eq. (3) normalized by $\Delta$ (thus $,y_{i}(t)=-x_{i}(t)=-\frac{1}{2}v_{i}(t)$), and where:

[TABLE]

[TABLE]

with: $\{z\}_{\text{frac}}\triangleq z-\lfloor z\rfloor,\chi\triangleq\omega t-\frac{3\pi}{2}$, $\kappa\triangleq\omega_{\text{FB}}^{2}/\omega^{2}$.

The eigenvalues of the monodromy matrix constructed as described in Sec. 3.1 for the compacton solution, cannot be computed analytically for the complete set of the parameter values (at least, we have not succeeded in doing that). However, it is possible to obtain asymptotic expansions of these eigenvalues for the asymptotic limits considered below.

3.3 Stability boundaries of the compacton mode – asymptotic results

The case of $\gamma\ll 1$. It is possible to obtain analytic approximations for the stability bounds of the compact solution by taking advantage of the compactness. For this sake, a chain of three two-mass elements is considered. The central element corresponds to the excited part of the lattice. The two neighbors are the minimum required to represent periodic boundary conditions symmetrically. The monodromy matrix for this three-element system, containing six masses and having the dimension of $12\times 12$, can be constructed analytically. Furthermore, the matrix can be expanded for small values of $\gamma$ – a third-order Taylor series proves sufficient. Order-balancing then requires that the near-critical dimensionless quarter-period of oscillations, $\hat{q}(\gamma)$, be expanded in powers of the stiffness ratio,$\gamma$, as well (it cannot be large since a finite upper existence bound holds from linear dynamics).

Next, exact analytical calculation of the spectral norm of the asymptotic expansion of the monodromy matrix in radicals is performed (the degeneracy renders the 12th-order polynomial characteristic equation exactly solvable). In the unstable case (identified a posteriori) two eigenvalues are equal to unity, there are four conjugate pairs on the unit circle and a pair of real reciprocals. Expanding the smooth portions of the radicals for the real larger-than-unity eigenvalue as four-term positive-power expansions in $\gamma$ (less proves insufficient) and calculating the limit, yields the following (critically-unstable-case) asymptotic expression:

[TABLE]

where $q_{0}$ emerges from the orders-balancing dimensionless quarter-period expansion, namely:

[TABLE]

The critical value rendering the spectral norm of the monodromy bounded by unity for $\gamma\ll 1$ is thus:

[TABLE]

which corresponds to the following stability-bounding asymptotic relation in the frequency–stiffness-ratio plane:

[TABLE]

Since this critical value corresponds to a transition from a pair of complex conjugate eigenvalues on the unit circle to a pair of real reciprocal eigenvalues, one would expect a pitchfork bifurcation in the underlying chain at this stability bound.

Once the critical value is obtained for the three-site construction, we verify that it remains the critical value for stability also for arbitrary (longer) closed chains. In other words, it holds that the obtained radical remains a root also for the higher-order characteristic equation emerging for the larger chain. Consequently, the obtained critical bound should be regarded not as a feature of the three-site chain, but rather as an asymptotic result true for an arbitrary chain (including an infinitely long one, where closed and open chains are equivalent).

For the case of vanishing shear stiffness, the classical cross-stitch lattice (with the addition of vibro-impact on-site nonlinearity) is reproduced, and for this case, the above result serves as a proof of linear stability of a (one-site) compacton for an arbitrarily long closed chain for all frequencies.

Figure 3 shows the agreement between the numerical calculation of the stability map for the compacton solution, and the asymptotic estimate.

The case of $\gamma\gg 1$. The $\gamma\gg 1$ case is the anti-continuum limit in which the chain decouples and one effectively has a single elementary cell including the two masses in vibro-impact potential, connected by a linear shear spring. The stability of the antisymmetric mode of vibration can be investigated analytically. It turns out that two eigenvalues of the corresponding monodromy matrix are equal to unity and two additional eigenvalues are real positive and have a product of unity. The smallest of these two is treated here for the sake of better presentation.

It appears that this eigenvalue depends only on a certain combination of the two parameters of the problem, namely $\gamma$ and $\hat{\omega}$, and not on both of them independently. A plot of the dependence of the critical eigenvalue on that combination reveals that the system is always unstable in this limit, as can be learned from Fig. 4. Moreover, calculation of the eigenvector of the monodromy matrix of the full chain (in the middle of the feasible range-of-change of the aforementioned combination of parameters) reveals strict localization (as one might expect by naïve understanding of the considered limit). This is shown in Fig. 5 .

3.4 Numeric stability analysis for a long cross-stitch lattice

Numerical spectral analysis of analytically constructed monodromy matrices for different stiffness-ratio/frequency pairs, reveals the existence of a growing stable stripe, with upper and lower bounds, for decreasing shear-spring stiffness. The upper (frequency) stability limit corresponds to the pitchfork bifurcation, and physically is related to the internal instability of the antisymmetric mode in the unit cell.

This conclusion is confirmed by almost perfect coincidence between the numeric simulation of the complete system and the analytic results for the single unit cell. The lower stability limit corresponds to the Neimark-Sacker bifurcation; this instability scenario is absent in the unit cell and should be attributed to interaction with the propagation zone of the lattice. The bounds, along with the corresponding asymptotic estimates and the flat band as a reference, are shown in Fig. 3.

3.5 The alternative localization regime – discrete breather (analytical DB-mode integration)

As was established above, each unit-cell in the lattice can admit the antisymmetric mode, corresponding to the compacton solution. However, it can also oscillate in the symmetric mode. In that case, no destructive interference occurs, and thus one expects to observe the conventional discrete breather. Similar breathers were treated, for instance, in [10]. Here, a more complete form is derived. The effective breather displacement, namely, $u_{n}(t)\triangleq x_{n}(t)+y_{n}(t)$ is given by a symmetric Fourier series (with the associated auxiliary quantities):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Delta$ is the amplitude of the breather, or half the distance between the limiters. In addition,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

(recall that $\hat{\omega}\triangleq\sqrt{\frac{m}{k}}\omega$). The initial conditions for the stability analysis of the breather mode are zero displacements with consistent initial velocities:

[TABLE]

and the following consistent (inhomogeneous) initial condition, with $t_{0}=\pi/(2\omega)$:

[TABLE]

3.6 Stability boundaries of the DB mode – asymptotic results

The case of $\gamma\ll 1$. Analytic estimation of the stability bounds of the breather mode in the $\omega$–$\gamma$ plane can be obtained asymptotically for $\gamma\ll 1$. To this end, once again, the previously discussed triple of elementary cells is considered. This is the smallest chain needed to represent a symmetric system with periodic boundary conditions. As in the case of the compacton mode, the stability criterion here prescribes an order of unity critical value on a combination of the parameters. One of these parameters is assumed asymptotically small, hence it is expected that the other parameter (here the frequency) would be of order different than unity in the critical case. Clearly, the frequency cannot be much smaller than the linear-regime limit, and thus a large critical value is expected. Unlike in the case of the compacton mode, where there is only one velocity value in the analysis, which (thus) cancels out, in the case of the breather, the exact initial velocity profile is required. For the asymptotic analysis considered here the high-frequency limit of Eq. (34) for $N_{s}=3$, as obtained by rigorous asymptotic treatment, is sufficient and reads, ($\hat{q}$ from Eq. (21)):

[TABLE]

Having now all the quantities required for the construction of the monodromy matrix (as functions of $\gamma$ and $\omega$), one can obtain it analytically. Taking $N_{s}=3$, just as for the compacton case, constructing the monodromy matrix and taking its two-variable second-order Taylor-series approximation with respect to the small parameters $\gamma$ and $\hat{q}$ (as an equivalent but alternative route to the one taken in the compacton case), and then computing the spectral norm exactly and further expanding the smooth parts of the radicals in three leading terms, and order-balancing, one obtains the real maximal eigenvalue of M (for the critically unstable case) as:

[TABLE]

where $q_{1}$ is given by:

[TABLE]

and

[TABLE]

The critical value rendering the spectral norm of the monodromy matrix bounded by unity is the same as for the compacton:

[TABLE]

This corresponds to the following stability-bounding asymptotic relation in the frequency–stiffness-ratio plane:

[TABLE]

One notes that in the $\gamma\to 0$ limit, the same curve in the $\omega$–$\gamma$ plane is obtained as the stability limit for either the compacton or the breather, with the compacton mode being stable below that line and the breather mode being stable above it, as arises from Eqs. (4,6,20,23,36,40).

The case of $\gamma\gg 1$. In the limit of large $\gamma$, the cross-stitch chain degenerates to the single two-mass ”molecule”. The stability of the symmetric mode of vibration of such a molecule can be examined analytically. It appears that once the assumption of a large $\gamma$ is introduced, allowing one to consider a single two-mass element, no further asymptotics is required, and the eigenvalues of the monodromy matrix can be obtained in closed-form. Two eigenvalues are identically unity, which represents the degeneracy in position and velocity arising from the symmetry. Two additional eigenvalues are generally complex conjugates with absolute values of unity. Those eigenvalues are the critical ones, and they are given by:

[TABLE]

This result is a first-order approximation in $\gamma^{-1}$. Clearly, higher terms would perturb those eigenvalues. For eigenvalues with nonzero imaginary parts, the perturbation could only change the arguments of the eigenvalues in the complex plane, since the symmetry would keep them complex conjugates and energy conservation would keep the product of their radii at the value of unity. Hence, they would both remain on the unit circle, and no instability would arise. However, for the special case where all four eigenvalues are equal to unity, a small perturbation can drive two of the eigenvalues apart on the real axis, keeping their product unitary. In this case, one eigenvalue would be real and greater than unity, corresponding to instability. Therefore, an infinite (countable) set of narrow instability tongues (turning to lines for $\gamma\to\infty$) would emerge on the $\omega$–$\gamma$ plane in the perturbed (chain) case, obeying the following asymptotic relationship:

[TABLE]

The quality of the asymptotic estimates for the stability bounds is shown in Figs. 6 and 7. Figure 7 also features a lower stability bound, corresponding to something that can be called double pitchfork bifurcation. In this case, both chains experience the pitchfork-type loss of stability simultaneously, which implies the destruction of periodicity for the breather. This differs from the regular pitchfork bifurcation emerging at the upper bound, where stability is lost due to internal resonance, where part of the energy is transmitted from the symmetric to the antisymmetric mode.

3.7 Stability bounds for the DB mode – numeric results

4 Numerical integration

To validate the results of the stability analysis numerically, we employ an algorithmically stable integration scheme, in conjunction with a smooth model representing the vibro-impact potential. The widely-known high-odd-power force term is used for modeling impact. The solver employed here is a one-step method based on the trapezoidal rule and a backward differentiation formula of order 2 with a free interpolant [44].

4.1 Simulation of the compacton in and below the stability region

For numerical integration, the system is presented in the following form:

[TABLE]

where $\hat{\textbf{A}}$ is defined in Eq. (9). The displacements vector is: $\textbf{w}=(x_{1},y_{1},x_{2},y_{2},...,x_{N_{s}},y_{N_{s}})^{\top}$, and the nonlinear impact-representing forcing vector is given by:

[TABLE]

where $\xi$ is the force exponent, reproducing impact when $\xi\to\infty$. As in [12], $\xi=300$ was chosen here. Figure 8 comprises, along with Fig. 3, perhaps the most important result in the present paper. It confirms the existence of a stable compacton in a nonlinear 1D lattice, by showing the dynamic boundedness of a perturbation to a perfectly compact solution. This is illustrated both by the time histories of the principal-site displacements and by the chain mode-amplitude profiles.

One can argue that the dashed line in Fig. 3 gives a complete stability picture for the asymptotic limit, the solid (red online) curve therein confirms that one can rely on this limit, when the parameter values are finite, and Fig. 8 confirms that linear stability analysis provides reliable results for actual small yet finite perturbations.

Figure 9 shows that when the frequency is smaller than the critical (lower) bound, delocalization occurs, through the Neimark-Sacker (NS) bifurcation, leading to, arguably, chaotic-like trajectories. The initial state there is the same slightly perturbed compacton as in the case exhibited in Fig. 8.

In both cases, the shear spring is taken to be an order of magnitude less rigid than the other springs. This is required for a stable solution to be feasible according to linear stability analysis.

4.2 Simulation of the discrete breather in and below the stability region

The results for the numerical integration of the discrete breather are given in Figs. 10 and 11.

Figure 10 confirms the stability of the breather for frequency above the lower threshold (noteworthy is the absence of an upper stability bound). Figure 11 illustrates the loss of stability of the breather mode for a frequency lower than the critical value. The long-time dynamics of the perturbed symmetric mode is not “chaotic” as the one emerging from the ’NS’ bifurcation shown in Fig. 9, but rather reminds (quasi) beating typical for local energy exchange between two nonlinear normal modes (NNMs). A similar picture is observed below, when crossing the ’PF’ instability boundary to the unstable side for the compacton mode.

4.3 Integration with compact initial conditions above the stable region – global dynamics

Asymptotic analysis reveals that for the breather:

[TABLE]

(where $\hat{q}$ is as defined in Eq. (21)), which implies the asymptotic compactness of the breather mode. Moreover, the compacton mode is strictly compact anyway for all frequencies. Consequently, for high frequencies both modes are compact. Therefore, the entire dynamics is compact and is thus expected to be equivalent to the dynamics of a single element for high kinetic-to-potential energy ratios.

In [45], a similar single element is analyzed and global dynamics is revealed through transformation to canonic variables (action–angle formalism). One result arising in [45] is that the one-element cross-stitch system has a saddle point corresponding to the antisymmetric mode. This means that for high energies and finite stiffness-ratio values, the compacton solution should be unstable. Furthermore, the aforementioned saddle point is a global feature, and the nearest minimum corresponds to the symmetric mode. However, this does not yet imply that starting at the compacton initial conditions for high enough energies, one would stabilize the system in the symmetric mode. The reason is that, as shown in [45], there is no equal-energy contour line passing through both the compacton-related saddle point and the breather-related minimum. Such stable global dynamic a route is hence inexistent. What does exist is a contour line passing through the saddle point, reaching an intermediate location and going back to the saddle point. The intermediate location in the canonic phase-space is close enough to the symmetric mode. This dynamics correlates well-enough with the results of the numerical integration of the chain, given in Figs. 12 and 13.

5 The second model system (1D chain) – two-masses in a massless box with smooth internal potential and linear links

Figure 14 presents the model setting for a chain of massless boxes with internal oscillators. For the required purposes, the interesting case would be that in which in every box the masses and spring characteristics are identical for the two oscillators (the reasoning is explained below). For such a case, the equations of motion for the system are as follows:

[TABLE]

where $\bar{u}_{n}$ and $\tilde{u}_{n}$ denote the displacements of the first and the second oscillator, respectively, for each box, and $u_{n}$ stands for the displacement of the $n^{\text{th}}$ box itself. The conservative case (no external forcing and no damping) is assumed. The parameters $k$ and $k_{1}$ stand for spring rigidities inside and outside each box, respectively, and $f^{NL}_{(x)}$ is an interaction force derived from a smooth anharmonic correction to on-site potential, representing local nonlinearity.

5.1 Linear dynamics

The set of equations given above can be written in normalized form, more convenient for mathematical treatment, which in the linear case becomes:

[TABLE]

(here $\kappa$ and $\gamma$ stand for spring rigidities normalized by mass, inside and outside each box, respectively).

As in the previous model system, due to the symmetry between the first and second equations in Eqs. (47), linear transformation of variables is employed, as follows:

[TABLE]

Combining Eqs. (47) and (48), the following equivalent system of equations is obtained:

[TABLE]

The first equation in (49) decouples from the rest of the equations. Moreover, one observes that the first equation in (49) is a local relation, and indeed oscillatory spectral analysis results in a dispersion relation for $v_{n}(t)$ that has the following form:

[TABLE]

in which $\omega_{0}$ stands for the frequency and $q$ denotes the wave number. Once again, due to the local configurational symmetry, a flat band emerges.

Next, proceeding to the two remaining equations in system (49), one finds that a plane wave oscillatory solution of this system is only possible in the case of synchronous motion of $u_{n}$ and $x_{n}$ (which is reasonable, since oscillatory motion requires mass inertia which is supplied to $u_{n}$ by $x_{n}$). For such motion, the corresponding dispersion relation takes the form:

[TABLE]

where $\delta\triangleq\kappa/\gamma$ and $\omega_{1}=\sqrt{\kappa}\hat{\omega}_{1}<\omega_{0}$ (implying no-crossing, unlike in Fig. 2), $\omega_{1}$ being the frequency. One can observe that in the limit $\delta\gg 1$, the dispersion relation of a simple linear chain is recovered, whereas in the limit $\delta\ll 1$, this dispersion relation becomes (asymptotically) flat too.

Another simple linear mechanical system that has a flat band structure has therefore been demonstrated. To conclude this part of the description, a plot of the oscillatory spectra of the system is presented in Fig. 15, illustrating the flat band in more detail.

5.2 Nonlinear dynamics: analysis of a single (elementary) cell

The basic system subjected to nonlinear analysis is presented in Fig. 16.

5.2.1 Analytical integration of the compact mode

In line with the well-known Fermi-Pasta-Ulam model, cubic nonlinearity is added to the internal interaction forces, producing the following equations of motion (where $p$ measures anharmonicity):

[TABLE]

Since the conservative case is considered, with no forcing or damping, the general expression for the energy of the model system is given by:

[TABLE]

Next, appropriate notation is introduced:

[TABLE]

[TABLE]

where the auxiliary quantities are given by:

[TABLE]

This can be used to represent the energy in the antisymmetric mode in terms of some amplitude $Y_{0}$ (later used, for convenience, in parametric analysis), which can be interpreted as the amplitude of oscillation in the symmetric mode with the same energy:

[TABLE]

The antisymmetric mode corresponds to zero motion of the center-of-mass of the system, which leads, following the definitions in Eq. (54), to the next relations:

[TABLE]

Taking the difference of the first two equations in Eqs. (52), expressed in terms of $y,v$ and $x_{0}$, after eliminating $x_{0}$ using the third, algebraic equation in Eqs. (52), and its first and second time derivatives, expressed in terms of the three aforementioned variables, leads to the following equation (required for stability analysis at a later instance):

[TABLE]

Substituting Eq. (58) into the expression for the total energy of the system given in Eq. (53), and using Eqs. (54), (56) and (57), one obtains the first integral of the antisymmetric mode of the system, in the following form:

[TABLE]

(here $Y_{0}(E_{0},k_{0},k,p)$ is perceived merely as convenient parametrization of the total energy in the system).

Further integration of Eq. (60) with a centered impulsive boundary condition ($\hat{v}(0)=0$), results in a closed form expression for the normalized antisymmetric mode displacement as a function of the dimensionless time, using the fifth Jacobi elliptic function ’sd’, as follows:

[TABLE]

One observes that for small amplitudes, the simple form $\hat{v}(t)=\sin(\omega t)$ is reproduced.

Equation (61) provides the explicit solution for the antisymmetic mode. In order to study its stability, Floquet-Hill method is applied, and Fourier-series representation of the solution is called for. The sine series representation of sd$(z|m)$ is given in [46]. In the employed parametrization it gives rise to the following series for the antisymmetric mode:

[TABLE]

where

[TABLE]

and $K(m)$ is the complete elliptic integral of the first kind (and $\delta_{a,b}$ is Kronecker’s delta).

5.2.2 Linear stability analysis of the compact mode

In order to examine the linear stability of the obtained solution, the equation of motion for the symmetric mode is required.

Finding from definition (54) that:

[TABLE]

and substituting this into the third equation in Eqs. (52), one obtains the (algebraic) equation of motion for the displacement of the box, as follows:

[TABLE]

Next, taking half of the sum of the first two equations in Eqs. (52) and substituting Eq. (64), the third of definitions (54) and the second time derivative of Eq. (65) into the result, eliminating the second time derivative of $v$ by isolating it from Eq. (59), one obtains the equation of motion for the symmetric mode (for mixed-mode dynamics), containing only $v$, $y$, their first time derivatives and the second time derivative of $y$, in terms of the dimensionless quantities defined in Eq. (62) and in the first row in Eq. (63). Consequently linearizing with respect to $y(t)$, one obtains the evolution equation for small perturbations to the antisymmetric mode.

Additional mathematical manipulations, described in [41], produce the explicit Hill form of the linear stability equation for the purely antisymmetric mode:

[TABLE]

with $\Omega_{v}$ and $\hat{v}(\hat{\tau})$ given by Eqs. (62) and (63), and with the Hill function, $h(\hat{\tau})$, being the expression multiplying $\tilde{y}(\hat{\tau})$ in Eq. (66).

Next, stability analysis of Eq. (66) is performed using Hill’s method. The details of this procedure, namely, rigorous deconvolution and the actual implementation of Hill’s method, for identification of instability tongues characterizing the antisymmetric mode, as well as further study of the nature of the occurring instabilities by use of Poincaré sections, are given in B.

Of the obtained results, the one mostly relevant for the consequent analysis of a nonlinear chain, namely a stability map in the parameter space (plane) for the representative element, is shown in Fig. 17.

The single-element analysis presented here is reproduced in part from [41], for better understanding of the dynamics in the context of a chain. This owes to the fact that for a chain, rigorous stability analysis is not within reach, other than when understanding the single element case as representative of the anti-continuum limit for a chain, taking $k/k_{0}\gg 1$. In that limit, as Fig. 17 shows, the compacton is linearly stable (see [41] for details).

6 Numerical integration for the chain

The existence of the compacton in a chain is evident from the fact of compactness itself. Since the compacton corresponds to the antisymmetric mode in the representative element, for which the box displacement is always zero (see Eq. (58)), implying no interaction between boxes, clearly the solution given in Eq. (61) should also hold for a chain – it is then simply the difference of the displacements of the two masses in a given box. In fact, several-site compactons can exist just as well. The only arising question would be that of stability. The case addressed in the present work would be, however, that of the stability of a single-site compacton.

The stability of a compacton in the chain shown in Fig. 14 can be examined by taking parameters inside and outside of the finite-width instability tongue, as calculated for a single element. This embodies the assumption that the stability limits for a chain are somehow correlated with those of a single element. The two points in the stiffness-ratio – energy plane chosen for stability analysis by direct integration of a chain are shown in Fig. 17. In order to perform the numerical integration from the perspective of stability of (nonlinear) normal modes (NNMs), first, the equations of motion given in Eq. (46) for the case of cubic nonlinearity in the internal interaction forces, are rewritten as a dynamic system, in which modal displacements can be integrated directly.

6.1 Equations of motion

The equations of motion for a chain of $N$ elements comprised of massless boxes occupied by two identical internal oscillators each, with linear-cubic interaction of each internal particle with the box, are given below. The equations are expressed in terms of the dimensionless quantities $\hat{v}_{i}$ and $\hat{y}_{i}$, representing the displacement difference of the internal oscillators and the center-of-mass displacement with respect to the box, for each element. In order to interpret the one-element stability results for the chain, the following relation is set: $k_{1}\triangleq k_{0}/2$.

[TABLE]

[TABLE]

[TABLE]

In order to avoid boundary effects, periodic boundary conditions are applied. In addition, in order to avoid singularity, the system (described by the matrix given in Eq. 69) is “statically condensed” by one degree of freedom (related to rigid-body motion). Compacton initial conditions are employed.

6.2 Integration results

Unstable case (moderate amplitude). The first result is obtained by taking the left point in Fig. 17, the one lying well within the instability tongue for a single element. The initial condition is of a slightly perturbed compacton. A stable integration algorithm with a bounded energy-error [44] is employed.

It appears, as shown in Fig. 18, that for the chosen parameters the stability of the compacton solution is lost, and small, yet finite-amplitude center-of-mass relative (to the boxes) displacements emerge for the ringed chain. The packets-pattern in Fig. 18 (top) is typical for energy exchange between different NNMs by the mechanism of parametric resonance (at least initially in the process of loss of stability). Fig. 18 (bottom) shows the long-time history confirming the persistence of the energy-exchange phenomenon after all transients are gone. This energy exchange between modes is similar to what is observed (for other parameter values) on a $y$–$y^{\prime}$ Poincaré section for a single element [41].

Figure 19 shows the flow close to a global bifurcation where KAM (Kolmogorov-Arnold-Moser) islands of the symmetric and antisymmetric modes “collide”, corresponding to the aforementioned significant (1:2 resonance for both modes) energy exchange (see [41]). Figure 20 shows the envelope (amplitude) of the symmetric mode displacement for the chain, which has the typical profile of a breather, with partial localization arising from the (cubic) on-site nonlinearity.

Stable case (moderate amplitude). The second result is obtained when taking the rightmost point in Fig. 17, the one lying well inside the stability region for a compacton in a representative element. Numerical integration is performed starting from a slightly perturbed compacton initial-condition.

Figure 21 shows that when starting from the same perturbation to the compacton as in the unstable-parameters’ case, the perturbation stays bounded, i.e. there is no delocalization. This is illustrated by the boundedness of the central-site history and maximal-displacement-time profile of the symmetric mode, and the antisymmetric-mode central-site displacement history, which exhibits stable periodic dynamics. Figures 18-21 show that, as predicted by linear stability analysis for a representative element, the nonlinearity may either retain or destroy the stability of the compacton solution (the sheer existence of the compacton obviously emerges from the configurational symmetry corresponding to the existence of a flat band in the dispersion spectrum of the system).

It also appears that periodic dynamics in a (closed) chain can correspond to a combination of a perfectly localized mode, namely the compacton, emerging due to configurational symmetry, and a partially localized mode, namely a breather mode, corresponding to localization-inducing on-site nonlinearity.

The following subsection extends the analysis to the case of a high-amplitude (high energy) regime, where an additional noteworthy phenomenon is observed.

6.2.1 High-amplitude regime – unstable case

The next result is obtained by increasing the (dimensionless) amplitude one order of magnitude and finding a corresponding stiffness ratio for which the representative element analysis shows the compacton solution to be unstable. This appears to correspond to highly rigid inter-site link potential, and repeating the integration, the obtained results remind those presented in Figs. 18-20. There appears to be the additional property that the breather solution possesses a higher level of spatial localization, which can be attributed to the higher energy contained in the system. This matter is elaborated on for the case of the non-smooth system in the analysis in Sec. 4.3.

Figure 22 shows the principal-site antisymmetric-mode displacement history for the perturbed case, on top of the result for the unperturbed case. The result exhibits, once again, the energy exchange feature indicative of resonance between the compacton and the symmetric mode. Figure 23 exhibits the symmetric mode displacement’s first energy-packet history, and the symmetric-mode maximal-displacement-time profile (the antisymmetric mode remains perfectly compact). One notes that the extent of energy exchange for the high-energy regime is almost full – considerably higher than the partial energy exchange observed in the moderate-amplitude regime.

The following subsection exhibits the integration results in the high-amplitude regime corresponding to a stiffness ratio well inside the stability region of the representative-element.

6.2.2 High-amplitude regime – stable case

Figures 24 and 25 show the integration results for the slightly perturbed compacton conditions (in solid blue online) on top of the unperturbed ones (in dashed red online).

The stability of the compacton is evident.

7 Concluding remarks

This paper demonstrates that, as one would expect, flat bands and compactons can be realized in purely mechanical lattices. However, mechanical lattices can possess elements that are hardly viable in optical lattices, namely, the impact constraints. This type of nonlinearity allows substantial simplification of the stability analysis; the latter can be accomplished within the framework of linear algebra. Consequently, it becomes possible to distinguish between the different physical reasons for compacton instability. The first one is related to the internal instability of the antisymmetric mode in the unit cell of the lattice. The other is related to resonance with the propagation spectrum. For the local mechanism, it was possible to obtain an analytic evaluation of the local instability threshold as a function of the system parameters. Moreover, it was observed that an unstable compacton can exchange energy with the exponentially localized discrete breather.

Similar considerations help to understand the stability characteristics of the lattice comprised of the massless boxes with the internal smooth nonlinear oscillators. There, unlike for the non-smooth system, the strict results entailing a complete stability picture are limited to the anti-continuum limit. One can argue that the advantage of the proposed framework beyond the aforementioned asymptotic limits is that once the monodromy matrix is constructed analytically, as it is done for the non-smooth system, the dependence of its eigenvalues on the parameters is proven to be regular, since only linear algebra is involved. This is in contrast with fully numerical Floquet analysis required for a general smoothly nonlinear system, which essentially involves numerical integration of a system formally not known to be integrable.

All lattices having flat bands inevitably possess complex unit-cells. Thus, one can expect that, generically, the stability of compactons will be affected by these two physical mechanisms – local instability in the cell, and a global one, due to interaction with propagating linear or weakly nonlinear waves. In the former case, it seems that one can, at least asymptotically, rely on the (much more simple) local stability analysis presented herein.

In addition to analytic stability threshold estimates, closed-form expressions were given here for the compacton solutions, for both systems.

Acknowledgments

The authors are grateful to the Israel Science Foundation (Grant No. 838/13) for financial support.

Appendix A Derivation of the saltation matrix for asynchronous two-impact jumps

First impact. For the first impact, the following relations hold, assuming small impact time (which owes to impact idealization) and small solution-perturbation (pertinent to linear stability analysis):

[TABLE]

where the smooth dynamics is linear and is represented by: $\dot{\textbf{x}}=\textbf{A}\textbf{x}$, the jump corresponding to a single mass impact is represented by the mapping $\textbf{g}_{1}(\textbf{x},t)$, the jump corresponding to synchronous impact of two masses at a two-mass site is represented by the mapping $\textbf{g}(\textbf{x},t)$, the periodic solution is x, the perturbed solution is $\bar{\textbf{x}}$, the perturbation in the solution is $\delta\textbf{x}$, the impact instance of the unperturbed solution is $T$, the first-impact instance of the perturbed solution is $\bar{T}$, and it is assumed (with no loss of generality) that one has $T>\bar{T}$, and the subscripts ‘$+$’ and ‘$-$’ denote limit-values before and after impacts.

Combining the eight equations in Eq. (A.1), and expanding the mappings $\textbf{g}_{1}$ and g as first-order multiple-variable Taylor series (which is justified for the rather general case of a smooth mapping representing impact), and retaining only first-order terms, the following relation is obtained:

[TABLE]

where $\delta{t}=T-\bar{T}$ is obtained as a function of $\delta{\textbf{x}}_{-}$ by use of the jump conditions, as follows. The jump conditions for the periodic and perturbed solutions as related to the first impact are given by:

[TABLE]

Substituting lines 5, 8 and 1 of Eqs. (A.1) into line 2 of Eq. (A.3), performing a first-order Taylor series expansion, and employing line 1 of Eq. (A.3), one obtains an equation relating $\delta{t}$ and $\delta{\textbf{x}}_{-}$, which reads:

[TABLE]

Substituting Eq. (A.4) into Eq. (A.2), one obtains the following relations:

[TABLE]

[TABLE]

Second impact. For the second impact, the following relations hold:

[TABLE]

where the jump corresponding to the second mass impact (for the perturbed periodic solution) is represented by the mapping $\textbf{g}_{2}(\textbf{x},t)$, the perturbation in the solution after both impacts occurred in both the perturbed and unperturbed solutions is $\delta\textbf{x}_{++}$, the second-impact instance of the perturbed solution is $\hat{T}$, and one has $\delta\bar{t}=\hat{T}-T>0$.

Combining the six equations in Eqs. (A.7), expanding the mapping $\textbf{g}_{2}$ as a first order multiple-variable Taylor series, retaining only first-order terms, and recalling line 6 of Eqs. (A.1), and Eq. (A.5), one obtains the following relation:

[TABLE]

where $\delta{\bar{t}}$ is obtained as a function of $\delta{\textbf{x}}_{-}$ by use of the jump conditions, as follows. The jump conditions for the periodic and perturbed solutions as related to the second impact are given by:

[TABLE]

Substituting lines 1, 2 and 6 of Eqs. (A.7), and Eq. (A.5) into line 2 of Eq. (A.9), performing a first-order Taylor series expansion, and employing line 1 of Eq. (A.9), one obtains the equation relating $\delta{\bar{t}}$ and $\delta{\textbf{x}}_{-}$, which reads:

[TABLE]

Substituting Eq. (A.10) into Eq. (A.8), one obtains the following relations:

[TABLE]

[TABLE]

The crux of the present derivation is embodied in Eq. (A.11). The linear stability analysis performed here assumes small perturbation to the solution prior to impact, and linear stability can only correspond to the perturbation staying small after the sequence of impacts has occurred. Consequently, a necessary condition for both the possibility of linear stability and the self-consistency of the linear stability analysis as a procedure is that the first two terms in the right-hand side of Eq. (A.11) cancel out (formally they can result in an infinitesimal quantity, yet they contain no parameter taken to the limit zero in the context of the linear stability analysis and thus can only result in a finite value that thus has to be an exact zero), namely:

[TABLE]

This consistency requirement can be termed “impact-integrability”, since it implies that a series of consecutive impacts produces the same result as a single synchronous impact. The collision between two particles and a fixed wall is always impact-integrable. However, in the case of a simultaneous many-body impact, integrability may hold only given some restriction on the problem parameters, solely under which slight asynchronicity will not produce finite change in the outcome. This intuitive observation is represented by a formal linear stability analysis as shown above.

After having been revealed in essence for the two-impacts case, the requirement can be straight-forwardly generalized to $n$ consecutive impacts, the induction procedure being trivial here. For a series of $n$ impacts, the “impact-integrability” requirement, or the necessary condition for stability, is:

[TABLE]

If indeed the system in question is impact-integrable, then it holds that by substituting Eq. (A.13) into the first term in the parentheses in the expression for $\textbf{S}_{2}$ in Eq. (A.12), it becomes apparent that the matrix $\textbf{S}_{2}$ is instantaneous and only depends on the value of x after the first impact. Then it becomes clear that $\textbf{S}_{2}$ is a legitimate saltation matrix and one has:

[TABLE]

for the two-impacts case, and:

[TABLE]

for the $n$-impacts case.

After establishing the impact-integrability requirement as the necessary condition for stability, it is left to note that in the setting assumed in Secs. 2-4, there are only impacts between a mass and a fixed wall, albeit multiple ones. Hence the aforementioned condition is satisfied automatically, and Eqs. (A.12), (A.13) and (A.15) apply.

Also noteworthy is the following. The partial temporal derivatives of the jump conditions and the jump mappings appearing in Eqs. (A.6) and (A.12) can generally be non-vanishing, as, say, for the case of a periodically externally moved wall and a periodically externally controlled restitution coefficient (as in magnetorheological composites). Hence the corresponding terms (with the associated signs) may be of relevance. However, in the case considered in Secs. 2-4, the walls are fixed and perfect restitution is assumed, for which the jump mapping is realized as a multiplication by a matrix. Consequently, Eqs. (13)-(15) hold.

Appendix B Details of the stability analysis of Eq. (66)

Here the procedure of the stability analysis of Eq. (66) is described. Asymptotic analysis of Eq. (66), resulting in the identification of the number of emerging instability tongues, derivation of starting points at the zero-amplitude limit for numerical instability-tongue boundaries calculation, as well as the derivation of leading-order asymptotic expansions for the boundaries of the first two instability tongues, using both Mathieu and higher-order approximations, can be found in [41].

B.1. Fourier decomposition by numerical deconvolution

In order to perform stability analysis by Hill’s determinants method, the Hill function in Eq. (66) has to be expressed as a Fourier series. To this end, first, the numerator and the denominator are decomposed, separately. Both the numerator and the denominator depend linearly on the square of the normalized displacement in the antisymmetric mode, $\hat{v}(\hat{\tau})$, which was already expressed in a sine series form in Eq. (62). Thus, the first thing one has to do is to perform convolution, since the Fourier series decomposition of the product of two functions (a function and itself, in the present case) equals the (discrete) convolution of the Fourier series decomposition of each of the functions. Multiplying, rearranging and adjusting the indices, one obtains the following expression:

[TABLE]

Next the definition of the Hill function, as arising from Eq. (66), is rearranged to a product form, $f_{1}(\hat{\tau})h(\hat{\tau})=f_{2}(\hat{\tau})$, and then $h(\hat{\tau})$ is formally expanded into a Fourier cosine series (as the ratio of two even functions, $h(\hat{\tau})$ is always even):

[TABLE]

Substituting Eqs. (B.1) and (B.2) into the aforementioned product-form equation for $h(\hat{\tau})$ that emerges from Eq. (66), expanding the resulting products of sums, using trigonometric identities to turn cosine products to single index-shifted cosines and performing the convolution, one finally obtains the product-form equation for $h(\hat{\tau})$ as a single Fourier cosine series that should be equal to zero at all times, namely,

[TABLE]

Fulfillment of this condition corresponds to the vanishing of the coefficients of this Fourier series. The equations representing the vanishing of these $total$ Fourier cosine series coefficients, formulated as equations for the Fourier cosine coefficients of $h(\hat{\tau})$, take the form of linear algebraic equations.

First, there is the condition $G_{0}=0$, producing the following explicit equation for $H_{0}$ in terms of $H_{n\geq 1}$:

[TABLE]

Second, taking the set of conditions $G_{n\geq 1}=0$, which contain both $H_{0}$ and $H_{n\geq 1}$, and eliminating $H_{0}$ from these equations by isolating it from Eq. (B.4), one obtains an infinite set of linear inhomogeneous equations for $H_{n\geq 1}$, which takes the following form in indicial notation:

[TABLE]

where

[TABLE]

Finally, Eqs. (B.5) are truncated, which is achieved by exchanging $\mathbb{N}$ with $\mathbb{N}\leq\hat{N}$ everywhere in these equations (this is the reason for the presence of Kronecker’s delta in the first term in the third line in Eq. (B.6)), and deconvolution is performed by solving the truncated system numerically, by a standard algorithm for the solution of linear systems, thus obtaining $H_{n\geq 1}$. The coefficient $H_{0}$ is then derived using Eq. (B.4).

The issue of convergence for this procedure is discussed in [41].

Having obtained the Fourier cosine coefficients of the Hill function, one performs the stability analysis by Hill’s (truncated) infinite determinants method. The result is shown in Fig. 17.

B.2. Asymptotic stability-limit estimates

The small-amplitude analytic estimates for the stability boundaries (given for energy in amplitude units, following Eq. (57)), are (see [41] for derivation):

[TABLE]

Plots of the first (finite-width) instability tongue boundaries for the antisymmetric mode, and the second tongue, which is (at least to second order) a collapsed, zero-width one, are presented in Fig. 17. Asymptotic approximations for the boundaries, as given in Eq. (B.7), are shown therein in dashes.

B.3. Poincaré sections

Figure 26 presents Poincaré sections of the flow in the vicinity of the second (degenerate) instability tongue. The emergence of KAM islands represents the effect of the loss of stability of the compact mode through parametric resonance.

The coordinates in the plots are the symmetric mode displacement and its (dimensionless) time derivative. Hence, energy localized entirely in the compact, antisymmetric mode is represented as a point in the center (or an oval around the center, in the perturbed case). As stability is lost, the oval around the center turns into four KAM islands, corresponding to the 1:2 resonance of the second instability tongue.

The first, finite-width instability tongue is standard manifestation of instability through parametric resonance – hence the derivation of the stability range boundaries seems sufficient. The second tongue appears to be a collapsed one, and therefore further investigation is called for, in order to understand the finite-amount topological implications (for parameter values at finite distance from the zero-width tongue) of the appearance of such a feature.

To this end Poincaré section analysis was employed, as a means to reveal the existence of instability through higher-order resonance manifested in the appearance of KAM tori (of course, the first instability tongue has its own second validation – by direct numerical integration, albeit for a chain, as performed in Sec. 6).

Figure 27, much like Fig. 19, shows resonance-originating KAM islands (increasing with energy), confirming the existence of periodic dynamics of energy exchange between the different modes, for higher energies, for the representative element. The same interpretation can be attributed to the non-divergent dynamics of energy exchange in the case of a chain, as observed in the histories of the principal-box (initial localization site) masses, as shown in Fig. 18.

The choice of a massless box in the conservative case allows, at least for a representative element, to describe regular (non-chaotic) dynamics completely, using Poincaré sections.

The emergence of KAM islands spanning from the center to the envelope implies the instability of the compact solution and corresponding significant energy exchange with the symmetric mode (this scenario corresponds to the high-energy points inside the first instability tongue). Chain dynamics can then be inferred – approximately for arbitrary link-stiffness, with the help of direct numerical integration, and exactly in the anti-continuum limit.