Time-Periodic Solutions of Driven-Damped Trimer Granular Crystals

In this work, we consider time-periodic structures of trimer granular crystals consisting of alternate chrome steel and tungsten carbide spherical particles yielding a spatial periodicity of three. 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 become 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.

In this work, we consider time-periodic structures of trimer granular crystals consisting of alternate chrome steel and tungsten carbide spherical particles yielding a spatial periodicity of three. 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 become 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.
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, e.g., in [24]. The broad and diverse span of fields where such structures have been of interest includes, among others, optical waveguide arrays or photorefractive crystals [25], micromechanical cantilever arrays [26], Josephson-junction ladders [27], layered antiferromagnetic crystals [28], halide-bridged transition metal complexes [29], dynamical models of the DNA double strand [30] and Bose-Einstein condensates in optical lattices [31].
In Fermi-Pasta-Ulam type settings (which are intimately connected to the realm of precompressed granular crystals), it was proven in [32,33] (see also the discussion in [24]) that small amplitude discrete breathers are absent in spatially homogeneous (i.e. monoatomic) chains. Instead, dark such states (those on top of a non-vanishing background) have been found therein [13,14]. 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 [8,10] (see also for recent experiments [34]) and a spatial periodicity of two (i.e. dimer lattices) [9,11]. Further recent experimental works explored solitary waves in trimers and higher periodicity chains [35,36], as well as the linear dispersion properties of such chains [37]. Motivated by these works, a theoretical study of breathers in granular crystals with higher order spatial periodicity (such as trimers and quadrimers) FIG. 1: (Color online) Experimental setup of a 21-particle granular chain composed of chrome steel and tungsten carbide beads in a 2:1 ratio. Actuation and sensing systems based on a piezoelectric force sensor and a laser Doppler vibrometer are also illustrated.
was recently conducted in [12]. Therein, it was demonstrated that breathers with a frequency in the highest gap appear to be more robust than their counterparts with frequency in the lower gaps.
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 [9][10][11] 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.
The paper is structured as follows. In Secs. II and III, we describe the experimental and theoretical setups respectively. The main results are presented in Sec. IV, 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, Sec. V provides our conclusions and discusses a number of future challenges. 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 N = 21 spherical beads. The beads are made out of chrome steel (S: gray particles in Fig. 1) and tungsten carbide (W: black particles) materials. See Table I 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 Fig. 1). We measure the pre-applied static force (F 0 = 10 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 mass.

II. EXPERIMENTAL SETUP
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.

III. THEORETICAL SETUP
The equation we use to model the experimental setup is a Fermi-Pasta-Ulam-type lattice with a Hertzian potential [1] leading to: where n = 1, . . . , N , u n = u n (t) ∈ R is the displacement of the n-th bead from its equilibrium position at time t ∈ R, M n is the mass of the n-th bead, δ 0,n is a precompression factor induced by the static force F 0 = A n δ 3/2 0,n and the bracket is defined by [x] + = max(0, x). The 3/2 power accounting for the nonlinearity of the model is a result of the sphere-to-sphere contact, i.e., the so-called Hertzian contact [38]. The form of the dissipation is a dash-pot, which has been utilized in the context of granular crystals in several previous works [10,14]. The strength of the dissipation is captured by the parameter τ , which serves as the sole parameter used to fit experimental data (τ = 2.1 ms in this study). The elastic coefficient A n depends on the interaction of bead n with bead n + 1 and for spherical point contacts has the form [39] where E n , ν n and r n are the Young's modulus, Poisson's ratio and the radius, respectively, of the n-th bead. The left boundary is an actuator and the right one is kept fixed, i.e., where a and f b represent the driving amplitude and frequency, respectively. Note that at the boundaries (i.e., n = 0 and n = N + 1), we consider flat surface-sphere contacts while the same material properties are assumed; thus, e.g., the elastic coefficient at the left boundary takes the form .

A. Linear regime and dispersion relation
Assuming dynamic strains that are small relative to the static precompression, i.e., Eq. (1) can be well approximated by its linearized form with K n = 3/2 A n δ 1/2 0,n corresponding to the linearized stiffnesses. Subsequently, Eq. (6) can be converted into a system of 2N first-order equations and written conveniently in matrix form aṡ with Y = u 1 , . . . , u N , v 1 (=u 1 ) , . . . , v N (=u N ) where O and I represent the N × N zero and identity matrices, respectively, and the N × N (tridiagonal) matrix C is given by Note that in the above linearization we applied fixed boundary conditions at both ends of the chain, i.e., u 0 = u N +1 = 0 (we account for the actuator in the following subsection). This equation is solved by Y n = y n e iωt , where ω corresponds to the angular frequency (with ω = 2πf ) and where, with (λ, y) corresponding to the eigenvalue-eigenvector pair, while λ = iω and y = (y 1 , . . . , y 2N ) T . The eigenfrequency spectrum f = −i λ 2π using the values of Table I and τ = 2.1 ms is shown in Fig. 2(a). On the other hand, the dispersion relation for the trimer chain configuration within an infinite lattice can be obtained using a Bloch wave ansatz [37]. We first re-write Eq. (6) in the following convenient form with j ∈ Z + , where we made use of the spatial periodicity of the trimer lattice, i.e., The Bloch wave ansatz has the form, where u 0 , v 0 and w 0 are the wave amplitudes, k is the wavenumber and α is the size (or the equilibrium length) of one unit cell of the lattice. Substitution of Eqs. (12) into Eqs. (11) yields     The non-zero 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 form which has (non-trivial) solutions at k = 0 and k = π/α (i.e., first Brillouin zone) given by (second lower optic), 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 theu n /τ term) in order to model the experimental setup. Therefore, the reported frequencies hereafter will correspond to the real part of these values, i.e., the part contributing to the dispersion properties of the solutions rather than their decay. In Table II, the predicted values of the cut-off frequencies are presented (up to two decimal places) using the values of the material parameters of Table I and τ = 2.1 ms. Finally, upon solving numerically the dispersion relation (cf. Eq. (14)), the frequency as a function of the wavenumber (κα) is presented in Fig. 2(b) together with the cut-off frequencies of Table II. Note that these cut-off frequencies are also plotted as horizontal dashed black lines in Fig. 2(a) for comparison. It can be discerned from both panels of Fig. 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, We also note the presence of two additional localized modes depicted in Fig. 2(a) between the first and second optical pass bands due to the fixed boundary conditions. These modes, denoted by black dots in Fig. 2(a), correspond to surface modes and their existence has been previously discussed, e.g., in [9].

B. 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 Eq. (7) in the form,Ẏ where the sole non-zero entry of F(ω b ) is at the (N + 1)th node and has the form F N +1 = a cos (ω b t) and the matrix A is given by Eq. (8). Equation (16) can be solved by introducing the ansatz with unknown (real) parameters a j and b j (j = 1, · · · , N ). Inserting Eq. (17) into Eq. (16) yields a system of linear equations which can be easily solved: where is the vector containing the unknown coefficients and F = a K0 M1 , 0, . . . , 0 . We refer to Eq. (17) as the asymptotic equilibrium of the linear problem (as dictated by the actuator), since all solutions of the linear problem approach it for t → ∞. Clearly, this state will be spatially localized if the forcing frequency f b 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 a, 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 timeperiodic 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 A for details). We use the asymptotic linear state (17) as an initial guess for the Newton iterations.

IV. MAIN RESULTS
Besides the linear limit considered above, another relevant situation is that of the Hamiltonian system, i.e., with 1/τ → 0 in Eq. (1) and a → 0 in Eq. (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 [24], including granular crystals with spatial periodicity [9]. In particular, Hamiltonian trimer granular chains were recently studied in [12]. 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 [9]. Furthermore, breathers of damped-driven dimer granular crystals (i.e., Eq. (1) with a spatial periodicity of two) were studied recently in [11]. 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 [11]. Beyond this critical point, no stable, periodic solutions were found to exist in the dimer case of [11] 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.

A. Driving in the second spectral gap
We first consider a fixed driving frequency of f b = 6.3087 kHz, which lies within the second spectral gap (see Table II). 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 Fig. 3). A continuation in driving amplitude of numerically exact time-periodic solutions (see the blue, red and green lines of Fig. 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 Eq. (17) (shown as a gray dashed line in Fig. 3(b)). Each solution making up this branch is in-phase with the boundary actuator (see Fig. 4(a)), and has lighter masses that are out-of-phase with respect to each other (see e.g., Fig. 7(a)). These solutions are asymptotically stable, which is also evident in the experiments (see Fig. 3(a) and the black markers with error bars in Fig. 3(b)). At a driving amplitude of a ≈ 22.75 nm, an unstable and stable branch of nonlinear surface breathers arise through a saddlenode bifurcation (see labels (b) and (c) of Fig. 3(b), respectively). At the bifurcation point a ≈ 22.75 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 a → 0 (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 a. As a is increased, the solutions constituting the unstable "separatrix" branch (b) resemble progressively more the ones of branch (a); see Fig. 4(b). For example, branch (b) progressively becomes in-phase with the actuator as a is increased. Indeed these two branches collide and annihilate at a ≈ 32.26 nm. On the other hand, the (stable, at least for a parametric interval in a) nonlinear surface breather (c) is out-of-phase with the actuator (see Fig. 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 Fig. 3(b) and Fig. 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 Fig. 3(b), whereas red dashed lines correspond to purely real instabilities (see also Fig. 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 [11]. Indeed, past the above saddle-node bifurcation, as the amplitude is further increased, for an additional narrow parametric regime, the computational dynamics (gray circles in Fig. 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 Fig. 3(c)), a very similar pattern is observed qualitatively, although the quantitative details appear to somewhat differ. Admittedly, as the more nonlinear regimes of the system's dynamics are accessed (as a 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. is similar to the one presented in Fig. 3(b), which was for a driving frequency in the second gap (f b = 6.3087 kHz). However, the bifurcation points occur for much larger values of the driving amplitude and the structure of the solutions themselves varies considerably in comparison to the second gap (see text). (c) Zoomed out version of panel (b). In this case, the jump of the experimental data is likely due to driving out of a controllable range, rather than to chaos, as in the case of the second gap breathers.

B. Driving in the first spectral gap
We now consider a fixed driving frequency of f b = 3.7718 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 Figs. 3 and 4 with Figs. 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 Fig. 7). The emergence of the nonlinear surface modes occurs for a much larger value of the driving amplitude (a ≈ 0.1778 µm) as well as the bifurcation of the state (a) dictated by the actuator with the relevant branch of the nonlinear surface modes (a ≈ 1.3 µ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 [10] 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 Fig. 5(b)). Another peculiar feature particular to this case is that the branch (a) loses its stability through a Neimark-Sacker bifurcation (a ≈ 0.9 µm), rather than through a saddle-node bifurcation with the nonlinear surface mode, although it regains stability shortly thereafter (a ≈ 0.9876 µ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.
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 f b ∈ [3.3550, 3.6850] will have a second harmonic that lies in the second optical band. Indeed, the breather solution depicted in Fig. 7(c) has a frequency f b = 3.4429, but the tail oscillates with twice that frequency (i.e., it is resonating with the linear mode at a frequency of 2f b ). 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., Fig. 6(d)), in agreement with what was found in [12].
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 Figs. 3 and 5. Note that the linear asymptotic equilibrium (shown by the dashed gray line) coincides with the theoretically predicted curve as soon as Eq. (5) holds (a feature that we also use as a consistency check for the numerics). In both cases (first or second gap) once the amplitude a of the actuator reaches a critical value (a ≈ 0.32016 µm for the second gap and a ≈ 0.82416 µm for the first gap), the experimental data experience jumps, possibly 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 Fig. 5), these driving amplitudes may be near or beyond the regime of (accurately) controllable experimentally accessible amplitudes.

C. 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 a = 0.24816 µm. For this driving amplitude, the response is highly nonlinear, where the ratio of dynamic to static compression ranges as high as e.g. |u n − u n+1 |/δ 0,n = 1.5. In the low amplitude (i.e., near linear) range, the peaks in the force response coincide with the eigenfrequency spectrum presented in     Fig. 8(a)). Similar to Figs. 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 f b ≈ 2.35 kHz and f b ≈ 2.7 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 f b ≈ 1 − 5.2 kHz in Fig. 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, e.g., for a recent example [37].
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 the driven oscillators; see, e.g., [40] 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 Fig. 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. Sec. IV A). Not surprisingly, the experimentally observed values depart from the theoretical prediction once stability is lost. However, the general rising and falling pattern is captured qualitatively, see Fig. 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. Fig. 8(b)), a large number of periods is required in order to converge to the exact periodic solution. In particular, it happens that the Floquet multipliers are inside the unit circle but very close to unity, 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, 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, which 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. Nonetheless, the overall experimental trends of bifurcation in Fig. 8(a) are in fair qualitative agreement with the theoretical predictions.

V. CONCLUSIONS & 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 that 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 to surface modes and generally rather complex and tortuous bifurcation diagrams. While similar structures are found in the 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. We also observed that the bifurcation structure is generally more complex in the case of breathers that 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 [13,14], 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 (see e.g. [41] for a recent example of a study on traveling waves in such settings) with both homogeneous, as well as heterogeneous compositions. Such studies are currently under consideration and will be reported in future publications.