Stability of compact breathers in translationally-invariant nonlinear chains with flat dispersion bands
aa r X i v : . [ n li n . PS ] A ug Stability of compact breathers in translationally-invariant nonlinear chains with flatdispersion bands
Nathan Perchikov, O.V. Gendelman
Faculty of Mechanical Engineering, Technion, Haifa 32000, Israel
Abstract
The paper addresses compact oscillatory states (compact breathers) in translationally-invariant lattices with flat disper-sion bands. The compact breathers appear in such systems even in the linear approximation. If the interactions arenonlinear, but comply with the flat-band symmetry, the compact breather solutions exist, but can lose their stability forcertain parameter values. As benchmark nonlinear potentials, we use the β -FPU (Fermi-Pasta-Ulam) and vibro-impactmodels. Loss of stability is numerically observed to occur through either pitchfork or Hopf bifurcations. The loss ofstability can occur through two qualitatively different mechanisms – through internal instability in the basic latticeelements, or through interaction of the compact breather with the linear passband of the lattice. The former scenariois more typical for high-amplitude breathers, and the latter – for low amplitudes. For the high-amplitude case, insightsinto the nature of compact-mode loss-of-stability are obtained by resorting to the limit of a piecewise-linear system,where interactions are represented by conservative impacts. This issue calls for detailed introspection into integrabilityof piecewise-linear (impacting) systems and their relation to the smooth system. An idea for a sensor based on thestudied mechanisms is suggested. Keywords: flat bands, compact breathers, cubic nonlinearity, resonance, vibro-impact limit, translational invariance
1. Introduction
Following the considerable interest of the recent years instudying breathers (using the terminology and presuppo-sitions introduced in [1]) and solitons in meta-materials,with special emphasis placed on the phenomena of reso-nance and spatial localization, this paper proposes a newmodel system with several advantageous features, and car-ries out theoretical analysis with possible engineering ap-plications.The main object of scientific investigation is a simplifiedmodel for a meta-material – a complex spatially-extendeddynamical system, in which a ‘large-scale’ excitation ap-plied to the boundary can have resonant interaction withspecific degrees-of-freedom. By fine-tuning the ‘structural’parameters, enhanced engineering functionality can be ob-tained, such as acoustic attenuation or shock absorbanceor, alternatively, high sensitivity and amplification of weakexcitation.Special interest is usually placed on optical meta-materials, such as the ones comprised of arrays of op-tical fibers, possibly with the nonlinear dielectric (Kerr)effect. In such systems, usually, identification of spa-tially localized nonlinear periodic solutions (modes), andanalysis of their stability or resonant interaction is per-formed after averaging, using the so-called discrete non-linear Shr¨odinger equation (DNLS), which is itself gener-ally non-integrable (in contrast to some known integrable discrete non-consistent analogues of the continuous non-linear Schr¨odinger equation). The focus on optical meta-materials is justified by a multitude of engineering applica-tions and by the physical-model validity, associated withthe fact that energetic losses in optical fibers are reason-ably negligible.Interesting work has also been done in the context ofacoustic (rather than optical) meta-materials. Althoughacoustic systems, especially ones operating in the nonlin-ear regime, may suffer from relatively prominent energydissipation effects, their study may be well justified. Thefirst reason would be, obviously, related to the sought ap-plication – acoustic noise attenuation, and in some casesalso acoustic sensing, can only be based on acoustic meta-materials. The second reason for one’s interest in acousticmeta-materials is mainly theoretical. It appears that formechanical systems, at least in some simplified geometries,one can exactly solve the equations of motion, obtaining,for example, spatially-localized temporally-periodic solu-tions. Moreover, for certain types of physically justifiableinteraction potentials, the stability problem for the afore-mentioned solutions appears to be tractable, at least par-tially (and not only for the averaged system).Addressing, in specific, the interesting question of spa-tial localization, it should be mentioned that of partic-ular interest would be perfectly-localized periodic solu-tions. Such solutions may be traveling wave envelopes,as in optical applications (the so-called compact solitons,
Preprint submitted to Chaos, Solitons & Fractals August 31, 2020 r compactons), or ‘standing’ perfectly-compact discretebreathers, as observed in acoustic settings. In both cases,perfect-compactness (in contrast to super-exponential butnon-perfect compactness) was shown to be possible to ob-tain when lattices with special local configurational sym-metry are used, ones for which linear spectral analysisreveals the emergence of flat bands (or curves for one-dimensional lattices).The terminology used for the dynamic regime centralfor the study reported in this paper is that of compactbreathers (CBs). The label ‘compact breather’ is used hereto make a distinction from a soliton (self-similar movingnonlinear wave with a localized profile) with perfect com-pactness (the tail-free case), such as the one considered byRosenau in the continuum-limit and named ‘compacton’[2]. In the same time, the relation to a breather owes tothe fact that a discrete medium is considered here and theobtained spatially-localized and temporally-periodic solu-tion is not moving with respect to the medium. Finally, theobserved nonlinear normal mode (NNM) appears strictlyspatially-compact (tail-free), and hence the term compact breather (CB).The study of systems with flat bands (FB) originatedin purely-linear solid-state quantum settings, where flat-band generating systems were studied in many differentgeometries, and the effect of symmetry-breaking fields wasexamined [3]. In subsequent publications, it was foundthat classical systems, whether optical or mechanical, canbe constructed, using the quantum analogy, such that inthe linear setting the system would have a flat band (orseveral flat bands), but then, when augmenting the inter-action potential by anharmonic additions, new noteworthyeffects would be observed.One of the observed effects is the fact that the spatially-detached, perfectly local, periodic mode of the system willstill exist even for anharmonic addition to the potential,given that this addition respects the original symmetriesof the system. Moreover, one can study the stability ofsuch compact localized states (CLSs, see [4]). Numericalstability-results for compact periodic solutions in physi-cally realizable optical models (with cubic ‘Kerr’ nonlin-earity) were given in [5]. A later work showed analyticstability-analysis results for a specific set of parameters fora similar optical setting [6]. Asymptotic stability-analysisresults for such optical systems were subsequently givenfor a family of solutions in [7].Analytic stability-analysis results for s specific set of pa-rameters and asymptotic results for a family of compactsolutions for the acoustic setting (with exact descriptionof the degrees of freedom) were given by the authors in[8]. The limitation of the latter work is the fact that thenonliniarity assumed therein was of the impact type, pro-vided by fixed, translationally non-invariant constraints.Similarly, in the studies dedicated to the optical setting,the model equations (DNLS) are formulated in a movingframework. However, with respect to that framework, notranslational invariance holds, and thus strictly mathemat- ically, the nonlinear dynamic difference equations in whicha compact periodic solution is sought, are not translation-ally invariant.Owing to all of the aforementioned and to the possi-ble applicability of translationally-invariant systems withcompact, marginally-stable, periodic modes in acousticsensor engineering (to give but one example), the presentwork opts to address the problem directly, by exploitingthe piecewise-linear interaction-potential limit, with theanalysis-advantages associated with it.Among the recent works dedicated to the studyof translationally-invariant nonlinear mechanical latticeswith compact solutions, it is worth mentioning [9], wherea prestressed granular one-dimensional chain is examined,and a continuous Nesterenko type equation of motion isconstructed consistently from the discrete lattice limit (forlarge wavelengths). Then, compact moving periodic solu-tions (compact solitons – compactons) are obtained andshown to be stable for certain parameters, both for thebright and the dark soliton cases. It should be noted thatcompactness is obtained not owing to local complexity andsymmetry, but due to taking the continuous limit.In [10], non-integrable DNLS models related to theAblowitz-Ladik equation and Bose-Hubbard systems withpower-law interaction nonlinearitis are examined, withsuper-exponentially localized (descrete) soliton solutionsobserved formally and their stability examined variation-ally.In [11], perfectly localized modes in photonic latticesdescribed by the discrete linear Schr¨odinger equation arestudied for the case of complex arrangements of one-dimensional lattices with one or multiple flat bands withnon-trivial edge effects. Noteworthy is the possibility thatthe model provides for the study of transport in such sys-tems.A very interesting recent study in a mechanical settingis presented in [12], where an array of pendula with beadshaving (generalized) Hertzian interaction (Newton’s cra-dle) is assumed. The slowly modulated amplitude of trav-eling solutions is examined through the associated dis-crete p -Schr¨odinger equation. Traveling waves with super-exponentially decaying tails are observed. The consid-ered system is translationally invariant, and moreover, thevibro-impact limit is addressed, albeit with no explicit re-sults presented regarding analytical stability investigation.In [13], a photonic lattice with two-dimensional Liebstructure and Kerr interaction nonlinearity is studied byaveraging, using the associated DNLS equation. Flatband-related compact discrete (mobile) solitons are observed(numerically) at the zero-power limit.In a very recent work, [14], a two-component Bose-Einstein condensate with cubic short-range interactionsand long-range magnetic coupling is studied, and perfectlycompact (tail-free) accelerating solitons are observed ana-lytically and numerically.The issue of localization in translationally-invariant non-linear chains was considered in a noteworthy earlier work,215], where the Klein-Gordon equation was studied anda rich resonance structure involving nontrivial transmis-sion properties was observed (associated with emergenceof discrete breathers).A very recent paper, [16], addresses the problem ofobtaining a two-dimensional lattice with configurationalsymmetries and translational invariance, characterized byan arbitrary number of flat bands, generated from compactlocalized states (CLSs) occupying an arbitrary number oflattice site. The generating algorithm is derived by solv-ing an inverse eigenvalue problem numerically employingchiral symmetry.Interesting results are presented in [17], where theauthors explore the spectral stability of Kuznetsov-Mabreathers, which are generalized Peregrine solitons. Thestability is analyzed using Floquet theory, where the pe-riod is successively increased until the limiting structure ofPerigrine solitons is reached. In the examined systems, flatdispersion modes are identified and observed to be relatedto compact nonlinear-case solutions with a parametricallyincreasing period. Loss of stability through pitchfork bi-furcation is demonstrated.Another work, [18], performs comparative stability anal-ysis of solitary traveling waves and discrete breathers inFermi-Pasta-Ulam (FPU) and Toda lattices, using energy-based criteria for the Hamiltonian case. Although perfectcompactness is not addressed there, the correspondencebetween stationary discrete breathers and traveling con-tinuous solitons, based on energy considerations, is studiedwithin the context of stability analysis, which is relevantfor the present work.A noteworthy study is reported on in [19], where a spe-cial one-dimensional lattice is constructed, which supportssolutions in the form of (smoothly propagating) ‘tail-free’traveling discrete breathers. The spatial compactness ofthe traveling discrete breathers in the specially-designedlattice arises due to lack of resonance of the localized non-linear mode of the discrete breather with phonon modes,contrary to the standard case reported on in β -FPU lat-tices. It should be noted that the spatial compactnessobserved for the special lattice is not perfect. Rather, in-stead of a constant-amplitude weak tail, which is absentfor continuous compactons (in waveguides with stronglynonlinear dispersion) but does emerge due to lattice dis-creteness in, say, β -FPU lattices, one finds that in thecase of the aforementioned special lattice, the tail is super-exponentially decaying in the co-traveling frame.In [20], a system of two linearly-coupled photonic chainswith cubic on-site nonlinearity is examined. Stationary(immobile) solitons are obtained in the discrete (DNLS)case and are found to be related to Peierls-Nabarro poten-tial pinning. In the continuous limit, the immobile modestransform to solitons, stable for subcritical amplitudes.The central goal of the present work is to give (ana-lytic) insight on the stability of perfectly-compact periodicmodes in conservative translationally-invariant nonlinearmechanical (quasi-)one-dimensional lattices. The structure of the paper is as follows. Sec-tion 2 presents a one-dimensional smoothly-nonlineartranslationally-invariant mechanical system admitting aflat band in the linear regime and compact solutions inits nonlinear extension; Section 3 presents the linear anal-ysis; Section 4 studies the (smoothly) nonlinear regime;Section 5 furthers the analysis of the nonlinear regime byexamining a nonsmooth interaction-potential analogue forthe large-amplitude limit; Section 6 examines a possibleapplication of the theoretical insights gained, and Section7 concludes.
2. The model system
As a model system, we employ a one-dimensional latticewith internal symmetry. The system is sketched in Fig. 1below, with linear properties specified by parameters andnonlinear augmentation denoted by the crossed-out springsymbols.The representative element (unit cell) of the systemcomprises a particle (of mass M ), to which two generallyanharmonic identical oscillators (with mass m each) areattached in parallel. The coupling between the two afore-mentioned oscillators is realized through common bound-ary conditions. For the symmetric mode, corresponding tosynchronous motion of the two oscillators, the system de-generates to a standard dimeric chain of alternating par-ticles of two types, with anharmonic potential. For thecase M = m , a trivial uniform generally nonlinear chainwith identical particles is recovered. The anti-synchronousmode corresponds to the chain decomposing to a collec-tion of non-interacting unit-cells. The local stability of themodes and their interaction, in view of possible traveling-wave (acoustic) perturbations, suggest a rich dynamic pic-ture. To start, linear analysis is performed next.
3. Linear analysis – dispersion bands
In the linear regime, neglecting boundary conditions,the system is assumed to satisfy the following equations ofmotion: M ¨ x n − k u n − x n ) − k v n − x n )++ k x n − u n − ) + k x n − v n − ) = 0 ,m ¨ u n − k x n +1 − u n ) + k u n − x n ) = 0 ,m ¨ v n − k x n +1 − v n ) + k v n − x n ) = 0 (1)Dispersion analysis, that is, assumption of (generally)different-amplitude equal-frequency planar waves for allthree vector displacements x , u , v , yields three dispersionrelations. The first one corresponds to different amplitudesof u and v , and reads (in units of p k/m ):ˆ ω ( q ) = 1 (2)3 ig. 1: Sketch of the (quasi-)one-dimensional (horizontal) model chain where q is the dimensionless wavenumber. This is a flatband . Two additional dispersion relations are obtainedunder the assumption of equal amplitudes of u and v , andthey are as follows:ˆ ω , ( q ) = r µ ± q (1 + µ ) − µ sin q √ µ (3)The following normalization is used:ˆ ω i ( q ) = ω i ( q ) p k/m , µ , M/ m (4)For the special case of µ = 1, one has:ˆ ω µ =12 , ( q ) = p ± cos ( q/
2) (5)For this special case, the three bands intersect at q = π .For µ >
1, there is a frequency gap between the acousticand optical bands. The dispersion relations areˆ ω ( q ) ≤ s µ − | − µ | µ = ( µ − / < , µ > , µ ≤ ω ( q ) ≥ s µ + | − µ | µ = ( , µ > µ − / > , µ ≤ | ∆ˆ ω | ≥ | µ − / − | .For µ <
1, there is intersection between the acoustic andflat bands. This intersection occurs for vanishing groupvelocity. This means that when energy is not transported,it may be transferred to a local mode.Figure 2 shows the dispersion bands for µ > , µ = 1and µ < q ˆ ω ( q ) q ˆ ω ( q ) q ˆ ω ( q ) Fig. 2:
Dispersion bands for the model chain for µ = 1 . µ = 1 (center) and µ = 0 . For the special case of µ = 1, the acoustic branch isnearly dispersionless, with the group velocity ( c g ) decreas-ing from 1 / √ q increases from 0 to π . Generally4ne has c (2) g (cid:12)(cid:12)(cid:12) q =0 = 12 √ µ , c (3) g (cid:12)(cid:12)(cid:12) q =0 = 0 ,c (2 , g (cid:12)(cid:12)(cid:12) q = π = ( ± , µ = 10 , µ = 1 (7)
4. The (smoothly) nonlinear regime
Here we assume potential anharmonicity by adding cu-bic force terms of the Fermi-Pasta-Ulam type. Neglectingboundary conditions, the system is assumed to satisfy thefollowing equations of motion: M ¨ x n − k u n − x n ) − p ( u n − x n ) − k v n − x n ) − p ( v n − x n ) + k x n − u n − ) + p ( x n − u n − ) ++ k x n − v n − ) + p ( x n − v n − ) = 0 ,m ¨ u n − k x n +1 − u n ) − p ( x n +1 − u n ) + k u n − x n )+ p ( u n − x n ) = 0 , m ¨ v n − k x n +1 − v n ) − p ( x n +1 − v n ) + k v n − x n ) + p ( v n − x n ) = 0 (8)(where we added a cubic force coefficient p > M and m ). In order to gain insight on local stability, a single ele-ment is analyzed first. An unconnected element as shownin Fig. 3 is assumed, for the problem to be solvable ana-lytically.
Fig. 3:
Sketch of a representative element of the model chain
The element is assumed to be disconnected, ‘float-ing’, and thus there is conservation of linear momentum(and position). In addition, there is energy conservation.Therefore, the solution lies on a three-dimensional man-ifold and can be described entirely on Poincar´e sections. Moreover, nonlinear normal modes (NNMs) can be iden-tified. For each NNM, the remaining integral can be ob-tained by direct time integration. After one NNM is ob-tained, its linear stability can be established by analysisof a single Hill equation using Hill’s determinants method.The equations of motion are as follows: x = − µ ( u + v ) , m ¨ u + k ( u − x ) + p ( u − x ) = 0 ,m ¨ v + k ( v − x ) + p ( v − x ) = 0 (9)which reduces to a system of two differential equations: m ¨ u + k µ [(2 µ + 1) u + v ] + p µ [(2 µ + 1) u + v ] = 0 m ¨ v + k µ [ u + (2 µ + 1) v ] + p µ [ u + (2 µ + 1) v ] = 0 (10)Introducing the sum of the displacements and their dif-ference, as follows: y , u + v, z , u − v (11)enables rewriting Eqs. 10 as evolution equations for thesymmetric and antisymmetric modes.The equation for the antisymmetric mode is m ¨ z + " k + 34 (cid:18) µµ (cid:19) p y z + 14 p z = 0 (12)The equation for the symmetric mode is m ¨ y + 1 + µµ (cid:18) k + 34 p z (cid:19) y + 14 (cid:18) µµ (cid:19) p y = 0 (13)The motivation of the present work is, in part, to in-vestigate (possibly resonant) interaction between CBs andphonons. Therefore, at this point, we are interested infinding the parameters corresponding to the loss of sta-bility of a CB due to resonance with a propagating wave,which is represented by the symmetric mode. Hence, weneed to obtain the compact antisymmetric mode itself andits stability equation. In order to obtain the antisymmetricmode, we take Eq. 12 to the y → m ¨¯ z + k ¯ z + 14 p ¯ z = 0 (14)where the bar denotes that it is the NNM solution.In order to be able to interpret the single-element anal-ysis results in the context of a chain, the coupling pa-rameters should be identified correctly. In the case of theantisymmetric mode, which is the case we are interested inwithin the context of stability analysis, we have y = x = 0.This means that, within the single element, each of thetwo masses is attached to a stationary point by a sin-gle linear-cubic force, with parameters k and p for each.However, in the chain depicted in Fig. 1, for stationarymasses M , each mass m is attached to stationary points5y two linear-cubic forces with parameters k and p foreach. Ordinarily, forces acting in a queue are not additive.However, since they attach the masses m to two station-ary points, they are equivalent to forces acting in parallel,and thus they are additive. Consequently, for a stable sta-tionary compact mode, the following relations between thepotential parameters can be identified: k = 2( k/
2) = k, p = 2 p (15)Therefore, the antisymmetric-mode equation of motionwould be: m ¨¯ z + k ¯ z + 12 p ¯ z = 0 (16)Introducing the following dimensionless quantities, ω , r km , τ , ω ( t − t | ¯ z =0 ) , () ′ , ddτ ,A , max t { ¯ z ( t ) } , α , pA k , ˆ z , ¯ zA (17)we obtain the dimensionless form of the first integral ofmotion (in Lagrangian description), similar to what is de-rived in [21]: (ˆ z ′ ) + ˆ z + 14 α (ˆ z −
1) = 1 (18)This equation has a solution that can be expressed usingthe fifth Jacobi elliptic function [21]:ˆ z ( τ ) = s α/
41 + α/ (cid:18)p α/ τ (cid:12)(cid:12)(cid:12)(cid:12) α/
21 + α/ (cid:19) (19)Moreover, ˆ z ( τ ) has a Fourier series representation, asfollows: ˆ z (ˆ τ ) = ∞ X n =1 Z n sin ( n ˆ τ ) (20)where ˆ τ = Ω z τ, ˆ m , αα + 2 , Ω z , π p α/ K ( ˆ m ) ,Z n , π s α/
41 + α/ p /αK ( ˆ m ) ( − n − δ n, N − cosh h nπK (1 − ˆ m )2 K ( ˆ m ) i (21)and K ( ˆ m ) is the complete elliptic integral of the first kind(and δ a,b is Kronecker’s delta).The linear stability of the derived antisymmetric modecan be examined by linearizing Eq. (13), which employingEqs. 15 and 17, takes the following form (where lineariza-tion yields the exchange of z by ¯ z ): d y (ˆ τ ) d ˆ τ + h (ˆ τ ) y (ˆ τ ) = 0 ,h (ˆ τ ) , µµ Ω − z (cid:20) α ˆ z (ˆ τ ) (cid:21) (22) (the use of Eq. 15 is justified since for small y , whichis assumed within the linearization framework, the box ofmass M is nearly stationary and the two springs are nearlyequivalent to being connected in parallel).The determination of the Hill function in Eq. 22 requiresthe presentation of the square of the antisymmetric modeas a Fourier series, which becomes:[ˆ z (ˆ τ )] = Z (2)0 ∞ X k =1 Z (2) k cos (2 k ˆ τ ) (23)where Z (2) k , " ∞ X n =1 , , Z n Z k + n − Z k δ k, N − − (1 − δ k, ) ∞ X n =1 , , Z k − n Z k + n δ k, N − (24)Thus, the corresponding Hill function has the followingFourier series representation: h (ˆ τ ) = c + ∞ X k =1 c k cos (2 k ˆ τ ) ,c = 1 + µµ Ω − z + 34 1 + µµ α Ω − z Z (2)0 ,c k = 32 1 + µµ α Ω − z Z (2) k (25)Consequently, we can plot a stability map in theamplitude–mass-ratio plane using Hill’s infinite determi-nants method. For α →
0, Eq. (19) yields a sine and thus:ˆ m → α , Z n → δ n, , Ω z → α , ˆ z (ˆ τ ) → sin ˆ τ ⇒ h (ˆ τ ) → µµ (cid:18) α (cid:19) −
34 1 + µµ α cos (2ˆ τ ) (26)This limit corresponds to the Mathieu equation, forwhich the first instability tongue is absent and only higher-order instability tongues exist, with the following zero-amplitude-limit critical mass ratios: µ α → = 1(1 + N ) − , N ∈ N (27)for which, as expected, µ cr <
1, and the value closest tounity is 1/3.Plotting the dispersion bands for this value shows a fi-nite region of possible intersection between the flat andacoustic bands due to weak nonlinearity, coinciding withthe expectation of parametric resonance as derived above.The dispersion bands plot for the aforementioned mass-ratio is shown in Fig. 4.6 ˆ ω ( q ) Fig. 4:
Dispersion bands for µ = 1 / The stability map for finite amplitudes is shown in Fig.5.
M/m A / ( k / p ) / Fig. 5:
Instability tongues (the first three) for the antisymmet-ric mode for the representative element of the model chain
Additional instability tongues, with decreasing widths,forming an infinite sequence, emerge to the left of the pre-sented range, for
M/m < /
8, at least for √ α <
2. Thedetails of the analysis used to produce Fig. 5 can be foundin [21]. The right boundary of the principal instabilitytongue has a vertical asymptote with µ cr+ → α →∞
1. The leftboundary has a vertical asymptote with µ cr − → α →∞ /
7. Inaddition, the left boundary (of the principal tongue) has alocal maximum with µ ′ ( α ) = 0 at α ≈ . , µ ≈ .
44. For µ >
1, when there is no intersection between the acous-tic and flat bands, the antisymmetric CB is linearly stablefor all amplitudes and could only be destroyed by finite,large-enough effects.The single-element model is valid whenever the anti-symmetric mode is stable, thus the white areas in Fig. 5imply (linear) stability of CBs in the chain. The shadedareas correspond to instability in the single-element model,and they are thus suspected to correspond to instability of CBs in the chain. However, since once instability oc-curs, the single-element model is no longer representativeof the chain, one should observe the behavior of finite per-turbations to CBs in a chain, not unlike when validatingthe implications of linear stability in a general case. Thefollowing section examines the dynamics of small pertur-bations to a single CB solution, assuming, consecutively,local (white) noise, and propagating (phonon or small-amplitude soliton) perturbations.
The equations of motion for a chain in dimensionlessform (where all displacements are relative to the CB am-plitude, A ), in modal coordinates, in terms of previouslydefined quantities, constitute the following dynamical sys-tem: ˆ x ′′ n = ˆ y n + ˆ y n − − x n µ − α µ (cid:8) (2ˆ x n − ˆ y n ) (cid:2) (2ˆ x n − ˆ y n ) + 3ˆ z n (cid:3) ++(2ˆ x n − ˆ y n − ) (cid:2) (2ˆ x n − ˆ y n − ) + 3ˆ z n − (cid:3)(cid:9) , ˆ y ′′ n = − ˆ y n + ˆ x n + ˆ x n +1 ++ α (cid:8) (2ˆ x n − ˆ y n ) (cid:2) (2ˆ x n − ˆ y n ) + 3ˆ z n (cid:3) ++ (2ˆ x n +1 − ˆ y n ) (cid:2) (2ˆ x n +1 − ˆ y n ) + 3ˆ z n (cid:3)(cid:9) , ˆ z ′′ n = − ˆ z n − α (cid:2) (2ˆ x n − ˆ y n ) + (2ˆ x n +1 − ˆ y n ) (cid:3) ˆ z n + 2ˆ z n , ∀ n ∈ N , ≤ n ≤ N − y ′′ N = 2ˆ x N − ˆ y N α (2ˆ x N − ˆ y N )[(2ˆ x N − ˆ y N ) + 3ˆ z N ]4 , ˆ z ′′ N = − ˆ z N − α x N − ˆ y N ) ˆ z N + ˆ z N x ′′ N can also be obtained from Eq.28. The equations for ˆ x ′′ can be replaced by a conservationlaw giving an expression for ˆ x , which for initial conditionsof a CB would be ˆ x = − N P n =2 ˆ x n − µ N P n =1 ˆ y n , yielding the7ollowing two remaining equations for N = 1:ˆ y ′′ = − µ µ ˆ y − N X n =3 ˆ x n − µ N X n =2 ˆ y n + α x − ˆ y )[(2ˆ x − ˆ y ) + 3ˆ z ] − α " N X n =2 (cid:18) x n + 1 µ ˆ y n (cid:19) + 1 + µµ ˆ y ×× " N X n =2 (cid:18) x n + 1 µ ˆ y n (cid:19) + 1 + µµ ˆ y + 3ˆ z , ˆ z ′′ = − ˆ z − α " N X n =2 (cid:18) x n + 1 µ ˆ y n (cid:19) + 1 + µµ ˆ y ++(2ˆ x − ˆ y ) ˆ z + 2ˆ z (30) Writing the state of the chain in vector form, as w =(ˆ x ⊤ , ˆ y ⊤ , ˆ z ⊤ , ˆ x ′⊤ , ˆ y ′⊤ , ˆ z ′ ) ⊤ , the evolution equation for thevariation of the solution from that of a CB at a givenlocation, under the assumption that the variation is small,can be written as: δ w ′ = A ( τ ) δ w (31)where the time-dependent matrix A ( τ ) can be expressedin a block-form as follows: A ( τ ) = (cid:20) N − I N − B ( τ ) N − (cid:21) (32)where the matrix d is a square zeros-matrix of dimension d , N is the number of representative-cells in the chain;the matrix I d is the identity matrix of dimension d , andthe matrix B ( τ ) is the gradient of the right-hand side ofthe equations of motion in Eqs. 28-30, estimated at theCB solution. The full expression for this matrix (for bothfree-ends and periodic boundary conditions) is given inAppendix A.The analysis for the smallest chain with a symmetricsingle CB shows the emergence of resonance overlay, asdepicted in Fig. 6.A CB stability map for a chain of N = 19 sites with free ends and excitation in the central site is given in Fig. 7.Stability regions are scarce. M/m A / ( k / p ) / Fig. 6:
Instability tongues for the (compact) antisymmetricmode for the model chain in a typical rectangular subspace sit-uated around α = 1 , M/m = 1, for N = 3 , N = 2, with pitch-fork instability bounds in blue (online) and Neimark-Sackerinstability bounds in red (online). Fig. 7:
Central-site CB stability map for N = 19 – pitchfork in-stability points denoted by (blue online) dots, Neimark-Sackerinstability points denoted by (red online) ‘+’ symbols. Figure 8 illustrates the stability of a central CB oflarge amplitude ( α = 100) with respect to a small zero-momentum perturbation introduced at the left boundary(by setting x ′ (1) = − y ′ (1)2 µ = 10 − , µ = ).Figure 9 shows the destruction of a central CB of largeamplitude ( α = 100) once it is hit with a small zero-momentum perturbation introduced at the left boundary(by setting x ′ (1) = − y ′ (1)2 µ = 10 − , µ = 1).It is thus demonstrated (in Figs. 8 and 9) that a stablelarge-amplitude CB exists in the considered chain for M = m , and that for other mass ratios, such as M = 2 m , forexample, the CB can be destroyed by phonons.8 z ( t ) -1-0.500.51 Fig. 8:
Central-site CB stability validation for N = 19 , µ = , α = 100 , x ′ (1) | t =0 = − y ′ (1) | t =0 = 10 − . The plots showpropagation maps of (top to bottom) x n ( t ), y n ( t ), z n ( t ), andalso z ( t ) t z ( t ) -1-0.500.51 Fig. 9:
Central-site CB destruction for N = 19 , µ = 1 , α =100 , x ′ (1) | t =0 = − y ′ ( N ) | t =0 = 10 − . The plots show propa-gation maps of (top to bottom) x n ( t ), y n ( t ), z n ( t ), and also z ( t ) .2.4. Compact breather stability – the moderate-amplitude regime Observation of Fig. 7 in a range of moderate amplitudesreveals an interesting picture. It appears that one canhave reasonable prediction of pitchfork instability. Figure10 shows that the boundaries of the two main pitchforkinstability tongues, as calculated for the N = 3 case, ac-curately represent the analogous regions for the N ≫ N = 19) case, implying that the associated instabilitymechanism is local. M/m A / ( k / p ) / Fig. 10:
Central-site CB stability map for N = 19 for mod-erate amplitudes. Pitchfork instability points are denoted by(blue online) circles. Neimark-Sacker instability points are de-noted by (red online) ‘+’ symbols. The dashed curves (cyanonline) bound the zone of resonance between symmetric-modepropagation frequencies and the main CB frequency. Solid(green online) curves represent the boundaries of the two mainpitchfork instability tongues as calculated for the N = 3 case.The dashed black curves at the lower left corner represent theboundaries of resonance of propagation frequencies (‘optical’phonons) with the second CB-frequency.
Moreover, the pitchfork-instability regions lying out-side the two aforementioned tongues seem to be relatedto resonance of the CB frequencies with linear phonons.Those regions seem to form quasi-fractal sets boundedby amplitude–mass-ratio curves corresponding to CB fre-quencies falling within the propagation (optical) range(as calculated in the spectral analysis section). Interest-ingly enough, these instability zones, which are related tothe CB resonance with phonons, seem to be manifestedthrough pitchfork instability. In addition, there seem toemerge Hopf-bifurcation related instability zones that can-not be predicted by linear-spectrum CB analysis. Thestrong dependence of the geometry of the latter on thenumber of elements in a chain (as can be learned from com-parison of Figs. 6 and 10), suggests that these zones canbe associated with resonance between the CB and weakly-nonlinear propagating waves. The dashed curves in Fig. 10 are obtained by comparing1- and 2-multiples of the frequency given in Eq. (21), usingEqs. (4) and (17), with the q = { , π } limits of Eq. 3,taking the ‘+’ sign for the µ ≤ µ − ( α ) = 4 π n h K (cid:16) α/ α +2 (cid:17)i α/ ,µ + ( α ) = π n α/ h K (cid:16) α/ α +2 (cid:17)i − − ; n = 1 , The stability map for large amplitudes and mass ratiosin the range 0 . < M/m < M/m ≈ . M/m ≈ . N = 3 case and the N = 19 case isless than one percent. Therefore, in regard with this first(aforementioned) large-amplitude seemingly-tractable fea-ture, there is no increase in pattern-complexity due to in-crease in the system size. In addition to this quantitativeobservation, it can be noted that even for a smaller sys-tem (that of a single representative element) it was alreadyshown – by rigorous asymptotic analysis of a rigorously de-rived Hill equation – that a sequence of instability-tonguescan be observed. For this sequence, one finds that theboundaries of the emerging instability tongues have ver-tical asymptotes, the quantitative description of which isprovided.This is qualitatively similar to the reality that can beobserved in Fig. 11, with its vertical finite-width stripecorresponding to stability.Therefore, it may be concluded that the first aforemen-tioned phenomenon, the finite gap between the two pitch-fork bifurcation-related instability tongues, appears to betractable.10 / ( p / k ) / M/m M/m
Fig. 11:
Central-site CB stability map for N = 19 for large amplitudes. Pitchfork instability points are denoted by (blue online)dots. Neimark-Sacker instability points are denoted by (red online) ‘+’ symbols (found only around M/m = 0 . M/m = 1)
This result is expected (or, at least, reasonable), giventhat the two instability tongues surrounding the aforemen-tioned gap are the two principal instability tongues, re-lated to the two local modes with which the CB modecan locally exchange stability (in the considered three-on-site-elements chain). In principle, resonance-relatedinstability-tongues are usually characterized by linearly-distinct boundaries. In the natural, α – µ parametriza-tion, these boundaries are indeed linearly distinct. Inthe amplitude–mass-ratio plane, however, one, obviously,observes a square-root-like initial increase in the tongueboundaries, as appears in Fig. 10, hiding the linear dis-tinction (between the boundaries).The second, pitchfork bifurcation-related (and thusmore ordered) distinct feature in the stability map in Fig.11 is an infinitesimal-width stability stripe, situated, as itappears, at M = m . This is a noteworthy phenomenonthat has the potential of being tractable. First of all, thestability at M = m can be observed already in the anal-ysis of the N = 3 case. Moreover, the same argumen-tation as before holds for the boundaries of the principallocal-instability-related tongues and their vertical asymp-totes (this time, however, for the right tongue and itsright boundary, rather than for the left tongue and its leftboundary, as for case of the finite-width stability stripe).However, here, further explanation would be of value,since a zero-width stability-gap is an unusual phenomenon.Usually, collision between boundaries of instability tonguesis related to degeneracies in the system. In the consid-ered case, a degeneracy clearly occurs for M = m , whentwo of the three nonlinear normal modes in the systembecome interchangeable. Such a degeneracy does not nec-essarily have to lead to collision of boundaries, let aloneto an asymptotic one. In fact, it only occurs here forlarge amplitudes, which suggests that it is not only themass symmetry that is at play. Rather, one should sus-pect that the phenomenon is related to combination be-tween mass-symmetry and high energy, or strong nonlin-earity. Strong nonlinearity in the interaction between themasses brings to mind the phenomenon of impact, andindeed one immediately notices that for M = m , impact-interaction between the masses has a special feature ofconserving periodicity of antisymmetric local vibration. One would therefore find it reasonable to consider a limit-case analogy, where nearest-neighbor interaction in thechain is modeled by linear springs connected in parallelto strongly-nonlinear springs. However, those strongly-nonlinear springs would not be smooth ones, modeled bya power-law dependence on displacement, but rather onesof the vibro-impact kind.
5. The vibro-impact limit
The (dimensionless) equation of motion for the antisym-metric compact periodic mode in the smooth system wouldbe: ˆ z ′′ + ˆ z + α z = 0 (34)For a smooth description of a linear system augmentedby a vibro-impact potential (instead of the quartic poten-tial), the equation of motion would have the form [22]ˆ z ′′ + ˆ z + α z α − = 0 (35)where α should be large.One observes that the two equations become identicalfor α = 4, which is just the upper limit of the verticalaxis in Fig. 10. Thus, the general form of the smoothversion of the vibro-impact system is equivalent to thatof the cubic system for α = 4. However, this choice isonly representative of actual impact if one agrees that 4 isa large enough number with respect to unity, which is acertain stretch. Nevertheless, an interesting result can beobtained if such stretched assumption is made.One would therefore opt to explain the “large”-amplitude limit of the stability map in Fig. 7 using the,perhaps, more tractable vibro-impact version of the con-sidered system (having in mind specifically the M = m case).We thus assume a modified system, in which the in-teractions are delivered by coupling consisting of linearsprings attached in parallel to vibro-impact ‘sleeves’, al-lowing maximum relative (normalized) displacement ofunity between the connected masses.11 ig. 12: Sketch of the 1D horizontal chain with linear links and vibro-impact ‘sleeves’ between masses m and M , where themaximum relative free displacement between mass m and any of its neighboring M masses is ∆ α The dispersion bands for this system are the same asfor the first system addressed in the beginning of this pa-per. However, the stability of a CB solution becomes aninteresting problem here, with reasonable inference on thestability of the original system, for the limit of large (butfinite) amplitudes.Figure 11 shows that there seems to be a unique isolated(‘punctured’) point of stability in the large-amplitude limitat M = m (within the examined mass-ratios range). Whatone can do relatively easily is to examine the M = m case,resorting to the system shown in Fig. 12.What should be done first is examination of the dynam-ics of a CB solution located, say, on the central site ofthe chain. The crucial instance is when the vibro-impactpotential comes into action. If for M = m the systemis ‘impact-integrable’, borrowing from a concept discussedin [8], then one can examine the linear stability of the CBusing Floquet theory, employing analytical construction ofthe monodromy matrix. One should start by assuming acentral-site CB, having the mass associated with displace-ment u n labeled as, say, element ‘1’; the mass associatedwith displacement v n labeled as, say, element ‘2’; the massassociated with displacement x n labeled as, say, element‘3’; and the mass associated with displacement x n +1 la-beled as, element ‘4’.It could then be assumed, without loss of generality, thatthe initial condition is u ′ n ( t = 0) = − V , v n ( t = 0) = V − |O ( ǫ ) | , x n ( t = 0) = |O ( ǫ ) | , x n +1 ( t = 0) = |O ( ǫ ) | .In line with this initial condition, one can assume, dueto symmetry, that after elements ‘1’ and ‘2’ pass the dis-tance ∆ α each, a sequence of impacts begins, with the firstimpact occurring between elements ‘1’ and ‘3’. Since themasses of elements ‘1’ and ‘3’ are identical (and even ifthey are only nearly identical), the result of the impactis that element ‘1’ obtains infinitesimal velocity, and ele-ment ‘3’ obtains approximately the velocity − V (that is V ‘leftwards’).The next impact has to involve element ‘2’. There aretwo options. One is impact between ‘2’ and ‘3’, and theother is impact between ‘2’ and ‘4’. In the former case, after impact between ‘2’ and ‘3’, which are two elementsof equal (or approximately equal) mass, moving with (ap-proximately) opposite velocities, the result of the impact isthe reversion of the directions of the two elements. In otherwords, after the second impact, the one between ‘2’ and‘3’, element ‘2’ obtains approximately the velocity − V .The third impact can then occur only between element‘3’, positioned to the left of element ‘1’ and now movingtowards element ‘1’, and element ‘1’, which has infinites-imal velocity. This impact, again, makes the elements ‘1’and ‘3’ swap their velocities (even if only approximately).Consequently, after this third impact, element ‘3’ hasinfinitesimal velocity, while element ‘1’ has approximatelythe velocity V (that is, V ‘rightwards’). Thus, after threeimpacts, the velocities of the two elements of the CB are re-versed and the neighboring elements remain approximatelyat rest. This possibility can be labeled as sequence ‘1-3,2-3,3-1’.The second possibility is that the second impact occursbetween elements ‘2’ and ‘4’. After this impact, the veloc-ity of element ‘2’ becomes infinitesimal, and the velocity ofelement ‘4’ becomes approximately V , due to the approx-imately equal masses of elements ‘2’ and ’4’. The thirdimpact in this version can be either between ‘4’ and ‘1’ orbetween ‘3’ and ‘2’. There is no significance to the order. Itcan be assumed that the impact between ‘2’ and ‘3’ occursfirst. This impact makes ‘3’ and ‘2’ approximately swaptheir velocities. Consequently, element ‘2’ starts movingwith velocity − V , and element ‘3’ resorts back to infinites-imal velocity. Finally, element ‘4’, which after the impactwith element ‘2’ moves with velocity V , now impacts el-ement ‘1’, which has infinitesimal velocity. Consequently,element ‘1’ obtains velocity V (‘rightwards’) and element‘4’ resorts back to infinitesimal velocity. In this option, asbefore, the sequence of four impacts reverses the veloci-ties of the two elements of the CB. This possibility can belabeled ‘1-3,2-4,3-2,4-1’ (or, alternatively, ‘1-3,2-4,4-1,3-2’,which is equivalent).This qualitative analysis shows that for M/m = 1 + O ( ǫ ), the CB with the vibro-impact nonlinearity is indeed12impact-integrable’. Consequently, the monodromy matrixcan be constructed using the saltation matrix. This salta-tion matrix would have to represent both scenarios, onein which there are three consecutive impact within a shortinstance, and the other one, in which there are four consec-utive impacts hidden in the instance of the reversal of theCB velocities. The ‘impact-integrability’ is encompassedin the fact that the result is the same whether the twoelements of the CB experience impact strictly simultane-ously or, rather, consecutively (in each of the two possiblescenarios discussed above), with infinitesimal time lags.The next step is to construct the two saltation matri-ces representing the sequences of impacts. According tothe result obtained in [8], under conditions of ‘impact-integrability’, the total saltation matrix representing theentire sequence of consecutive impacts can be constructedas the product of saltation matrices representing the in-dividual impacts, in the order of their occurrence. Theexpressions for the saltation matrices for the two scenar-ios, taking advantage of the derivations in [8] (computedfor the assumption of perturbations carefully taken to thezero limit), are given in Appendix B (along with some ad-ditional details pertinent to the application of Floquet the-ory to vibro-impact systems). Numerical spectral analysis of the monodromy matrix(analytically-constructed as detailed in Appendix B) forthe two possible impact scenarios shows interesting results.It appears that both for the three-masses impacts andthe four-masses impacts, there is a critical (normalized)frequency corresponding to (subcritical) pitchfork bifurca-tion. For lower frequencies there are two real eigenvalues,one of which is larger than unity. For higher frequenciesall eigenvalues lie on the unit circle. The value of this(first, as shown below) critical normalized frequency is ex-actly the square root of three. This value is obtained bynumerical analysis, in which consecutive digits in the deci-mal representation of this irrational number were obtainedby carefully approaching the bifurcation. This numbercorresponds precisely to the maximum frequency of thepropagation spectrum shown in the linear analysis in thebeginning of the present paper.This correspondence is logical but not trivial. Indeed,one may logically anticipate that for CB frequencies be-low the largest propagation frequency, there may be res-onance between the antisymmetric compact mode and asmall propagating perturbation. This is reasonable. Thenontrivial part is that we see here an example of a sit-uation where the pitchfork bifurcation is related not topurely-local instability, but to something else. Rather, thebifurcation is related to instability corresponding to res-onance between a local mode and an ‘anti-local’ mode.This ‘anti-local’ mode is the planar-wave propagating per-turbation. Quantitatively, this result is summarized in the equation below, where use is made of Eq. 3).ˆ ω (PF)cr = lim µ → / max q { ˆ ω prop ( µ, q ) } == lim µ → / r µµ = s / / √ ω (NS)cr = 2 . N = 19, to comply with the cubic-nonlinearity results).One observes the two real eigenvalues just outside theunit circle in the top plot in Fig. 13.One observes two real eigenvalues just outside the unitcircle in the top plot in Fig. 14, corresponding to the pitch-fork bifurcation in the four-impacts case. Those eigenval-ues move to the unit circle in the plot on the right. Inaddition, in both plots in Fig. 14 there are two pairs ofcomplex conjugate eigenvalues outside the unit circle, illus-trating the Neimark-Sacker-bifurcation-related instability.Those eigenvalues vanish for normalized frequencies higherthan the second critical (normalized) frequency as given inEq. (37). As shown above, pitchfork-bifurcation-related instabil-ity in the vibro-impact system, associated with resonanceof the compact nonlinear mode with the propagation spec-trum, can be avoided for (normalized) frequencies abovethe value of √
3. The same is true for the system with13 FB -1 -0.5 0 0.5 1-0.500.5 FB Re λ Im λ Re λ Im λ Fig. 13:
Monodromy matrix eigenvalues λ for ω/ω FB = 1 . ω/ω FB = 2 . the cubic nonlinearity. Substituting µ = 1 / α corre-sponding to which the principal frequency of the CB for M = m in the cubic-nonlinearity case becomes equal tothe upper limit of the (normalized) propagation frequen-cies spectrum, namely, √
3. Solving this equation for α ,one obtains the value α c = 5 .
49. This value is about 30to 35 percents larger than the value 4, associated with theequivalence between the cubic system and a smooth systemrepresenting a vibro-impact potential under the assump-tion that 4 is a large-enough number. For our purposes,an error of 30 to 35 percents is tolerable, being about halfan order of discrepancy.Next, we can write the energy of the system (as it movesin the compact mode) for the cubic nonlinearity, as follows: E cu = 12 kA (cid:16) α (cid:17) (38)On the other hand, for the vibro-impact limit, the an-tisymmetric mode in compact periodic dynamics admitsdisplacement which is sinusoidal in time until first impact.The associated time between impacts is related to inversenormalized frequency. This allows one to express the en-ergy of the system as follows (assuming that the maximum -2 -1 0 1-1-0.500.51 FB -1 -0.5 0 0.5 1-0.500.5 FB Re λ Im λ Re λ Im λ Fig. 14:
Monodromy matrix eigenvalues λ for ω/ω FB = 1 . ω/ω FB = 2 . displacement before impact is equal to 1 in units normal-ized the same way as for the cubic system, that is by theamplitude, A ): E v-i = 12 kA (cid:0) π ˆ ω − (cid:1) (39)Requiring energy equivalence between the cubic andthe vibro-impact systems and substituting the critical(squared normalized amplitude) value α = 5 .
49, one ob-tains the critical normalized frequency above which reso-nance of the compact mode with the propagation spectrumis avoided – if drawing from the (somewhat stretched)equivalence between the cubic and the vibro-impact non-linearity types: ˆ ω v-i/cu c = 2 .
22 (40)This is the value above which resonance with the propa-gation spectrum can be avoided both for the cubic and forthe vibro-impact system, if one requires their equivalence.In contrast, if the two cases are treated independently,the corresponding value of the normalized frequency is √ ≈ . α = 4 ≫
1, which leads to the 30 to 35 percents14rror in terms of α c , also leads to about 30 percents of er-ror (0.2836, to be precise) in terms of ˆ ω c . This error is notnegligible, but it is also not very large, only half an order,meaning that the calculations are more correct than not,and that some of the essence of the equivalence is indeedcaptured. Clearly, this error could have been smaller hadwe solved a system with a quintic force-term added to thelinear interaction, instead of the cubic term. In that case,instead of the elliptic Jacobi function, a different functionwould have emerged, with a different frequency-amplituderelation, and the critical amplitude for avoiding resonancewith the propagation spectrum would have had a differentvalue. This value, used in the energy equivalence relation,would have yielded back, most probably, a value of thecritical frequency much closer to √
3. But the essence ofthe matter is qualitatively captured already in the ana-lyzed case. The crucial point in drawing the equivalencewith the vibro-impact system is the possibility to explainthe isolated stability stripe around M = m in the param-eter plane. This isolated stability stripe has been shownhere to emerge from considerations of impact-integrabilityfor the implicitly-related nonsmooth limit.In order to conclude the analysis, it would be valuableto look whether the smooth system has additional zero-width vertical stripes of stability in the parameter plane,for values of M/m other than 1, at least in the consideredrange
M/m ∈ [0 . , α ≫ , M = m region isexclusively stable within its neighborhood. The α ≫ high-energy limit in the smoothsystem, as can be learned from Eq. (38). Without havingformal equivalence in the equations of motion for α ≫
1, itmay still be interesting to address the question of stabilityfor the M = m case in the high-energy limit for the vibro-impact system. From Eq. (39) we see that the only wayto obtain high energy in the vibro-impact system is takingthe ˆ ω ≫ energy scales as ˆ ω . Thus,the limit in which the vibro-impact system is energetically-equivalent to a high-energy smooth system for M = m (forthe case where the compact antisymmetric periodic solu-tion is stable), is the high-frequency limit. As previouslydiscussed, above the derived critical frequency, the vibro-impact system for M = m does indeed have its compactperiodic antisymmetric solution linearly stable. In thissense, the large-amplitude-limit stable stripe at M = m in the smooth system, associated with the compact peri-odic antisymmetric mode, can be considered qualitativelyexplained – with the help of the nonsmooth-limit analysis. The existence of additional zero-width stability stripesin the parameter plane for the case of smooth nonlinearityis addressed here by referring to the related vibro-impactsystem. For the latter system, the parameter plane rele-vant for the stability diagram is represented by the coor-dinates µ and ˆ ω . In order for stable regions, as wouldemerge from application of Floquet theory, to exist inthis plane, first, there need to exist µ -values guaranteeingimpact-integrability. It was already shown that impact-integrability is guaranteed for M/m = 1. The questionis, are there additional
M/m values in the range [0 . , m positioned between twomasses M ).To answer this, we, as before, denote the momenta ofthe ‘upper’ mass m , the ‘lower’ mass m , the ‘left’ mass M and the ‘right’ mass M by p , p , p and p , respectively.Next, it is assumed that one can express p in terms ofthe other three momenta using conservation of total linearmomentum for the four masses. Then, expressing p (andthus also p ) in terms of p and p using conservation ofenergy, we can obtain a two-dimensional mapping for thestate vector p = { p , p } ⊤ . This mapping can representall the possible impacts occurring between the four masses.To this end, the mapping would have to be defined bycases, with a different set of µ -dependent coefficients foreach of the four possible impacts (1-3,3-2,2-4,4-1). Clearly,this mapping would be nonlinear, due to the nonlinearityof the dependence of the Hamiltonian on the momenta.Now, even the one-dimensional nonlinear mapping canyield chaotic trajectories, let alone a two-dimensional one.Therefore, whether the exact-inversion mapping appliedto the couplet {− , } ⊤ exists for a specific value of µ ,would depend on the number of impacts occurring, andtheir order. This is not an easy problem to solve, andhad we wanted to find a positive answer, we would havebeen forced to go the hard way of direct enumeration ofthe possible sequences. For each value of µ , starting from p = { p , p } ⊤ = {− , } ⊤ , we would have been forcedto check all possible combinations of the feasible impacts,examining the final value of the vector p in each case andlooking for the results that yield { , − } ⊤ .Luckily, we do not have to follow this difficult route.As numerical analysis of the smooth system shows, thereappear to be no additional impact-integrable µ -values inthe examined range. This means that if the vibro-impactlimit is any indicator, the instantaneous impact-dynamicsof the four masses for the compact periodic antisymmetricmode is not impact-integrable. This means that not all possible trajectories in the phase-space of this system areperiodic. This implies that at least some trajectories arenon-periodic. This, in turn, means that it would be suf-ficient to find a non-periodic trajectory in a sub-space ofthe phase-space of the four-mass system, for µ in the ex-amined range, to prove the absence of additional isolated15ero-width stable stripes in the parameter plane. Thissub-space could be associated with specific initial condi-tions manifested as a specific infinitesimal perturbation ofthe antisymmetric compact mode.A reasonable choice would be perturbation of the initialconditions that would lead to a sequence of impacts in-volving not all four but only three masses. Then, insteadof four different feasible impacts (on the order of occur-rence of which one has to wonder), there would be onlytwo ‘types’ of impact, occurring one after the other.This is a reasonable scenario. If we assume that theinitial leftward velocity of particle 1 is slightly larger thanthe rightward velocity of particle 2, then the first impactwould be between particles 1 and 3. Then, not only for thecase of M = m , as mentioned earlier, but, in fact, for anymass ratio, the second impact would be between particles3 and 2.The reason for this is that after the first impact, thedistances between the impacting parts of particles 2 and 4and those of particles 2 and 3 are the same. However, thespeed of approach of particles 3 and 2 is larger than thatof particles 2 and 4. The latter fact owes to the negativevelocity that particle 3 would have after the first impact(for any mass ratio).After the impact between particles 3 and 2, the speed ofparticle 2 would decrease. If the decrease is sufficient, im-pact with particle 4 could be avoided. In addition, particle4 has to have initial negative infinitesimal displacement,large enough for impact between particles 1 and 4 to beavoided. In the same time, that displacement should besmall enough for impact between particles 2 and 4 to beavoided.In the following, explicit kinetic and kinematic calcula-tions are performed for the case of perturbation in the ini-tial conditions which guarantees the occurrence of impactsbetween particles 1, 2 and 3 alone (with no participationof particle 4 in the impacts-sequence). The following perturbed antisymmetric initial condi-tions just before the first impact are assumed for the 4particles: u (0) n = − ∆ , ˙ u (0) n = − v [1 + O ( ǫ )] , v (0) n = ∆(1 − ǫ ) , ˙ v (0) n = v , x (0) n = 0 , ˙ x (0) n = 0 , x (0) n +1 = − ζ, ˙ x (0) n +1 = 0 (41)Since the sequence of impacts is quasi-instantaneous,the coordinates do not change during the process, and thekinetics is fully described by the momenta.Assuming ideal impacts, one can use conservation oflinear momentum and energy for any two impacts (sincein the perturbed antisymmetric case there are only im-pacts involving two particles, simultaneously), to obtainthe transformation mapping for any two-dimensional ve- locity vector before and after impact: v ′ , = 1 − ˆ µ µ v , + 2ˆ µ µ v ,v ′ = 21 + ˆ µ v , − − ˆ µ µ v , ˆ µ , Mm (42)here the notation is changed and all the velocities are de-noted by v , with the subscript referring to the particleindex (unless it is zero, in which case the value refers tothe nominal speed of a particle of mass m of the CB justbefore impact). The subscript in the expressions belowrefers to the serial number of the latest impact that hadalready occurred.The relations given above are true for impacts betweenparticle 3 (having mass M ) and either particle 1 or 2 (eachhaving mass m ). Particle 4 (with mass M ) is assumed notto interact by impact, and this assumption is validated inthe following subsection, dedicated to the kinematics.Having the initial conditions (just before the first im-pact) and the transformation relations given above, wecan ‘integrate’ to obtain the velocities after impact. Do-ing this we bear in mind that the impact dynamics isquasi-instantaneous, in the sense that the correspondingdisplacements are infinitesimal. Consequently, linear elas-tic coupling (the springs) has no effect on the dynamics.As aforementioned, the first impact is between particles 1and 3. Therefore, the velocities of those particles after im-pact become (neglecting ǫ with respect to 1 in the kineticcalculations here and onward): v (1)1 = − − ˆ µ µ v , v (1)2 = v , v (1)3 = −
21 + ˆ µ v (43)For 0 < ζ ≪ , ˆ µ >
1, particle 1 changes the directionof its velocity before experiencing impact with particle 4.Consequently, particle 1 starts moving rightwards, and atthe considered quasi-instance it will not experience impact(through the ‘sleeve’, as suggested in the sketch in Fig. 12)with particle 4. Therefore, we have shown that one canavoid impact between particles 1 and 4 by proper choice ofinitial perturbations, for the case of ˆ µ >
1. After the firstimpact, for ˆ µ >
1, particle 1 is still moving leftwards. Thesecond impact involves particles 2 and 3. The velocitiesafter this second impact would be as follows: v (2)1 = − − ˆ µ µ v ,v (2)2 = 1 − µ − ˆ µ (1 + ˆ µ ) v ,v (2)3 = 4(1 + ˆ µ ) v (44)As far as the kinetics of particle 2 is concerned, it is evi-dent that impact between particles 2 and 4 can be avoidedat the considered quasi-instance if the velocity of particle2 changes its direction. Negative velocity can be obtainedfor particle 2 after its first impact with particle 3 under the16ondition ˆ µ = M/m > √ − ≈ . M/m ∈ [0 . ,
2] is thatof the velocity of particle 1. It is clear from Eq. (44) thatparticle 3 will move toward particle 1 after the impact withparticle 2, and thus the third impact will involve particles1 and 3, and produce velocities as follows: v (3)1 = − µ + ˆ µ − ˆ µ (1 + ˆ µ ) v ,v (3)2 = 1 − µ − ˆ µ (1 + ˆ µ ) v ,v (3)3 = − − ˆ µ )(3 + ˆ µ )(1 + ˆ µ ) v (45)Solution of the cubic equation in the numerator inthe first expression given above shows that direction-reversal of the velocity of the first particle after its sec-ond impact with particle 3 is obtained for ˆ µ = M/m ∈ [0 . , . µ = M/m ∈ [0 . , µ = M/m ∈ [0 . , . µ = M/m ∈ [0 . , . v (4)1 = − µ + ˆ µ − ˆ µ (1 + ˆ µ ) v ,v (4)2 = (1 − ˆ µ )(1 − µ − µ − ˆ µ )(1 + ˆ µ ) v ,v (4)3 = 8 1 − µ − ˆ µ (1 + ˆ µ ) v (46)For ˆ µ < − µ − µ − ˆ µ .This cubic function is negative for ˆ µ > . µ = M/m ∈ [0 . , . µ = 1, exactly three impact (in total) reverse the initial velocities of particles 1 and 2 and bring the velocityof particle 3 back to zero.It has thus been shown that in the parameter range an-alyzed in the smooth case, small perturbation of initial an-tisymmetric compact conditions exists, for which particle4 does not participate in the quasi-instantaneous impact-dynamics – at least from the kinetics perspective (the ve-locities involved). The next subsection examines the con-sistency of the aforementioned assumption of no impactsinvolving particle 4 from the perspective of the kinematics– the positions of the particles. We assume, with no loss of generality, and remainingaccurate up to a second order correction in ǫ , that one canstart with a relaxed system at rest and supply velocities v (0) and − v (0) + O ( ǫ ) to particles 1 and 2, respectively,at positions q (0) = q (0) = 0. Consequently, after cer-tain time t , particle 1 will reach the position − ∆, andparticle 2 will reach the position ∆(1 − ǫ ). The velocitiesof the particles just before the first impact would be − v (of non-specified value) and v − O ( ǫ ), for particles 1 and2, respectively. In this framework, the first impact occursbetween particles 1 and 3. The requirement for avoidingimpact between particles 1 and 4 up to this instance is sim-ply − q (0) = ζ >
0. The second impact occurs betweenparticles 2 and 3. The time from the first impact to thesecond impact is the time it takes the right boundary ofthe sleeve of particle 3 to cover its distance from particle2, moving in the velocity of approach between them: δt = ∆ ǫv (1)2 − v (1)3 = 1 + ˆ µ µ ∆ v ǫ + O ( ǫ ) (47)During this time, the position of particle 1 changes to q (2)1 = − ∆ − v (1)1 δt = − ∆ − − ˆ µ µ ∆ ǫ + O ( ǫ ) (48)This means that up to this instance (or always if ˆ µ > ζ = − ˆ µ µ ∆ ǫ + ∆ ǫ / . At the time of the firstimpact with particle 3, the position of particle 2 is q (2)2 = ∆ − ∆ ǫ + 1 + ˆ µ µ ∆ ǫ + O ( ǫ ) (49)Taking into account the introduced displacement of par-ticle 4 needed to guarantee avoidance of impact with par-ticle 1, the requirement for no impact between particles 2and 4, at this instance, is guaranteed by definition (ˆ µ > µ ∈ [0 . , µ ∈ [0 . , − δt = q (2)1 − q (2)3 v (2)3 − v (2)1 = (1 + ˆ µ ) (3 + ˆ µ )(5 − ˆ µ ) ∆ v ǫ + O ( ǫ ) (50)From the first impact of particles 2 and 3 till the secondimpact of particles 1 and 3, the time interval given aboveallows particle 1 to shift further to the left (for ˆ µ < q (3)1 = q (2)1 + v (2)1 δt = − ∆ − − ˆ µ µ ∆ ǫ − − ˆ µ µ (1 + ˆ µ ) − ˆ µ ∆ ǫ + O ( ǫ ) (51)For self-consistency of the assumption of no impact be-tween particles 1 and 4 for ˆ µ <
1, the absolute value ofthe displacement of particle 4 just before the first (overall)impact, has to satisfy ζ = 2 − µ − ˆ µ ∆ ǫ + ∆ ǫ / (52)The reason for the presence of the second, fractional-order, term is that there has to be a positive addition oforder higher than unity to maintain the tightest possiblebound and allow the maximum freedom for particle 2. Inthe same time, it is necessary to have the additional termstronger than (possible) second-order corrections, emerg-ing from expansion of rational functions originating fromlinear corrections to the initial conditions for the veloci-ties. The chosen power is simply a symmetric compromisefor the exponent, and the coefficient of unity is a symmet-ric positive choice of a coefficient of order unity with nomagnitude-significance.The last validation that has to be made now concernsthe position of particle 2 at the time of its second impactwith particle 3. This position can be calculated accordingto formulas similar to the ones used in previous calcula-tions. First, one needs to calculate the position of particle2 at the instance of the third (overall) impact: q (3)2 = q (2)2 + v (2)2 δt = ∆ − ∆ ǫ + 1 + ˆ µ µ ∆ ǫ ++ 1 + ˆ µ µ − µ − ˆ µ − ˆ µ ∆ ǫ + O ( ǫ ) (53)Next, the time interval between the third and the fourthimpact should be derived, using the same strategy as be-fore: δt = q (3)1 + 2∆ − q (3)2 v (3)2 − v (3)3 =1 + ˆ µ µ µ + ˆ µ − ˆ µ (1 + ˆ µ ) − µ − µ − ˆ µ ∆ ǫv + O ( ǫ ) (54) A positive value for the time interval is guaranteed for7 − µ − µ − ˆ µ >
0, which holds for ˆ µ ∈ [0 , . second impact between particles 2 and 3 necessaryfor reversing the direction of particle 2 (the range ˆ µ ∈ [0 . , . q (4)2 = q (3)2 + v (3)2 δt =∆ − ∆ ǫ + 1 + ˆ µ µ ∆ ǫ + 1 + ˆ µ µ − µ − ˆ µ − ˆ µ ∆ ǫ ++ 1 + ˆ µ µ − µ − ˆ µ − ˆ µ µ + 5ˆ µ + ˆ µ − µ − µ − ˆ µ ∆ ǫ + O ( ǫ ) (55)The self-consistency requirement for the absence of im-pact between particles 2 and 4 can then be expressed bythe inequality q (4)2 + ζ < ∆. This leads to the followingrequirement: 1 + ˆ µ µ + 1 + ˆ µ µ − µ − ˆ µ − ˆ µ +1 + ˆ µ µ − µ − ˆ µ − ˆ µ µ + 5ˆ µ + ˆ µ − µ − µ − ˆ µ + 2 − µ − ˆ µ < µ ∈ [0 . , . second impact between particles 2 and 3, necessaryfor reversing the direction of particle 2. In addition, one re-calls the aforementioned requirement 7 − µ − µ − ˆ µ > − ˆ µ >
0. All these inequalities,when superimposed, lead to the following equivalent self-consistency condition necessary for avoiding impact be-tween particles 2 and 4: ˆ µ < µ ∈ [0 . , . finite margin , clearlythe ǫ / correction in ζ is insufficient for enabling contactbetween particles 2 and 4. Therefore, it is established thatthere exists an infinitesimal perturbation of compact an-tisymmetric initial conditions for the vibro-impact-limitsystem, for which the dynamics involves only particles 1,2, and 3. This dynamics consists in possibly an infinite se-quence of consecutive impacts of particle 3 with particles1 and 2, eventually leading (or not) to exact reversal ofthe velocities of the particles assigned to them just beforethe first impact. If the aforementioned sequence consistsof a finite number of steps, then the 3-element subsys-tem can be considered impact-integrable. Otherwise, thesubsystem can be considered not impact-integrable. Now,if a slightly-perturbed periodic solution is non-integrable,even for a specific (feasible) perturbation, then the exactperiodic solution in question is, clearly, unstable.18hat is now left to do in order to answer the stabilityquestion for the compact periodic solution definitively, forthe assumed parameter range, is to find the parametervalues corresponding to impact-integrability of the three-body problem (for the examined system). For those massratios, a perturbation would have to be found that wouldlead to four-body dynamics. For this four-body dynamics,it would have to be shown that there is either no impact-integrability, or no frequency guaranteeing the stability ofthe compact periodic solution (as should be learned bysaltation matrix analysis). It is clear that perturbed compact antisymmetric ini-tial conditions involving three of the masses cannot leadto quasi-instantaneous velocities-reversal with a single im-pact. The smallest number of sufficient impacts would be2. From Eq. (44), the triplet { v (2)1 , v (2)2 , v (2)3 } becomes thetriplet { v , − v , } for ˆ µ → ∞ . The next natural numberof impacts to yield velocities-reversal is 3, which, as canbe found from Eq. (45), is obtained for ˆ µ = 1.In order to find the values of ˆ µ corresponding tovelocities-reversal for 4,5,6 etc. impacts, the three-bodyproblem needs to be further integrated, continuing the se-quence presented in Eqs. (43)-(46).First, one observes that the triplet { v (4)1 , v (4)2 , v (4)3 } be-comes the triplet { v , − v , } for ˆ µ = √ − ≈ . µ ∈ [0 . , after the fifth impact): v (5)1 = − µ − µ − µ − ˆ µ + ˆ µ (1 + ˆ µ ) v ,v (5)2 = (1 − ˆ µ )(1 − µ − µ − ˆ µ )(1 + ˆ µ ) v ,v (5)3 = − − µ − µ + 4ˆ µ + ˆ µ (1 + ˆ µ ) v (58)The triplet { v (5)1 , v (5)2 , v (5)3 } becomes the triplet { v , − v , } for ˆ µ = √ − ≈ . µ ∈ [0 . , after the sixth impact): v (6)1 = − µ − µ − µ − ˆ µ + ˆ µ (1 + ˆ µ ) v ,v (6)2 = 1 − µ + 85ˆ µ + 48ˆ µ (1 + ˆ µ ) v − µ + 12ˆ µ + ˆ µ (1 + ˆ µ ) v ,v (6)3 = 4 (1 − ˆ µ )(3 − µ − µ − µ )(1 + ˆ µ ) v (59) The triplet { v (6)1 , v (6)2 , v (6)3 } becomes the triplet { v , − v , } for ˆ µ = √ − ≈ . µ ∈ [0 . , v (7)1 = − µ − µ − µ (1 + ˆ µ ) v ++ 125ˆ µ + 51ˆ µ + ˆ µ − ˆ µ (1 + ˆ µ ) v ,v (7)2 = 1 − µ + 85ˆ µ + 48ˆ µ (1 + ˆ µ ) v − µ + 12ˆ µ + ˆ µ (1 + ˆ µ ) v ,v (7)3 = − − µ + 49ˆ µ + 76ˆ µ (1 + ˆ µ ) v − µ − µ − ˆ µ (1 + ˆ µ ) v (60)It turns out that the triplet { v (7)1 , v (7)2 , v (7)3 } becomes thetriplet { v , − v , } for no real positive value of ˆ µ .Finally, in order to see whether velocities-reversal can beobtained for exactly 8 impacts, we apply Eq. (42) to Eq.(60) for the eighth impact, which would occur betweenparticles 3 and 2. The following velocities are obtained( after the eighth impact): v (8)1 = − µ − µ − µ (1 + ˆ µ ) v ++ 125ˆ µ + 51ˆ µ + ˆ µ − ˆ µ (1 + ˆ µ ) v ,v (8)2 = 1 − µ + 364ˆ µ − µ − µ (1 + ˆ µ ) v ++ − µ + 44ˆ µ + 16ˆ µ + ˆ µ (1 + ˆ µ ) v ,v (8)3 = 16 1 − µ + 21ˆ µ + 20ˆ µ (1 + ˆ µ ) v −
16 5ˆ µ + 6ˆ µ + ˆ µ (1 + ˆ µ ) v (61)The triplet { v (8)1 , v (8)2 , v (8)3 } becomes the triplet { v , − v , } for ˆ µ = p − √ − ≈ .
08. This value isalready outside the examined range of ˆ µ ∈ [0 . , µ j , where j is the number of impacts requiredfor exact velocities-reversal and ˆ µ j is the corresponding19ass ratio:ˆ µ → ∞ , ˆ µ = 1 , ˆ µ = √ − ≈ . , ˆ µ = √ − ≈ . , ˆ µ = 2 √ − ≈ . , ˆ µ = q − √ − ≈ .
08 (62)It would be reasonable to conclude by informal induc-tion that ˆ µ j is (at least weakly) monotonically decreasingwith j when sufficiently small. We can also address thismatter more rigorously asymptotically . Performing first-order expansion of Eq. (42) for particle 1 for ˆ µ ≪
1, onegets: v ( k +2)1 = v ( k )1 + 2ˆ µ h v ( k +1)3 − v ( k )1 i (63)where k denotes the impact number. Formally, one can‘integrate’ this difference equation starting from the initialcondition, which to first order is v (0)1 = − v , up to thevelocity after the last impact between particles 1 and 3,which for impact-integrability should be v ( j )1 = v . Thisleads to the following relation: v = ˆ µ ⌊ ( j − / ⌋ X k =1 , , h v ( k +1)3 − v ( k )1 i == ˆ µ ⌊ ( j − / ⌋ X k =1 , , v ( k +1)3 − ˆ µ ⌊ ( j − / ⌋ X k =1 , , v ( k )1 (64)The sequence v ( k )1 goes from − v to v , passing throughzero. Moreover, we look for values of ˆ µ corresponding toimpact-integrability. Therefore, all the assumed sequencesare finite and periodic, and thus time-reversal symmetryapplies, and hence v ( k )1 is also symmetric on a period or antisymmetric on half a period (which is what we con-sider in the velocity-reversal problem). Consequently, thesequence v ( k )1 is antisymmetric with respect to the “mid-dle” impact. Therefore, the last sum in Eq. (64) amountsto zero for large-enough values of j , for which the sequence v ( k )1 is smooth enough for the sum to approximate an (odd-function) integral.Regarding, v ( k )3 , we can say the following. For ˆ µ → k , one has v ( k )3 = 4( k/ v . As we consider periodicimpact-sequences, one can assume that for large-enoughvalues of j , for k close-enough to j , symmetric behaviorwould emerge, namely: v ( k )3 = 4[( j − k ) / v . For non-zero values of ˆ µ , v ( k )3 should initially increase with each even impact. Then, when v ( k )1 and v ( k )2 would reach theirintermediate values close to zero, all the energy would bestored in v ( k )3 , and thus, in the middle of the sequence, oneshould expect a maximum value of v ( k )3 , after which theenergy will again ‘flow’ to v ( k )1 and v ( k )2 , and v ( k )3 will returnto its (initial) smallest value. This implies that there is amaximum value of v ( k )3 and that one can estimate it from energy considerations, namely, assuming that particle 3contains all (or almost all) the energy: v (max)3 ∼ v p ˆ µ/ v and v only slightly).Substituting this result, along with the understandingregarding the sum over v ( k )1 into Eq. (64), one obtains thefollowing relation: v ∼ ˆ µ v p ˆ µ/ ⌊ ( j − / ⌋ X k =1 , , v ( k +1)3 v (max)3 ∼∼ v p µ (cid:22) j − (cid:23) ρ, ρ , h v ( | k even)3 i avg v (max)3 (66)In order to estimate the last quotient in Eq. (66), onehas to know the distribution v ( k +1)3 for odd increasing val-ues of k (the ‘current’ impact number). For a large numberof impacts needed for completing half a period, a smoothapproximation for this distribution would be a concavefunction, linear at the edges (initial and final values of k ), and symmetric around the middle. In principle, forsuch a distribution, the value of the aforementioned quo-tient should be between 1/2 (for a triangular function) and1 (for a nearly uniform distribution – except the sharpchange at the edges of the range). Any convexity in thedistribution would imply increase in the forcing, whichis unreasonable, since the force is exerted by particles 1and 2, whose momenta decrease from the beginning of therange to its middle. Thus, the feasible range would be ρ ∈ [1 / , µ ≪
1, in-sight into the problem is given in Appendix C. It is shownthere that a good estimate for the velocity distribution ofparticle 3 at even-valued impact numbers is the sinusoidal distribution on the positive half period. Consequently,the corresponding aforementioned quotient can be takenas ρ ∼ /π .Using this result and assuming large values of j (manyimpacts for velocity reversal – the assumption of large j is crucial for the calculation of the estimate for ρ ), oneobtains the following asymptotic relation:ˆ µ ∼ j ≫ π j (67)which implies inversely quadratic decay of ˆ µ with j .Extrapolating back to 8 impacts, the highest number ofimpacts for which we have exact integration results, oneobtains from Eq. (67) the following estimate:ˆ µ ∼ π ≈ .
08 (68)20hich is fairly close to the exact value given in Eq. (62).This gives validation to the relevance of the asymptoticanalysis presented herein.Hence, one can conclude that, indeed, for j >
8, one hasˆ µ j < ˆ µ , and thus the only ‘impact-integrable’ values ofˆ µ falling into the chosen range of [0 . ,
2] are the 4 valueslisted in Eq. (62).It has thus been shown that there are only 4 (isolated)points (mass-ratio values) that should be suspected asstable in the discrete-limit system. If the analogy withthe smooth (linear-cubic) system is worthwhile, as arguedabove, then there are only 4 points ‘suspected’ to corre-spond to ‘isolated’ stability of the compact antisymmetricsolution in the smooth system for large amplitudes. Oneof those points, namely, the one for which M = m , wasalready shown to be stable.Now, instead of scanning the entire parameter domain,one can assume a value of the amplitude at the upperbound of the examined range, namely α = 400, take N =19 representative elements, and perform numerical Floquetanalysis only for the cases of ˆ µ , M/m = ˆ µ , ˆ µ , ˆ µ , ˆ µ .The case of M/m = ˆ µ = 1 is already known to correspondto stability (for high-enough frequencies), but its resultsare shown in the figures below for clarity anyway, alongwith the results for the other three mass-ratio values.Clearly, it could have been advantageous to examinefour-particle dynamics for the four relevant mass-ratio val-ues, however, this task may be rather cumbersome. It issufficient, though, that the discrete-system analogy dis-tilled a small number of isolated mass-ratio values to ex-amine in the smooth large-amplitude case. When thereare specific values to check the stability for in the param-eter plane, the tendency of the stable domain to be quasi-fractal no longer poses an obstacle for analysis, whichrenders the problem tractable. The eigenvalues of themonodromy matrices of the compact periodic solutions for α = 400 , N = 19 and the aforementioned mass ratios (ob-tained by numerical implementation of Floquet theory forthe smooth system) are presented in Figs. 15 and 16.One notes that of the four mass ratios, stability of thecompact antisymmetric solution is observed for an isolatedvalue of the mass ratio only for M/m = 1, as suggested bycruder analysis shown already in Fig. 11.To conclude, it may be asserted that it was possible toexplain the stability pattern of the compact antisymmet-ric periodic dynamic solution in a translationally-invariantnonlinear one-dimensional lattice, at least in what con-cerns stability in the context of pitchfork bifurcations. Itis not surprising that pitchfork bifurcations appear to betractable in the analyzed system, in contrast to Neimark-Sacker bifurcations. The reason for the tractability canbe the fact that pitchfork bifurcations are normally asso-ciated with a symmetry group of higher symmetry, andhigher-symmetry systems (modes, in this case) are nat-urally more tractable (this is, of course, the reason forthe power of symmetry considerations when obtaining so-lutions in physics in general). Indeed, the two instabil- Re λ Im λ Re λ Im λ Fig. 15:
Monodromy matrix eigenvalues λ for the compact an-tisymmetric periodic mode for N = 19 , α = 400 for M/m = 1(top) and
M/m = √ − smooth system ity cases that were explained here were related to reso-nance between the compact antisymmetric nonlinear nor-mal mode and either phonons (least-localized modes), orthe (complementary) symmetric nonlinear normal mode,which is (exponentially) localized. In both cases, the cor-responding bifurcation was of the pitchfork type. Boththe strong (exponential) and the vanishing (planar-wave)localization are examples of extreme levels of localization,and as such, they can be reasonably-well argued to be as-sociated with high (spatial) symmetry.Thus, the proposed translationally-invariant nonlinearone-dimensional mechanical lattice shown here to possessstable compact modes can be considered well-enough stud-ied, at least for a first look into such a system. Higher-order effects can be addressed in future work.The following section suggests a possible application forthe studied phenomena.21 Re λ Im λ Re λ Im λ Fig. 16:
Monodromy matrix eigenvalues λ for the compact an-tisymmetric periodic mode for N = 19 , α = 400 for M/m = √ − Mm = √ − smooth system
6. Application – a moving acoustic sensing devise
A natural application of the studied system could be inthe field of acoustic sensor engineering. It may be nec-essary to detect a small acoustic perturbation to ambientconditions, such as, for example, emanation from workingequipment. It may be important to detect the spatial loca-tion of the perturbation with high accuracy. To this end,one might look for a sensor with the effect of amplifiedresolution.It may be advantageous to construct an arrangement ofchains as the one shown in Fig. 1, in the form of a paral-lel array (of such chains). The idea is to excite compactantisymmetric periodic modes localized at different posi-tions along the chain. For example, positioning the CBsite at locations shifted a given number of representativeelements in the direction of the chains, such that this shiftis increased by a fixed amount from one chain to the nextone. This design may have advantages. If the magnitude of a single shift along the chain, when passing from one chainto the next one, is an order of magnitude larger than thespatial distance between two adjacent chains, then thereis already an order of magnitude gain in the amplificationfactor of the sensor.The concept of the device is that when a small (local-ized) acoustic perturbation reaches a specific chain in thearray, it can destroy the CB-mode excited (in advance) inthat chain. The underlying assumption is that the param-eter (the mass ratio, for example) was chosen such that theCB is marginally-stable (so that it would require a finite,even if small, perturbation to cause instability).By observing the location of the destroyed CB, it wouldbe possible to know the position of the external perturba-tion in question (see Fig. 17). The high accuracy of thesensing would be based on the large distance between pre-existing CBs in the chain direction, as compared to thedistance between the chains (which is the actual spatialresolution).As for the mechanism of ‘reading’ the device, clearly itshould not be problematic to determine the CB of whichposition (along the chain) was destroyed. Due to the rela-tively large intervals between the CB sites, a low-resolutionsensor can be used to ‘read’ the proposed device, and inconjunction with the calibration (associating the CB posi-tions with the chain positions in the array – which could bean affine relation), an overall high-resolution sensor couldbe obtained.Possible challenges in the engineering of the sensor couldbe the assembling of the chains in a casing, such that in-teraction between the chains is avoided. A second crucialissue is damping. The solution to the second challenge maybe the addition of small background periodic excitation.Advanced-stage analysis of the device may requirestudying an array of the chains as suggested, with theaddition of weak strongly-nonlinear coupling between thechains and introduction of small forcing and damping. Astability chart for the overall device with the mentionedadditions would have to be constructed.Two points are importing to note. The first is that thetranslational invariance of the lattice is crucial for detect-ing moving acoustic perturbations, especially if a moving(tracking) device is sought. The second point is that in or-der to avoid false positives, a threshold excitation shouldbe designed for, meaning that the system cannot simplybe tuned to be unstable. Rather marginal stability shouldbe sought, assuming that nonlinear instability is likely tobe associated with marginal linear stability. To this end,it is crucial for the stability bounds to be tractable. Hencethe large emphasis put on the matter in the present study,which proved to be reasonably successful, at least for high-symmetry instabilities.22 δ x Fig. 17:
Schematic sketch of the suggestion for the sensor –the spikes represent antisymmetric ( z -mode) discrete CBs lo-cated on shifted localization-sites in an array of translationally-invariant chains, with a schematic spatially-resolved external-mode incident acoustic perturbation
7. Conclusions
The present paper proposes a system consistent of a one-dimensional nonlinear non-integrable classical lattice de-scribing the exact dynamics of the fundamental degrees offreedom (rather than averaged quantities) in the mechan-ical setting. The system is comprised of ‘on-site‘ complexblocks containing several degrees of freedom each, charac-terized by the existence of symmetry between the degreesof freedom within each block.The aforementioned symmetry is one that allows spatialdecoupling of a single block from the rest of the lattice,corresponding to antisymmetric (momentumless) mode ofmotion. In linear analysis (for the case of nonlinearityhaving a weakly-nonlinear limit), the system shows threedispersion curves, an acoustic branch, an optical branch,and a flat band, which has a tangent point with one of thetwo other bands.Nonlinear analysis was performed by assuming a β -FPUtype interaction potential. The analysis of a single repre-sentative cell, assuming conservation of momentum andenergy, can be performed by rigorous implementation ofHill’s method. Standard instability tongues emerge. Aqualitatively similar picture of instability tongues is ob-tained for a chain of three blocks, by numerical integra-tion in the framework of the Floquet theory, with exactdetection of the two boundaries for each tongue. A verylarge number of overlapping tongues is observed. Twoof the tongues, related to the pitchfork bifurcation, arereproduced also for the case of a chain comprised of 19blocks, where Floquet theory is, again, implemented usingnumerical integration. This time, however, the instability-tongues’ boundaries are not detected exactly. Rather,the parameter plane (mass-ratio–amplitude) is divided bya dense grid and each node is determined to be either stable or unstable. A very complicated instability pat-tern emerges. Apart from the two principal pitchfork-bifurcation-related tongues evident already for the three-blocks chain, an additional, comb-like, array of pitchfork-instability-related tongues is observed. This comb-like ar-ray of instability tongues falls into a region the bound-aries of which can be theoretically derived. In fact, thoseboundaries appear to be related to resonance between thecompact antisymmetric nonlinear normal mode (associ-ated with a flat dispersion band) and the linear propaga-tion spectrum. This renders the stability map tractable forlow amplitudes, at least for the case of the high-symmetrypitchfork-bifurcation-related instability.For the case of large amplitudes, almost the entire pa-rameter plane appears to become unstable, with two ex-ceptions. One is a narrow but finite stability stripe, relatedto a gap between the two principal pitchfork-bifurcation-related instability tongues emerging already for the chainwith three elements. The second exception corresponds tostability for an isolated value of the mass ratio.A major part of the paper is dedicated to the effortto explain the aforementioned emerging isolated stability-preserving mass-ratio value. To this end, the smoothlynonlinear large-amplitude regime is associated with a non-smooth system with conservative impact interaction (andlinear dynamics between impacts). We use a constructionof logical arguments, energy considerations, and direct in-tegration of a dynamic mapping equation. In addition,we employ the saltation matrix realization of the Floquettheory for piecewise linear systems, applying asymptoticconsiderations of various sorts. The overall effort turnsout to be relatively successful and the isolated-parameter-value-related stability is explained.Finally, based on the property of translational invarianceof the system, in which internal symmetry and the associ-ated flat-band produce a perfectly-localized periodic mode,and relying on the more-or less tractable stability pictureof the system, possible practical implementation is sug-gested. This implementation consists of an acoustic mo-tion sensor (where translational invariance is important),for which high spatial resolution of detecting acoustic per-turbations can be obtained as follows. Spatially-close localperturbations are related to well-distant ‘standing waves’in an array of chains. In such a device, the perfect com-pactness of the standing waves (before they are destructedby the perturbations, while those are being detected), canallow ‘cleaner’ detection (since even after destruction, thecompactness holds for some time). Furthermore, the anal-ysis capabilities, as presented in this work, allow betterdesign and fine-tuning of the aforementioned (still hypo-thetical) device (for example, by tuning the parameters tomarginal stability for then relying on possible associatednonlinear instability, to obtain an excitation threshold).The present work lies within the emerging body of liter-ature of the recent years related to such concepts as appli-cations of flat bands, dynamic localization, meta-materialsand resonance-based sensors.23he main contribution of this work is the imple-mentation of a perfectly-compact-breather solution in atranslationally-invariant nonlinear mechanical system, andthe use of the vibro-impact system analogy for study ofotherwise hardly-tractable stability characteristics. Acknowledgments
The financial help of the Israel Science Foundation,Grant No. 1696/17, is gratefully acknowledged.
References [1] S. Aubry, Breathers in nonlinear lattices: Existence, linear sta-bility and quantization, Physica D103 (1997) 201-250.[2] P. Rosenau and A. Zilburg, On a strictly compact discretebreather in a Klein-Gordon model, Physics Letters A 379 (43-44) (2015) 2811–2816.[3] D. Leykam, S. Flach, O. Bahat-Treidel, and A. Desyatnikov, Flatband states: disorder and nonlinearity, Physical Review B 88(2013) 224203.[4] W. Maimaiti, A. Andreanov, H. C. Park, O. Gendelman, andS. Flach, Compact localized states and flatband generators inone dimension, Physical Review B 95 (2017) 115135.[5] M. Johansson, U. Naether and R. A. Vicencio, Compactificationtuning for nonlinear localized modes in sawtooth lattices, Phys-ical Review E 92 (2015) 032912.[6] K. Zegadlo, N. Dror, N. V. Hung, M. Trippenbach, and B. A.Malomed, Single and double linear and nonlinear flatband chains:spectra and modes, Physical Review E 96 (2017) 012204.[7] C. Danieli, A. Maluckov, and S. Flach, Compact descretebreathers on flat-band networks, Low Temperature Physics 44(2018) 678-687[8] N. Perchikov and O. V. Gendelman, Flat bands and compactonsin mechanical lattices, Physical Review E 96 (2017) 052208.[9] A. Sergyeyev, S. Skurativskyi, and V. Vladimirov, Compactonsolutions and(non)integrability of nonlinear evolutionary PDEsassociated with a chain of prestressed granules, Nonlinear Anal-ysis: Real World Applications 47 (2109) 68–84.[10] J. Cuevas-Maraver, P. G. Kevrekidis, B. A. Malomed, andL. Gao, Solitary waves in the Ablowitz–Ladik equation withpower-law nonlinearity, Journal of Physics A: Mathamatical andTheoretical 52 (2019) 065202.[11] G. C´aceres-Aravena and R. A. Vicencio, Perfect localiza-tion on flat band binary one-dimansional photonic lattices,arXiv:1903.00377.[12] G. James, Travelling breathers and solitary waves in stronglynonlinear lattices, Philosophical Transactions of the Royal Soci-ety A 376 (2018) 1–25.[13] B. Real and R. A. Vicencio, Controlled mobility of compact dis-crete solitons in nonlinear Lieb photonic lattices, Physical ReviewA 98 (2018) 053845.[14] J. Qin, Z. Liang, B. A. Malomed, and G. Dong, Tail-free self-accelerating solitons and vortices, Physical Review A 99 (2019)023610.[15] S. W. Kim and S. Kim, Fano resonances in translationally-invariant nonlinear chains, Physical Review B 63 (2000) 212301.[16] W. Maimaiti, S. Flach, and A. Andreanov, Universal d = 1 flatband generator from compact localized states, Physical ReviewB 99 (2019) 125129.[17] J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, N. I.Karachalios, M. Haragus, and G. James, Floquet analysis ofKuznetsov-MA breathers: A path toward spectral stability ofrogue waves, Physical Review E 96 (2017) 012202.[18] J. Cuevas-Maraver, P. G. Kevrekidis, A. Vainchtein, andH. Xu, Unifying perspective: Solitary traveling waves as discretebreathers in Hamiltonian lattices and energy criteria for theirstability, Physical Review E 96 (2017) 032214. [19] Y. Doi and K. Yoshimura, Symmetric Potential Lattice andSmooth Propagation of Tail-Free Discrete Breathers, PhysicalReview Letters 117 (2016) 014101.[20] S. V. Suchkov, B. A. Malomed, S. V. Dmitriev, and Y. S.Kivshar, Solitons in a chain of parity-time-invariant dimers,Physical Review E 84 (2011) 046609.[21] N. Perchikov, O. V. Gendelman, Nonlinear dynamics of hiddenmodes in a system with internal symmetry, JSV 377 (2016) 185–215.[22] N. Perchikov, O. V. Gendelman, Dynamics and stability of a dis-crete breather in a harmonically excited chain with vibro-impacton-site potential, Physica D 292 (2015) 8–28. Appendix A. Details of the Floquet-theory analy-sis of the CB for the smooth system
The matrix B ( τ ) appearing in Eq. (32) can be presentedin a block-form as: B ( τ ) = − B ( xx ) B ( xy ) N − ,N B ( yx ) − B ( yy ) N N,N − N − B ( zz ) (A.1)where p,q is a rectangular zeros-matrix of p rows and q columns. The components of the four defining rectangu-lar matrices (with occasional commas between subscriptsadded for clarity) are given by B ( xx ) ij = β ( xx ) i δ ij , β ( xx ) i , − ∂ ˆ x ′′ i +1 ∂ ˆ x i +1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 =1 µ (cid:20) α ˆ z ( τ )( δ N ,i +1 + δ N ,i ) (cid:21) , ≤ i ≤ N − ,B ( xy ) ij = ∂ ˆ x ′′ i +1 ∂ ˆ y j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 14 µ [ δ i +1 ,j + δ ij ++ 32 α ˆ z ( τ )( δ N ,i +1 δ i +1 ,j + δ N ,i δ ij ) (cid:21) , ≤ j ≤ N ; B ( yx ) ij = ∂ ˆ y ′′ i ∂ ˆ x j +1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = (cid:20) α ˆ z ( τ ) δ N ,i (cid:21) ×× ( δ i,j +1 + δ i,j ) , ≤ i ≤ N − , ≤ j ≤ N − B ( yx ) N,j = ∂ ˆ y ′′ N ∂ ˆ x j +1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 == (cid:20) α ˆ z ( τ ) δ N ,N (cid:21) δ N,j +1 ; B ( yy ) ij = − ∂ ˆ y ′′ i ∂ ˆ y j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = (cid:20) α ˆ z ( τ ) δ N i (cid:21) δ ij , ≤ i ≤ N − , ≤ j ≤ N ; B ( yy ) N,j = − ∂ ˆ y ′′ N ∂ ˆ y j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 == (cid:20)
12 + 34 α ˆ z ( τ ) δ N ,N (cid:21) δ N,j , ≤ j ≤ N (A.2)24 ( zz ) ij = − ∂ ˆ z ′′ i ∂ ˆ z j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = (cid:20) α ˆ z ( τ ) δ N i (cid:21) δ ij , ≤ i ≤ N − , ≤ j ≤ N ; B ( zz ) N,j = − ∂ ˆ z ′′ N ∂ ˆ z j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = (cid:20)
12 + 34 α ˆ z ( τ ) δ N ,N (cid:21) δ N,j , ≤ j ≤ N ; (A.3)The remaining expressions for the components of B ( yx )1 ,j and B ( yy )1 ,j are given by B ( yx )12 = ∂ ˆ y ′′ ∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 0 , B ( yx )1 ,j = ∂ ˆ y ′′ ∂ ˆ x j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = − − α ˆ z ( τ ) δ N , ∀ j > , B ( yy )11 = − ∂ ˆ y ′′ ∂ ˆ y (cid:12)(cid:12)(cid:12)(cid:12) ˆ x , ˆ y =0 = 1 + 2 µ µ (cid:20) α ˆ z ( τ ) δ N , (cid:21) , B ( yy )1 ,j = − ∂ ˆ y ′′ ∂ ˆ y j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 12 µ (cid:20) α ˆ z ( τ ) δ N , (cid:21) ∀ j ≥ δ ab or δ a,b is Kronecker’s delta, and ˆ z ( τ ) is the (pe-riodic) CB solution at the CB localization-site, which isidentical to the solution given in Eqs. (17) and (19). A.1. Periodic boundary conditions
For periodic boundary conditions (and the assumptionof zero total momentum and displacement, as before), Eq.(29) should be replaced by the following expressions:ˆ y ′′ N = − µ µ ˆ y N − N − X n =2 ˆ x n − µ N − X n =1 ˆ y n ++ α x N − ˆ y N )[(2ˆ x N − ˆ y N ) + 3ˆ z N ] − α N X n =2 ˆ x n + 1 µ N − X n =1 ˆ y n + 1 + µµ ˆ y N ! ×× N X n =2 ˆ x n + 1 µ N − X n =1 ˆ y n + 1 + µµ ˆ y N ! + 3ˆ z N ˆ z ′′ N = − ˆ z N − α (cid:2) (2ˆ x N − ˆ y N ) + N X n =2 ˆ x n + 1 µ N − X n =1 ˆ y n + 1 + µµ ˆ y N ! ˆ z N + 2ˆ z N (A.5)Then, the fourth and sixth rows in Eqs. (A.2) and thesecond row in Eq. (A.3) should be replaced by the fol-lowing expressions (much simplified, since in a ring, theCB site can be chosen arbitrarily and there is no point in setting N = N ): B ( yx ) N,j = ∂ ˆ y ′′ N ∂ ˆ x j +1 (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = − , ≤ j ≤ N − B ( yx ) N,N − = ∂ ˆ y ′′ N ∂ ˆ x N (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 0 B ( yy ) N,j = − ∂ ˆ y ′′ N ∂ ˆ y j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 12 µ , ≤ j ≤ N − B ( yy ) N,N = − ∂ ˆ y ′′ N ∂ ˆ y N (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = 1 + 2 µ µB ( zz ) N,j = − ∂ ˆ z ′′ N ∂ ˆ z j (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ y =0 = δ N,j , ≤ j ≤ N (A.6) Appendix B. Details of the Floquet-theory analy-sis of the CB – nonsmooth system
The matrices required for analytical construction of themonodromy are as follows: C j,l , " I ˆ N ˆ N ˆ N I ˆ N + ∆ ˆ C ( j,l ) , A ∗ , (cid:20) ˆ N I ˆ N ˆ A 0 ˆ N (cid:21) , ˆ A , − µ I N − µ G N − ,N µ G N − ,N ( G N − ,N ) ⊤ − I N N ( G N − ,N ) ⊤ N − I N (B.1)where N is a matrix of zeros of size N , I N is the identitymatrix of rank N, ˆ N , N − G ij , δ ij + δ i,j − .∆ ˆ C (1 , ik = ∆ ˆ C (3 , ik , − δ ik ( δ i,n − + δ i, n − )+ δ i,n − δ k, n − + δ n − ,k δ i, n − +∆ ˆ C (2 , ik = ∆ ˆ C (3 , ik , − δ ik ( δ i,n − + δ i, n − )+ δ i,n − δ k, n − + δ n − ,k δ i, n − +∆ ˆ C (1 , ik , − δ ik ( δ i,n + δ i, n − )++ δ i,n δ k, n − + δ n ,k δ i, n − ∆ ˆ C (2 , ik , − δ ik ( δ i,n + δ i, n − )++ δ i,n δ k, n − + δ n ,k δ i, n − (B.2) S j,l = C j,l + ( A ∗ C j,l − C j,l A ∗ ) D j,l (B.3) D , , Lx [ n (1 , ] ⊤ [ n (1 , ] ⊤ A ∗ Lx , D , , C , Lx [ n (2 , ] ⊤ [ n (2 , ] ⊤ A ∗ C , Lx , D , , C , C , Lx [ n (1 , ] ⊤ [ n (1 , ] ⊤ A ∗ C , C , Lx , D , , C , Lx [ n (2 , ] ⊤ [ n (2 , ] ⊤ A ∗ C , Lx , D , , C , C , Lx [ n (2 , ] ⊤ [ n (2 , ] ⊤ A ∗ C , C , Lx , D , , C , C , C , Lx [ n (1 , ] ⊤ [ n (1 , ] ⊤ A ∗ C , C , C , Lx (B.4)25here n (1 , i = δ i, n − − δ i, n − , n (2 , i = δ i, n − − δ i, n , n (1 , i = δ i, n − − δ i, n +1 , and n (2 , i = δ i, n − δ i, n +1 , for i = 1 , , .., N and δ a,b is Kronecker’s delta. M ( p ) = [ LS ( p ) L ] , L = exp (cid:16) π ω A ∗ (cid:17) , S (3) = S , S , S , , S (4) = S , S , S , S , (B.5)where ˆ ω = ω/ p k/m is the normalized frequency of the pe-riodic solution ( ω being the frequency), which should beabove the flat-band frequency (ˆ ω FB = 1). The initial con-dition is x = ( q ˆ N , p ˆ N ) ⊤ , q i = 0 , p i = δ i, n − − δ i, n − .The (absolute) value of the initial velocity of the elementsof the CB appears explicitly only in Eq. (B.4), whereit contracts, and thus one can understand x as a nor-malized quantity. The initial velocity still influences thevalue of the period of the CB, represented by ˆ ω . Equation(B.1) assumes that one (say, the leftmost) mass M is eitherfixed or its displacement and velocity are known from totalmomentum conservation (‘free-free’ boundary conditions).Therefore, the number of independent displacements is as-sumed to be 3 N −
1. Also, the second boundary needs acondition, say that of a free boundary. Unlike for the caseof moderate amplitudes, where boundary conditions affectthe emergence of Hopf instabilities, it appears that forlarge amplitudes the Hopf instabilities nearly vanish andonly pitchfork instabilities remain, for which the specificchoice of boundary conditions is less important, for a longchain.Moreover, the choice of boundary conditions does not af-fect the emergence of a stable CB for large amplitudes forthe ratio
M/m = 1 (the value for which there is reason toconstruct the monodromy matrix in the nonsmooth case).Therefore, Eq. (B.1) is corrected to represent the simplestset of boundary conditions – ‘fixed-free’. It is assumedthat x is known and needs not an equation. Thus, thequantity G N − ,N is understood to be having its first rowcorresponding to the second of Eqs. (1). The same goes forthe equations defining C j,l , D j,l and S j,l . Accordingly, thematrix G N,N − has its first column corresponding to thevariable x . This covers the ‘left’ (fixed) boundary condi-tion. Furthermore, if one assumes that the ‘right’ bound-ary condition is free, then the matrix I N in the diagonalblocks of ˆ A is defined as I ij = δ ij ∀ i < N, I NN = 1 / µ = 1 / A symmetric for an infinite chain or periodic boundary con-ditions). An odd value of N should be taken, preferablywith n = ( N + 1) /
2, for maximum distance of the CBfrom the boundaries.
Appendix C. Analysis of the three-body dynamicsin the smooth limit for ˆ µ ≪ Rearranging Eq. (63), we get the following two-impactsfinite-difference equation: v ( k +2)1 − v ( k )1 = 2ˆ µ h v ( k +1)3 − v ( k )1 i (C.1) The expression for v ( k +1)3 can be obtained from Eq. (42),by taking the ˆ µ ≪ v ( k +1)3 = 2 v ( k )2 − v ( k )3 + 2ˆ µ h v ( k )3 − v ( k )2 i (C.2)Substituting Eq. (C.2) into Eq. (C.1), one gets, tofirst order, the following explicit two-impacts differenceequation: v ( k +2)1 − v ( k )1 = 2ˆ µ h v ( k )2 − v ( k )3 − v ( k )1 i (C.3)Next, in order to obtain v ( k +2)3 , one resorts again to Eq.(42) in the ˆ µ ≪ v ( k +2)3 = 2 v ( k )1 − v ( k +1)3 + 2ˆ µ h v ( k +1)3 − v ( k )1 i (C.4)Substituting Eq. (C.2) into Eq. (C.4) results, to firstorder in ˆ µ , in the following explicit two-impacts differenceequation:: v ( k +2)3 − v ( k )3 = 2 v ( k )1 − v ( k )2 ++2ˆ µ h v ( k )2 − v ( k )3 − v ( k )1 i (C.5)Recalling the equation of conservation of linear momen-tum for the three-body system, namely, v ( k )2 = − v ( k )1 − ˆ µv ( k )3 (C.6)allows eliminating v ( k )2 , and writing Eqs. (C.3) and (C.5)as a closed system: v ( k +2)1 − v ( k )1 = − µ h v ( k )1 + v ( k )3 i v ( k +2)3 − v ( k )3 = 4 v ( k )1 − µ h v ( k )1 + v ( k )3 i (C.7)(omitting terms nonlinear in ˆ µ ).In Eqs. (C.7), the updated variables refer to the valuesof the velocities after impacts involving particles 1 and 3.Hence, the change due to two impacts in the velocity ofparticle 3 is negative as long as the velocity of particle 1is negative. For the positive points in the history of thevelocity of particle 3, one would have to use the solutionof the system in Eqs. (C.7), along with Eqs. (C.2) and(C.6) – this is what the calculation in Eq. (66) requires.The sequence v ( k )3 is sign-changing (between consecutiveentries in the sequence). However, if only odd values of k are observed, then a negative sequence emerges, with nosign-changes within it.Now, the changes between elements in this sequence arefinite, as well as the differences between the instances.Therefore, the left-hand sides in Eqs. (C.7) can be un-derstood as numerical derivatives with respect to the in-teger variable ¯ k , where k = 2¯ k −
1, such that v ′ , ∼ [ v (¯ k +1)1 , − v (¯ k )1 , ] / [(¯ k + 1) − ¯ k ], where the prime implies dif-ferentiation with respect to increments (of magnitude 1)of ¯ k . Then one gets the following relations: v ′ ∼ − µ (3 v + v ) , v ′ ∼ v − µ (4 v + v ) (C.8)26the superscripts are omitted henceforth for clarity).Defining ˜ v (¯ k )1 , ∼ v ( k )1 , , one can rewrite the above relationsas two coupled first-order differential equations, which canbe decoupled analytically into a single second-order linearequation and an auxiliary relation, as follows:˜ v ′′ + 8ˆ µ ˜ v ′ + 8ˆ µ ˜ v = 0 , ˜ v = (cid:18)
14 + ˆ µ (cid:19) ˜ v ′ + ˆ µ v (C.9)(omitting terms nonlinear in ˆ µ ).The first (second-order) equation in Eqs. (C.9) can bereadily solved, yielding the solution˜ v (¯ k )3 = V e − µ ¯ k sin (cid:16)p µ ¯ k + φ (cid:17) (C.10)(where for the sake of the well-posedness of the derivative,¯ k can be hypothetically extended over the field of reals).Consequently, the velocity associated with particle 1would to leading order in ˆ µ be˜ v (¯ k )1 = √ µ V e − µ ¯ k cos (cid:16)p µ ¯ k + φ (cid:17) (C.11)Therefore, using the definitions introduced above andemploying Eqs. (C.6) and (C.2), one would have, to lead-ing order in ˆ µ , v ( k even)1 ∼ √ µ V e − µ ¯ k cos hp µ ¯ k + φ i == v ( k odd)1 , v ( k odd)2 → ˆ µ ≪ − v ( k odd)1 ⇒ v ( k even)3 → ˆ µ ≪ − v ( k odd)1 − v ( k odd)3 → ˆ µ ≪ − v ( k odd)3 ∼∼ V e − µ ¯ k sin hp µ ¯ k + φ − π i (C.12)which implies that, to leading order in ˆ µ , each next impactmerely changes the sign of the velocity of particle 3 aftera given impact (rather than also doubling its magnitude,as in the strict ˆ µ = 0 case – a sign of singularity of theintegrability, in the sense of velocity-reversal, with ˆ µ as itapproaches zero).Now, we see that the sequence of velocities of particle3 after an even number of impacts is governed by two nu-merical time scales, one associated with exponential decayon a typical time of τ e ∼ (4ˆ µ ) − , and the other associatedwith sinusoidal change with a period of τ s ∼ π (8ˆ µ ) − / ,which for ˆ µ ≪ τ e ).Since we are interested in integrable (periodic) solutions,the relevant time scale of consideration would be (half)a period (namely, τ s / e − π √ µ , which forˆ µ ≪ − π √ µ , and leads to a two-termsexpression for v ( k even)3 . In this two-terms expression, thesecond term is negligible with respect to the first term inthe limit of ˆ µ ≪
1. Therefore, to leading order in ˆ µ , onehas: v ( k even)3 ∼ V sin hp µ ¯ k + φ − π i (C.13) Next, enforcing the initial condition, v (0)3 = 0, leadsto the determination of the phase as φ = π − √ µ andthis, along with determining the amplitude from Eq. (65),yields the following expression (correct to first-order) forthe sequence of velocities of particle 3 after an even numberof impacts: v ( k )3 ∼ v p ˆ µ/ hp µk i , ∀ k even (C.14)(where one recalls that k = 2¯ k −
1, as defined earlier).Clearly, the obtained velocity evolution is sinusoidal , as‘assumed’ in the derivation of Eq. (67). Indeed, the sinefunction is concave on the half-period where it is positive;it is symmetric around the middle of this range; and it islinear in the vicinity of the edges of this range, starting andending as zero – all the requirements from the distributionas dictated by periodic dynamics in the ˆ µ ≪ ρ from Eq. (66)), amounts to 2 /π , which isindeed the value used for the derivation of Eq. (67).It is noteworthy that another way to validate Eq. (67) isto require that the argument of the sine in Eq. (C.14) be-comes exactly π just when k = j and ˆ µ takes the value sug-gested by Eq. (67) (and one bears in mind that v ( k = j )3 = 0is a condition required for velocity reversal).This is not trivial, since the functional form of ˆ µ j inEq. (67) was obtained from energy conservation and kine-matic considerations for the discrete system representedby the mapping. However, here, the same functional formof inverse square dependence is obtained due to the struc-ture of the coefficients of a second-order linear differentialequation derived by taking ‘the smooth’ limit (convertingthe mapping into a flow). The fact that identical func-tional forms were obtained, as aforementioned, is valida-tion of the self-consistency of the process of resorting tothe smooth-system limit (for the ˆ µ ≪ ρ , which is required forobtaining Eq. (67). Therefore, the analysis presented inthis appendix, although it may appear as (partially) basedon cyclic logic, does in fact provide independent validationnecessary for checking overall self-consistency.As a last note, one observes that taking the strict limitˆ µ → µ ≪
1) in Eq. (C.14), leads tothe result v ( k even)3 = 2 kv for small values of k . Acomplementary-symmetric result is also obtained for k → π/ √ µµ