#### Abstract

We consider time-periodic structures of granular crystals consisting of alternate chrome steel (S) and tungsten carbide (W) spherical particles where each unit cell follows the pattern of a 2 : 1 trimer: S-W-S. The configuration at the left boundary is driven by a harmonic in-time actuation with given amplitude and frequency while the right one is a fixed wall. Similar to the case of a dimer chain, the combination of dissipation, driving of the boundary, and intrinsic nonlinearity leads to complex dynamics. For fixed driving frequencies in each of the spectral gaps, we find that the nonlinear surface modes and the states dictated by the linear drive collide in a saddle-node bifurcation as the driving amplitude is increased, beyond which the dynamics of the system becomes chaotic. While the bifurcation structure is similar for solutions within the first and second gap, those in the first gap appear to be less robust. We also conduct a continuation in driving frequency, where it is apparent that the nonlinearity of the system results in a complex bifurcation diagram, involving an intricate set of loops of branches, especially within the spectral gap. The theoretical findings are qualitatively corroborated by the experimental full-field visualization of the time-periodic structures.

#### 1. Introduction

Granular chains, which consist of closely packed arrays of particles that interact elastically, have proven over the last several decades to be an ideal testbed to theoretically and experimentally study novel principles of nonlinear dynamics [1–3]. Examples include, but are not limited to, solitary waves [1, 2, 4, 5] and dispersive shocks [6–8], as well as bright and dark discrete breathers [9–15]. Beyond such fundamental aspects, their extreme tunability makes granular crystals relevant for numerous applications such as shock and energy absorbing layers [16–19], actuating devices [20], acoustic lenses [21], acoustic diodes [12], and switches [22], as well as sound scramblers [23, 24].

Our emphasis in the present work will be on coherent nonlinear waveforms that are time-periodic. A special instance of this is when the spatial profile is localized, in which case the structure is termed a discrete breather. The study of discrete breathers has been a topic of intense theoretical and experimental interest during the 25 years since their theoretical inception, as has been summarized, for example, in [25]. The broad and diverse span of fields where such structures have been of interest includes, among others, optical waveguide arrays or photorefractive crystals [26], micromechanical cantilever arrays [27], Josephson-junction ladders [28, 29], layered antiferromagnetic crystals [30, 31], halide-bridged transition metal complexes [32], dynamical models of the DNA double strand [33], and Bose-Einstein condensates in optical lattices [34].

In Fermi-Pasta-Ulam type settings (which are intimately connected to the case of precompressed granular crystals in the weakly nonlinear regime, see, e.g., (1)) it was proven in [35, 36] (see also the discussion in [25]) that small amplitude discrete breathers are absent in spatially homogeneous (i.e., monoatomic) chains. Instead, dark breathers (those on top of a non-vanishing background) have been found therein [15]. It is for that reason that the first theoretical and experimental investigations of breathers with a vanishing background (i.e., bright breathers) have taken place in granular chains with some degree of spatial heterogeneity, which plays a critical role in the emergence of such patterns. Examples include chains with defects [9, 10, 12] (see also [37] for recent experiments) and a spatial periodicity of two (i.e., dimer lattices) [11, 13]. Motivated by these works, a theoretical study of breathers in granular crystals with higher order spatial periodicity (such as trimers and quadrimers) was recently conducted in [14]. Therein, it was demonstrated that breathers with a frequency in the highest finite gap appear to be more robust than their counterparts with frequency in the lower gaps. Another relevant study is [38], where experimental and numerical analysis of band gap effects related to transient behavior of driven dissipative granular chains were studied.

The goal of the present work is the systematic study of time-periodic solutions (including breathers) of trimer granular crystals with frequency in the first or second gap, as well as in the acoustic and optical bands. In particular, we investigate the robustness of the breathers experimentally using a full-field visualization technique based on laser Doppler vibrometry. This is a significant improvement over the aforementioned experimental observation of bright breathers [11–13] where force sensors are placed at isolated particles within the granular chain, which does not allow a full-field representation of the breather. We complement this study with a detailed theoretical probing of the more realistic damped-driven variant of the pertinent model. Our extensive analysis of such modes consists of the study of their family under continuations in both the amplitude and the frequency parameters of the external drive and a detailed comparison of the findings between numerically exact (up to a prescribed tolerance) periodic orbit solutions and experimentally traced counterparts thereof. We also note that there are limited studies reporting the full-field experimental visualization of time-periodic solutions in granular crystals. Indeed, the only work to this effect that the authors are aware of is [15], which deals with the much simpler setting of a monomer granular crystal, where surface breathers are not addressed. Although the theme of heterogeneous granular crystals has been receiving considerable theoretical and even experimental attention recently, most of it has been focused on the realm of traveling waves [39–42], while the emphasis herein will be on the time-periodic, exponentially localized discrete breather states.

The paper is structured as follows. In Sections 2 and 3, we describe the experimental and theoretical setups, respectively. The main results are presented in Section 4, where time-periodic solutions with frequency in the first/second gaps and in the spectral bands are studied in both of these setups and compared accordingly. Finally, Section 5 provides our conclusions and discusses a number of future challenges.

#### 2. Experimental Setup

Figure 1 shows a schematic of the experimental setup consisting of a granular chain and a laser Doppler vibrometer. In this study, we consider a granular chain composed of spherical beads. The beads are made out of chrome steel (S: gray particles in Figure 1) and tungsten carbide (W: black particles) materials. See Table 1 for nominal values of the material parameters used hereafter. The granular chain has a spatial periodicity of three particles, and each unit cell follows the pattern of a 2 : 1 trimer: S-W-S. The spheres are supported by four polytetrafluoroethylene (PTFE) rods, which allow axial vibrations of particles with minimal friction, while restricting their lateral movements. The granular chain is compressed by a moving wall at one end of the chain that applies static force in a controllable manner via a linear stage (see Figure 1). We measure the preapplied static force ( N in this study) by using a static force sensor mounted on the moving wall. We assume that this moving wall is stationary throughout our analysis, since it exhibits orders-of-magnitude larger inertia compared to the particle’s masses.

The granular chain is driven by a piezoelectric actuator positioned on the other side of the chain. We impose actuation signals of chosen amplitude and frequency through an external function generator and a power amplifier. The dynamics of individual particles are scoped by a laser Doppler vibrometer (LDV, Polytec, OFV-534), which is capable of measuring particles’ velocities with a resolution of 0.02 *μ*m/s/Hz^{1/2}. The LDV scans the granular chain through the automatic sliding rail and measures the vibrational motions of each particle three times for statistical purposes. We obtain the full-field map of the granular chain’s dynamics by synchronizing and reconstructing the acquired data.

#### 3. Theoretical Setup

The equation we use to model the experimental setup is a Fermi-Pasta-Ulam-type lattice with a Hertzian potential [1] leading towhere , is the displacement of the th bead from its equilibrium position at time , is the mass of the th bead, is a precompression factor induced by the static force , and the bracket is defined by . The 3/2 power accounting for the nonlinearity of the model is a result of the sphere-to-sphere contact, that is, the so-called Hertzian contact [43]. We model the dissipation as a dash-pot, which can be interpreted as the friction between individual grains and PTFE rods. This form of dissipation has been utilized in the context of granular crystals in several previous works [12, 15], but it is worth noting that works such as [6, 7] considered the internal friction caused by contact interaction between grains. The strength of the dissipation is captured by the parameter , which serves as the sole parameter used to fit experimental data (ms in this study). The elastic coefficient depends on the interaction of bead with bead and for spherical point contacts has the form [44]where , , and are the Young’s modulus, Poisson’s ratio, and the radius, respectively, of the th bead. The left boundary is an actuator and the right one is kept fixed; that is,where and represent the driving amplitude and frequency, respectively. Note that at the boundaries (i.e., and ), we consider flat surface-sphere contacts while the same material properties are assumed; thus, for example, the elastic coefficient at the left boundary takes the form

##### 3.1. Linear Regime and Dispersion Relation

Assuming the dynamic strains that are small relative to the static precompression, that is,equation (1) can then be well approximated by its linearized formwith corresponding to the linearized stiffnesses. Subsequently, (6) can be converted into a system of first-order equations and written conveniently in matrix form aswith andwhere and represent the zero and identity matrices, respectively, and the (tridiagonal) matrix is given byNote that in the above linearization we applied fixed boundary conditions at both ends of the chain; that is, (we account for the actuator in the following subsection). This equation is solved by , where corresponds to the angular frequency (with ) and wherewith corresponding to the eigenvalue-eigenvector pair, while and . The eigenfrequency spectrum using the values of Table 1 and ms is shown in Figure 2(a).

**(a)**

**(b)**

On the other hand, the dispersion relation for the trimer chain configuration within an infinite lattice can be obtained using a Bloch wave ansatz [45]. We first rewrite (6) in the following convenient formwith , where we made use of the spatial periodicity of the trimer lattice; that is, and and and . The Bloch wave ansatz has the formwhere , , and are the wave amplitudes, is the wavenumber, and is the size (or the equilibrium length) of one unit cell of the lattice. Substitution of (12a), (12b), and (12c) into (11a), (11b), and (11c) yieldsThe nonzero solution condition of the matrix-system (13), or dispersion relation (i.e., the vanishing of the determinant of the above homogeneous linear system), has the formwhich has (nontrivial) solutions at and (i.e., first Brillouin zone) given bywhere . The solutions (15) correspond to the* cut-off frequencies* of the spectral bands. The above values are (in principle) complex since we consider dissipative dynamics in the system (1) (embodied by the term) in order to model the experimental setup. Therefore, the reported frequencies hereafter will correspond to the* real* part of these values, that is, the part contributing to the dispersion properties of the solutions rather than their decay. In Table 2, the predicted values of the cut-off frequencies are presented (up to two decimal places) using the values of the material parameters of Table 1 and ms.

Finally, upon solving numerically the dispersion relation (cf. (14)), the frequency as a function of the dimensionless wavenumber () is presented in Figure 2(b) together with the cut-off frequencies of Table 2. Note that these cut-off frequencies are also plotted as horizontal dashed black lines in Figure 2(a) for comparison. It can be discerned from both panels of Figure 2 that there exist two finite gaps in the frequency spectrum (in general their number is one less than the period of the granular crystal), namely, and , together with one semi-infinite gap . We also note the presence of two additional localized modes depicted in Figure 2(a) between the first and second optical pass bands due to the fixed boundary conditions. These modes, denoted by black dots in Figure 2(a), correspond to surface modes and their existence has been previously discussed, for example, in [11].

##### 3.2. Steady-State Analysis in the Linear Regime

In order to account for the actuation in the linear problem we add an external forcing term to (7) in the formwhere the sole nonzero entry of is at the th node and has the form and the matrix is given by (8). An inhomogeneous solution of (16) can be obtained by introducing the ansatzwith unknown (real) parameters and (). Inserting (17) into (16) yields a system of linear equations which can be easily solved:where with given by (9), while is the vector containing the unknown coefficients and . We refer to (17) as the asymptotic equilibrium of the linear problem (as dictated by the actuator), since all solutions of the linear problem approach it for . Clearly, this state will be spatially localized if the forcing frequency lies within a spectral gap; otherwise it will be spatially extended (with the former corresponding to a surface breather, which we discuss in the following section). For small driving amplitude , we find that the long-term dynamics of the system approaches an asymptotic state that is reasonably well approximated by (17). However, for larger driving amplitudes the effect of the nonlinearity becomes significant, and we can no longer rely on the linear analysis.

In order to understand the dynamics in the high-amplitude regime, we perform a bifurcation analysis of time-periodic solutions of the nonlinear problem (1) using a Newton-Raphson method that yields exact time-periodic solutions (to a prescribed tolerance) and their linear stability through the computation of Floquet multipliers (see appendix for details). We use the asymptotic linear state (17) as an initial guess for the Newton iterations.

#### 4. Main Results

Besides the linear limit considered above, another relevant situation is that of the Hamiltonian system, that is, with in (1) and in (3). It is well known that time-periodic solutions that are localized in space (i.e., breathers) exist in a host of discrete Hamiltonian systems [25], including granular crystals with spatial periodicity [11]. In particular, Hamiltonian trimer granular chains were recently studied in [14]. There, it was observed that breathers with frequency in the second spectral gap are more robust than those with frequency in the first spectral gap. Features that appear to enhance the stability of the higher gap breathers are (i) tails avoiding resonances with the spectral bands and (ii) lighter beads oscillating out-of-phase. In that sense, the breathers found in the second spectral gap of trimers are somewhat reminiscent of breathers found in the gap of dimer lattices [11]. Furthermore, breathers of damped-driven dimer granular crystals (i.e., (1) with a spatial periodicity of two) were studied recently in [13]. In that setting, the breathers become* surface breathers* since they are localized at the surface, rather than the center of the chain. Yet, if one translates the surface breather to the center of the chain, it bears a strong resemblance to a “bulk” breather. Thus, nonlinearity, periodicity, and discreteness enabled two classes of relevant states: The nonlinear surface breathers and ones tuned to the external actuator (i.e., those proximal to the asymptotic linear state dictated by the actuator). These two waveforms were observed to collide and disappear in a limit cycle saddle-node bifurcation as the actuation amplitude was increased [13]. Beyond this critical point, no stable, periodic solutions were found to exist in the dimer case of [13] and the dynamics were found to “jump” to a chaotic branch. We aim to identify similar features in the case of the trimer with a particular emphasis on the differences arising due to the higher order periodicity of the system. To that end, we first present results on surface breathers with frequency in the second gap and perform parameter continuation in driving amplitude to draw comparisons to the bifurcation structure in dimer granular lattices. Following that, we investigate breathers in the first spectral gap and compare them to their second gap breather counterparts. Finally we identify both localized and spatially extended states as the driving frequency is varied through the entire range of spectral values covering both gaps and the three pass bands between which they arise.

##### 4.1. Driving in the Second Spectral Gap

We first consider a fixed driving frequency of kHz, which lies within the second spectral gap (see Table 2). For low driving amplitude there is a single time-periodic state (i.e., the one proximal to the driven linear state), which the experimentally observed dynamics follows (see panels (a) and (b) in Figure 3). A continuation in driving amplitude of numerically exact time-periodic solutions (see the blue, red, and green lines of Figure 3(b)) reveals the existence of three branches of solutions for a range of driving amplitudes. The branch indicated by label (A) is proximal to the linear driven state given by (17) (shown as a gray dashed line in Figure 3(b)). Each solution making up this branch is in-phase with the boundary actuator (see Figure 4(a)) and has lighter masses that are out-of-phase with respect to each other (see, e.g., Figure 7(a)). These solutions are asymptotically stable, which is also evident in the experiments (see Figure 3(a) and the black markers with error bars in Figure 3(b)). At a driving amplitude of nm, an unstable and stable branch of nonlinear surface breathers arise through a saddle-node bifurcation (see labels (B) and (C) of Figure 3(b), resp.). At the bifurcation point nm, the profile strongly resembles that of the corresponding Hamiltonian breather (when shifted to the center of the chain), hence the name surface breather. It is important to note that the presence of dissipation does not allow for this branch to bifurcate near (where this bifurcation would occur in the absence of dissipation). Instead, the need of the drive to overcome the dissipation ensures that the bifurcation will emerge at a finite value of . As is increased, the solutions constituting the unstable “separatrix” branch (b) resemble progressively more the ones of branch (a); see Figure 4(b). For example, branch (b) progressively becomes in-phase with the actuator as is increased. Indeed these two branches collide and annihilate at nm. On the other hand, the (stable, at least for a parametric interval in ) nonlinear surface breather (c) is out-of-phase with the actuator (see Figure 4(c)). This solution loses its stability through a Neimark-Sacker bifurcation, which is the result of a Floquet multiplier (lying off the real line) acquiring modulus greater than unity (see label (D) of Figures 3(b) and 4(d)). Such a Floquet multiplier indicates concurrent growth and oscillatory dynamics of perturbations and thus the instability is deemed as an oscillatory one. Solutions with an oscillatory instability are marked in green in Figure 3(b), whereas red dashed lines correspond to purely real instabilities (see also Figure 4(b) as a case example) and solid blue lines denote asymptotically stable regions. The reason for making the distinction between real and oscillatory instabilities is that quasi periodicity and chaos often lurk in regimes in parameter space where solutions possess such instabilities [13]. Indeed, past the above mentioned saddle-node bifurcation, as the amplitude is further increased, for an additional narrow parametric regime, the computational dynamics (gray circles in Figure 3(c)) follows the upper nonlinear surface mode of branch (c); yet, once this branch becomes unstable the dynamics appears to reach a chaotic state. In the experimental dynamics (black squares in Figure 3(c)), a very similar pattern is observed qualitatively, although the quantitative details appear to differ somewhat. Admittedly, as more nonlinear regimes of the system’s dynamics are accessed (as increases), such disparities are progressively more likely due to the opening of gaps between the spheres and the limited applicability (in such regimes) of the simple Hertzian contact law. See also the discussion in Section 4.3 for more details on possible explanations for disagreement in high amplitude regimes.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

##### 4.2. Driving in the First Spectral Gap

We now consider a fixed driving frequency of kHz, which lies within the first spectral gap. Qualitatively, the breather solutions and their bifurcation structure are similar to those in the second gap (compare Figures 3 and 4 with Figures 5 and 6). However, there are several differences, which we highlight here. For example, the lighter masses now oscillate (nearly) in-phase, rather that out-of-phase (for comparison, see panels (a) and (b) in Figure 7). The emergence of the nonlinear surface modes occurs for a much larger value of the driving amplitude ( *μ*m) as well as the bifurcation of the state (a) dictated by the actuator with the relevant branch of the nonlinear surface modes ( *μ*m). Indeed, the latter turns out to be out-of-range for experimentally controllable driving amplitudes given the stroke of the piezoelectric actuator and the power by the electric amplifier in our experimental setup. Thus, applications such as bifurcation based rectification [12] are not suitable for breather frequencies in the first gap in the present setting. Although the range of amplitudes yielding stable solutions is much larger, these solutions are still effectively linear (see the dashed gray line in Figure 5(b)). Another peculiar feature particular to this case is that the branch (a) loses its stability through a Neimark-Sacker bifurcation ( *μ*m), rather than through a saddle-node bifurcation with the nonlinear surface mode, although it regains stability shortly thereafter ( *μ*m). Interestingly, however, the corresponding experimental branch deviates from the theoretical (near-linear) branch close to this destabilization point, apparently leading to large amplitude, chaotic behavior thereafter.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

As an additional feature worth mentioning, there are breathers within the first gap that resonate with the linear modes. For example, any breather with a frequency will have a second harmonic that lies in the second optical band. Indeed, the breather solution depicted in Figure 7(c) has a frequency , but the tail oscillates with twice that frequency (i.e., it is resonating with the linear mode at a frequency of ). Finally, we note that the magnitude of the instabilities tends to be much larger in the first spectral gap, due possibly to their spatial structure (see, e.g., Figure 6(d)), in agreement with what was found in [14].

While a detailed probing of the spatial profiles and the corresponding bifurcation structure in the first and second gap reveal several differences, a common theme is that the experimental data points generally follow rather accurately the theoretically predicted curve. In particular, the agreement between experiment and theory is rather satisfactory as long as the data points are within the stable parametric regions. This is clearly depicted in panel (b) of Figures 3 and 5. Note that the linear asymptotic equilibrium (shown by the dashed gray line and corresponding to ) becomes closer to the theoretically predicted curve for small relative displacements (see, e.g., (5)), a feature that we also use as a consistency check for the numerics. In both cases (first or second gap) once the amplitude of the actuator reaches a critical value ( *μ*m for the second gap and *μ*m for the first gap), the experimental data experience jumps, seemingly leading to chaotic behavior. Importantly, these jumps arise close to the destabilization points of the respective branches. Nevertheless, it should also be pointed out that, in some of these regimes (especially so in the case of Figure 5), these driving amplitudes may be near or beyond the regime of (accurately) controllable experimentally accessible amplitudes. See the discussion in Section 4.3 for more details on possible explanations for disagreement in high amplitude regimes.

##### 4.3. A Full Driving Frequency Continuation

Our previous considerations suggest that the driving frequency plays an important role in the observed dynamics of the time-periodic solutions. Figure 8 complements those results by showing a frequency continuation for a fixed driving amplitude *μ*m. For this driving amplitude, the response is highly nonlinear, where the ratio of dynamic to static compression ranges as high as, for example, . In the low amplitude (i.e., near linear) range, the peaks in the force response coincide with the eigenfrequency spectrum presented in Figure 2(a), with all solutions being stable. For the chosen value of driving amplitude, the features of the diagram at a low driving frequency are similar to their linear counterparts (see, e.g., the smooth resonant peaks in the acoustic band of Figure 8(a)). Similar to Figures 3 and 5, the experimentally obtained values match the numerically obtained, dynamically stable response quite well, with a few notable exceptions: (i) in narrow parametric regions around kHz and kHz where the experimental points experience jumps and (ii) the experimentally observed velocities for the same frequency are upshifted when compared with the theoretical curve (see, e.g., the region kHz in Figure 8(a) where the peaks do not align exactly). Possible explanations for such discrepancies are given below, but we also note that discrepancies in experimental and theoretical dispersion curves of granular crystals have also been previously discussed; see, for a recent example, [45].

**(a)**

**(b)**

In general, as the driving amplitude is increased, the peaks in the bifurcation diagram bend (in some cases, they bend into the gap) in a way reminiscent of the familiar fold-over event of driven oscillators; see, for example, [46] for a recent study extending this to chains. This, in turn, causes a series of intricate bifurcations and a loss of stability (see, e.g., the second gap in Figure 8(a)). In particular, large regions of oscillatory instabilities are born, and hence, we expect regions of quasi periodicity and of potentially chaotic response (see, e.g., Section 4.1). Not surprisingly, the experimentally observed values depart from the theoretical prediction once stability is lost. However, the general rising and falling pattern are captured qualitatively, see Figure 8(b). Nevertheless, it is worthwhile to note that the latter figure contains a very elaborate set of loops and a high multiplicity of corresponding (unstable) periodic orbit families for which no clear physical intuition is apparent at present. More troubling, however, is the fact that there are also some stable regions where the experimental points are off the theoretical curve. A systematic study using numerical simulations revealed that, within these narrow stable regions at hand (see, e.g., Figure 8(b)), a large number of periods are required in order to converge to the exact periodic solution. In particular, it happens that the Floquet multipliers are inside the unit circle, but also very close to it, suggesting that the experimental time window of 50 ms used is not sufficiently long to observe experimentally the periodic structures predicted theoretically by our model. Thus, the reported experimental observations in this regime appear to be capturing a transient stage and not the eventual asymptotic form of the dynamical profile. A final comment on the possible discrepancy between the model and experiments is due to the fact that in this highly nonlinear regime (of large driving amplitudes), the beads come out of contact (the dynamic force exceeds the static force), in which case other dynamic effects, such as particles’ rotations, could have a significant role. Such features are not described in our simple model (1). In fact, the contact points of the particles cannot be perfectly aligned in experiments due to the clearance between the beads and the support rods. This is likely to result in dynamic buckling of the chain under the strong excitations, which are considered in some of the cases above. Here also, it is plausible that additional forms of dissipation affect the particle motions, leading the experimental findings (in such settings) to underpredict the velocities observed in direct numerical simulations. Nonetheless, the overall experimental trends of bifurcation in Figure 8(a) are in fair qualitative agreement with the theoretical predictions.

#### 5. Conclusions and Future Challenges

We have presented a natural extension of earlier considerations of (i) Hamiltonian breathers in trimer chains and (ii) surface breathers in damped-driven dimers, by studying a damped-driven trimer granular crystal through a combination of experimental and theoretical/computational approaches. We found that the breathers with frequency in the second gap are in analogy with those in the sole gap of the damped-driven dimer, in which the interplay of damped-driven dynamics with nonlinearity and spatial discreteness gives rise to saddle-node bifurcations of time-periodic solutions and turning points beyond which there is no stable ordered dynamics as well as surface modes and generally rather complex and tortuous bifurcation diagrams. While similar structures are found in the first gap, they appear to be less robust (given the magnitude of their instabilities) and can resonate with the higher-order linear bands, a feature interesting in its own right. A continuation in driving frequency revealed that the nonlinearity causes the resonant peaks to bend, possibility into a spectral gap. However, this nonlinear bending also causes the solutions to lose stability in various ways, such as Neimark-Sacker and saddle-node bifurcations. More importantly, all of these features were validated experimentally through laser Doppler vibrometry, allowing, for the first time to our knowledge, the full-field visualization of surface breathers in granular crystals.

Several interesting questions remain unexplored, including, for example, the mechanism leading to the emergence of the observed chaotic dynamics and the more accurate quantitative description of that regime. We also observed that the bifurcation structure is generally more complex in the case of breathers that are resonant with the linear modes, which bears further probing. Other possible avenues for future research include the investigation of dark breathers in such chains. Recall that, as in [15], dark breathers may bifurcate from the top edge of the acoustic, first, or second optical band, under suitable drive. Lastly, it would be quite interesting to generalize considerations of breathers and surface modes in two-dimensional hexagonal lattices with both homogeneous and heterogeneous compositions. Such studies are currently under consideration and will be reported in future publications.

#### Appendix

#### Numerical Methods and Linear Stability

In this Appendix we shortly discuss the numerical methods employed for finding exact time-periodic solutions of (1) together with the parametric continuation techniques for identifying families of such solutions as either the value of the actuation frequency (or, equivalently, the period ) or amplitude changes. The interested reader can find more information along these directions, for example, in [47, 48] (among others) and references therein.

In order to find time-periodic solutions of (1), we convert the latter into a system of first order ODEs written in compact formwith , while and represent the -dimensional position and velocity vectors, respectively. Next, we define the Poincaré map: where is the initial condition and is the result of integrating (A.1) forward in time until (using, e.g., the DOP853 [49] or ODE [50] time integrators). Then, a periodic solution with period (i.e., satisfying the property ) will be a root of the map . To this end, Newton’s method is applied to the map and thus yields the following numerical iteration schemetogether with the Jacobian of the map , namely, where is the identity matrix; is the solution to the variational problem: with initial data , while is the Jacobian of the equations of motion (A.1) evaluated at the point .

Subsequently, the iteration scheme (A.2) is fed by a suitable initial guess (for fixed values of the actuation amplitude and frequency ) and applied until a user-prescribed tolerance criterion is satisfied. Thus, a time-periodic solution is obtained upon successful convergence (and with high accuracy) together with the eigenvalues or Floquet multipliers (denoted by ) of the* monodromy matrix * which convey important information about the stability of the solution in question. In particular, a periodic solution is deemed asymptotically stable (or a stable limit cycle) if there are no Floquet multipliers outside the unit circle whereas a periodic solution with Floquet multipliers lying outside the unit circle is deemed to be linearly unstable. Note that the iteration scheme (A.2) can be initialized using the linear asymptotic equilibrium given by (17). Alternatively, one may start with zero initial data and directly integrate the equations of motion (A.1) forward in time. Then, (A.2) can be initialized with the output waveform of the time integrator at a particular time.

Having obtained an exact periodic solution to (1) for given values of and , we perform parametric continuations over these parameters using the computer software AUTO [51], which employs a pseudo-arclength continuation and thus enables the computation of branches past turning points. This way, we are able to trace entire families of periodic solutions and study their corresponding stability characteristics in terms of the Floquet multipliers.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

E. G. Charalampidis gratefully acknowledges financial support from the FP7-People IRSES-606096: “*Topological Solitons, from Field Theory to Cosmos*”. P. G. Kevrekidis acknowledges support from the National Science Foundation (NSF) under Grants CMMI-1000337 and DMS-1312856, from the ERC and FP7-People under Grant IRSES-606096, and from the US-AFOSR under Grant FA9550-12-1-0332. P. G. Kevrekidis’s work at Los Alamos is supported in part by the U.S. Department of Energy. C. Chong was partially supported by the ETH Zurich Foundation through the Seed Project ESC-A 06-14. F. Li and J. Yang thank the support of the NSF (CMMI-1414748) and the US-ONR (N000141410388). E. G. Charalampidis and C. Chong thank M. O. Williams (PACM, Princeton University) for insight regarding the AUTO codes used for the bifurcation analysis performed in this work.