# Time-Periodic Solutions of Driven-Damped Trimer Granular Crystals

Authors: E. G. Charalampidis; F. Li; C. Chong; J. Yang; P. G. Kevrekidis
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 Sections2 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

Figure1 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 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 (F0=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 masses.Table 1
Parameter values of the trimer granular chain (S-W-S).
Material
Parameter
E [N/m2]
r [m]
ν
ρ [kg/m3]
Chrome steel (S)
200
×
1
0
9
9.525
×
1
0
-
3
0.3
7780
Tungsten carbide (W)
628
×
1
0
9
9.525
×
1
0
-
3
0.28
14980Figure 1
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.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/Hz1/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 to(1)u¨n=An-1Mnδ0,n-1+un-1-un+3/2-AnMnδ0,n+un-un+1+3/2-u˙nτ,where n=1,…,N, un=un(t)∈R is the displacement of the nth bead from its equilibrium position at time t∈R, Mn is the mass of the nth bead, δ0,n is a precompression factor induced by the static force F0=Anδ0,n3/2, 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, 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 (τ=2.1ms in this study). The elastic coefficient An depends on the interaction of bead n with bead n+1 and for spherical point contacts has the form [44](2)An=4EnEn+1rnrn+1/rn+rn+13En+11-νn2+3En1-νn+12,where En, νn, and rn are the Young’s modulus, Poisson’s ratio, and the radius, respectively, of the nth bead. The left boundary is an actuator and the right one is kept fixed; that is,(3)u0t=acos⁡2πfbt,uN+1t=0,where a and fb 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, for example, the elastic coefficient at the left boundary takes the form(4)A0=2E1r131-ν12.


### 3.1. Linear Regime and Dispersion Relation

Assuming the dynamic strains that are small relative to the static precompression, that is,(5)un-un+1δ0,n≪1,equation (1) can then be well approximated by its linearized form(6)u¨n=Kn-1Mnun-1-un-KnMnun-un+1-u˙nτ,with Kn=3/2Anδ0,n1/2 corresponding to the linearized stiffnesses. Subsequently, (6) can be converted into a system of 2N first-order equations and written conveniently in matrix form as(7)Y˙=AY,with Y=(u1,…,uN,v1=u˙1,…,vN=u˙N)T and(8)A=OIC-1/τI,where O and I represent the N×N zero and identity matrices, respectively, and the N×N (tridiagonal) matrix C is given by(9)
C
=
-
K
0
+
K
1
M
1
K
1
M
1
0
⋯
0
K
1
M
2
-
K
1
+
K
2
M
2
K
2
M
2
⋮
0
⋱
⋱
⋱
0
⋮
K
N
-
2
M
N
-
1
-
K
N
-
2
+
K
N
-
1
M
N
-
1
K
N
-
1
M
N
-
1
0
⋯
0
K
N
-
1
M
N
-
K
N
-
1
+
K
N
M
N
.Note that in the above linearization we applied fixed boundary conditions at both ends of the chain; that is, u0=uN+1=0 (we account for the actuator in the following subsection). This equation is solved by Yn=yneiωt, where ω corresponds to the angular frequency (with ω=2πf) and where(10)Ay=λy,with λ,y corresponding to the eigenvalue-eigenvector pair, while λ=iω and y=y1,…,y2NT. The eigenfrequency spectrum f=-i(λ/2π) using the values of Table 1 and τ=2.1 ms is shown in Figure 2(a).Figure 2
Eigenfrequency spectrum of the trimer granular chain for parameter values given in Table1 and τ=2.1ms. (a) Spectrum corresponding to (10) for N=21 beads. Note that the horizontal dashed black lines correspond to the cut-off frequencies listed in Table 2. (b) Numerically obtained frequencies of the infinite lattice using (14) as a function of the wavenumber (kα).
(a)
(b)Table 2
Calculated cut-off frequencies (in kHz) using (15).
f
c
,
5
f
c
,
4
f
c
,
3
f
c
,
2
f
c
,
6
3.04
4.23
6.04
6.71
7.37On 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 form(11a)u¨3j-2=K3M1u3j-3-u3j-2-K1M1u3j-2-u3j-1-u˙3j-2τ,(11b)u¨3j-1=K1M2u3j-2-u3j-1-K1M2u3j-1-u3j-u˙3j-1τ,(11c)u¨3j=K1M1u3j-1-u3j-K3M1u3j-u3j+1-u˙3jτ,with j∈Z+, where we made use of the spatial periodicity of the trimer lattice; that is, M3j=M3j-2≡M1(=M3) and K3j-2=K3j-1≡K1(=K2) and M3j-1≡M2 and K3j≡K3. The Bloch wave ansatz has the form(12a)u3j-2=η1eikαj+ωt,(12b)u3j-1=η2eikαj+ωt,(12c)u3j=η3eikαj+ωt,where η1, η2, and η3 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 (12a), (12b), and (12c) into (11a), (11b), and (11c) yields(13)
ω
2
-
K
1
+
K
3
M
1
-
i
ω
τ
K
1
M
1
K
3
M
1
e
-
i
k
α
K
1
M
2
ω
2
-
2
K
1
M
2
-
i
ω
τ
K
1
M
2
K
3
M
1
e
i
k
α
K
1
M
1
ω
2
-
K
1
+
K
3
M
1
-
i
ω
τ
η
1
η
2
η
3
=
0
0
0
.The 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 form(14)
M
1
2
M
2
ω
6
-
i
3
M
1
2
M
2
τ
ω
5
-
M
1
2
K
3
M
2
+
K
1
M
1
+
M
2
+
3
M
1
M
2
τ
2
ω
4
+
i
M
1
τ
4
K
3
M
2
+
K
1
M
1
+
M
2
+
M
1
M
2
τ
2
ω
3
-
4
K
1
2
K
3
sin
2
⁡
k
α
2
+
K
1
K
1
+
2
K
3
2
M
1
+
M
2
+
2
M
1
τ
2
K
3
M
2
+
K
1
M
1
+
M
2
ω
2
-
i
K
1
τ
K
1
+
2
K
3
2
M
1
+
M
2
ω
=
0
,which has (nontrivial) solutions at k=0 and k=π/α (i.e., first Brillouin zone) given by(15)fc,1k=0=i12πτlower acoustic,fc,2k=0=±12πK1+2K3M1-14τ2+i14πτsecond lower optic,fc,3k=0=±12πK12M1+M2M1M2-14τ2+i14πτfirst upper optic,fc,4k=πα=±12πK1M1-14τ2+i14πτfirst lower optic,fc,5k=πα=±12πK12M1+M2+2K3M2-f~2M1M2-14τ2+i14πτupper acoustic,fc,6k=πα=±12πK12M1+M2+2K3M2+f~2M1M2-14τ2+i14πτsecond upper optic,where f~=K12M1+M2+2K3M22-16K1K3M1M2. The solutions (15) correspond to thecut-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 u˙n/τ term) in order to model the experimental setup. Therefore, the reported frequencies hereafter will correspond to thereal 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 τ=2.1ms.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, [fc,5,fc,4] and [fc,3,fc,2], together with one semi-infinite gap [fc,6,∞). 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 form(16)Y˙=AY+Ft,where the sole nonzero entry of F is at the (N+1)th node and has the form FN+1=a(K0/M1)cosωbt and the matrix A is given by (8). An inhomogeneous solution of (16) can be obtained by introducing the ansatz(17)Y=a1cos⁡ωbt+b1sin⁡ωbt⋮aNcos⁡ωbt+bNsin⁡ωbt-a1ωbsin⁡ωbt+b1ωbcos⁡ωbt⋮-aNωbsin⁡ωbt+bNωbcos⁡ωbt,with unknown (real) parameters aj and bj (j=1,…,N). Inserting (17) into (16) yields a system of linear equations which can be easily solved:(18)A~X=F~,where(19)A~=C~ωτI-ωτIC~,C~=-C+ωb2I with C given by (9), while X=a1,…,aN,b1,…,bN is the vector containing the unknown coefficients and F~=a(K0/M1),0,…,0. 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 t→∞. Clearly, this state will be spatially localized if the forcing frequency fb 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 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.


## 3.1. Linear Regime and Dispersion Relation

Assuming the dynamic strains that are small relative to the static precompression, that is,(5)un-un+1δ0,n≪1,equation (1) can then be well approximated by its linearized form(6)u¨n=Kn-1Mnun-1-un-KnMnun-un+1-u˙nτ,with Kn=3/2Anδ0,n1/2 corresponding to the linearized stiffnesses. Subsequently, (6) can be converted into a system of 2N first-order equations and written conveniently in matrix form as(7)Y˙=AY,with Y=(u1,…,uN,v1=u˙1,…,vN=u˙N)T and(8)A=OIC-1/τI,where O and I represent the N×N zero and identity matrices, respectively, and the N×N (tridiagonal) matrix C is given by(9)
C
=
-
K
0
+
K
1
M
1
K
1
M
1
0
⋯
0
K
1
M
2
-
K
1
+
K
2
M
2
K
2
M
2
⋮
0
⋱
⋱
⋱
0
⋮
K
N
-
2
M
N
-
1
-
K
N
-
2
+
K
N
-
1
M
N
-
1
K
N
-
1
M
N
-
1
0
⋯
0
K
N
-
1
M
N
-
K
N
-
1
+
K
N
M
N
.Note that in the above linearization we applied fixed boundary conditions at both ends of the chain; that is, u0=uN+1=0 (we account for the actuator in the following subsection). This equation is solved by Yn=yneiωt, where ω corresponds to the angular frequency (with ω=2πf) and where(10)Ay=λy,with λ,y corresponding to the eigenvalue-eigenvector pair, while λ=iω and y=y1,…,y2NT. The eigenfrequency spectrum f=-i(λ/2π) using the values of Table 1 and τ=2.1 ms is shown in Figure 2(a).Figure 2
Eigenfrequency spectrum of the trimer granular chain for parameter values given in Table1 and τ=2.1ms. (a) Spectrum corresponding to (10) for N=21 beads. Note that the horizontal dashed black lines correspond to the cut-off frequencies listed in Table 2. (b) Numerically obtained frequencies of the infinite lattice using (14) as a function of the wavenumber (kα).
(a)
(b)Table 2
Calculated cut-off frequencies (in kHz) using (15).
f
c
,
5
f
c
,
4
f
c
,
3
f
c
,
2
f
c
,
6
3.04
4.23
6.04
6.71
7.37On 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 form(11a)u¨3j-2=K3M1u3j-3-u3j-2-K1M1u3j-2-u3j-1-u˙3j-2τ,(11b)u¨3j-1=K1M2u3j-2-u3j-1-K1M2u3j-1-u3j-u˙3j-1τ,(11c)u¨3j=K1M1u3j-1-u3j-K3M1u3j-u3j+1-u˙3jτ,with j∈Z+, where we made use of the spatial periodicity of the trimer lattice; that is, M3j=M3j-2≡M1(=M3) and K3j-2=K3j-1≡K1(=K2) and M3j-1≡M2 and K3j≡K3. The Bloch wave ansatz has the form(12a)u3j-2=η1eikαj+ωt,(12b)u3j-1=η2eikαj+ωt,(12c)u3j=η3eikαj+ωt,where η1, η2, and η3 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 (12a), (12b), and (12c) into (11a), (11b), and (11c) yields(13)
ω
2
-
K
1
+
K
3
M
1
-
i
ω
τ
K
1
M
1
K
3
M
1
e
-
i
k
α
K
1
M
2
ω
2
-
2
K
1
M
2
-
i
ω
τ
K
1
M
2
K
3
M
1
e
i
k
α
K
1
M
1
ω
2
-
K
1
+
K
3
M
1
-
i
ω
τ
η
1
η
2
η
3
=
0
0
0
.The 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 form(14)
M
1
2
M
2
ω
6
-
i
3
M
1
2
M
2
τ
ω
5
-
M
1
2
K
3
M
2
+
K
1
M
1
+
M
2
+
3
M
1
M
2
τ
2
ω
4
+
i
M
1
τ
4
K
3
M
2
+
K
1
M
1
+
M
2
+
M
1
M
2
τ
2
ω
3
-
4
K
1
2
K
3
sin
2
⁡
k
α
2
+
K
1
K
1
+
2
K
3
2
M
1
+
M
2
+
2
M
1
τ
2
K
3
M
2
+
K
1
M
1
+
M
2
ω
2
-
i
K
1
τ
K
1
+
2
K
3
2
M
1
+
M
2
ω
=
0
,which has (nontrivial) solutions at k=0 and k=π/α (i.e., first Brillouin zone) given by(15)fc,1k=0=i12πτlower acoustic,fc,2k=0=±12πK1+2K3M1-14τ2+i14πτsecond lower optic,fc,3k=0=±12πK12M1+M2M1M2-14τ2+i14πτfirst upper optic,fc,4k=πα=±12πK1M1-14τ2+i14πτfirst lower optic,fc,5k=πα=±12πK12M1+M2+2K3M2-f~2M1M2-14τ2+i14πτupper acoustic,fc,6k=πα=±12πK12M1+M2+2K3M2+f~2M1M2-14τ2+i14πτsecond upper optic,where f~=K12M1+M2+2K3M22-16K1K3M1M2. The solutions (15) correspond to thecut-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 u˙n/τ term) in order to model the experimental setup. Therefore, the reported frequencies hereafter will correspond to thereal 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 τ=2.1ms.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, [fc,5,fc,4] and [fc,3,fc,2], together with one semi-infinite gap [fc,6,∞). 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 form(16)Y˙=AY+Ft,where the sole nonzero entry of F is at the (N+1)th node and has the form FN+1=a(K0/M1)cosωbt and the matrix A is given by (8). An inhomogeneous solution of (16) can be obtained by introducing the ansatz(17)Y=a1cos⁡ωbt+b1sin⁡ωbt⋮aNcos⁡ωbt+bNsin⁡ωbt-a1ωbsin⁡ωbt+b1ωbcos⁡ωbt⋮-aNωbsin⁡ωbt+bNωbcos⁡ωbt,with unknown (real) parameters aj and bj (j=1,…,N). Inserting (17) into (16) yields a system of linear equations which can be easily solved:(18)A~X=F~,where(19)A~=C~ωτI-ωτIC~,C~=-C+ωb2I with C given by (9), while X=a1,…,aN,b1,…,bN is the vector containing the unknown coefficients and F~=a(K0/M1),0,…,0. 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 t→∞. Clearly, this state will be spatially localized if the forcing frequency fb 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 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, with1/τ→0 in (1) and a→0 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 becomesurface 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 offb=6.3087kHz, 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 a≈22.75 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 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 Figure 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 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 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. See also the discussion in Section 4.3 for more details on possible explanations for disagreement in high amplitude regimes.Figure 3
Comparison between experimental and theoretical results for a driving frequencyfb=6.3087 kHz, which lies in the second spectral gap. (a) Profile of solution for a driving amplitude of a=32.16 nm. Experimental values (green shaded areas) are shown as a superposition of the velocities measured over the 50 ms time window of recorded data. The numerically predicted extrema (over one period) of the exact solution obtained via Newton’s method are shown with blue (open) circles. (b) The quantity v4 as a function of the amplitude a of the drive. For the numerically exact solutions, v4 is the maximum velocity of the 4th bead over one period of motion. Note that the blue segments correspond to stable parametric regions while the red and green ones correspond to real and oscillatory unstable parametric regions, respectively. Black squares with error bars represent the experimentally measured mean values of the quantity v4 over three experimental runs, where v4 was calculated as the maximum measured over the 50 ms time window of recorded data. Finally, the dashed gray line is the maximum velocity (over one period) of bead 4 as predicted by the linear solution of (17). (c) Zoomed out version of panel (b). The gray circles correspond to numerical results that were obtained by dynamically running the equations of motion and monitoring the maximum velocity v4 over the experimental time window. The experimental data in this panel are shown again by black squares.
(a)
(b)
(c)Figure 4
Displacement profiles and corresponding Floquet multipliers of time-periodic solutions obtained at a driving frequency offb=6.3087 kHz for a driving amplitude of (a) a=0.2 nm, (b) a=30 nm, (c) a=40 nm, and (d) a=0.1 μm. These solutions correspond to the (A)–(D) labels of Figure 3(b). The lighter masses (i.e., the Steel beads) oscillate (nearly) out-of-phase and are shown as blue diamonds whereas the heavier masses (the Tungsten Carbide ones) are shown as orange squares.
(a)
(b)
(c)
(d)


### 4.2. Driving in the First Spectral Gap

We now consider a fixed driving frequency offb=3.7718kHz, 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 (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 [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 (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.Figure 5
Same as Figure3 but for a driving frequency of fb=3.7718 kHz, which lies in the first spectral gap. (a) Profile of solution for a driving amplitude of a=0.24816 μm. (b) The bifurcation diagram for fb=3.7718 kHz (shown here) is similar to the one presented in Figure 3(b), which was for a driving frequency in the second gap (fb=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 (see text), rather than to chaos, as in the case of the second gap breathers.
(a)
(b)
(c)Figure 6
Displacement profiles and corresponding Floquet multipliers of time-periodic solutions obtained at a driving frequency offb=3.7718 kHz for a driving amplitude of (a) a=50 nm, (b) a=0.94 μm, (c) a=1.2 μm, and (d) a=0.94 μm. These solutions correspond to the (A)–(D) labels of Figure 5(b). The lighter masses (i.e., the Steel beads) oscillate in-phase and are shown as blue diamonds whereas the heavier masses (the Tungsten Carbide ones) are shown as orange squares.
(a)
(b)
(c)
(d)Figure 7
(a) Displacement versus time over one period of motion of two adjacent steel massesu3(t) and u4(t) for a solution in the second spectral gap (with fb=6.3087 kHz and a=40 nm). In this case the motion is (nearly) out-of-phase. (b) Same as (a), but for a solution in the first spectral gap (with fb=3.7718 kHz and a=50 nm), in which case the motion is in-phase. (c) Velocities versus time for the motion of steel masses v2(t) and v15(t) for a solution with frequency fb=3.4429 kHz and amplitude a=50 nm. Solid lines correspond to numerical results while the dash-dotted lines correspond to experimental data (mean of the three runs). Notice that the tail oscillates with twice the frequency of the primary nodes near the left boundary, where 2fb lies in the second optical band.
(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 frequencyfb∈[3.3550,3.6850] will have a second harmonic that lies in the second optical band. Indeed, the breather solution depicted in Figure 7(c) has a frequency fb=3.4429, but the tail oscillates with twice that frequency (i.e., it is resonating with the linear mode at a frequency of 2fb). 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 Figures3 and 5. Note that the linear asymptotic equilibrium (shown by the dashed gray line and corresponding to |un-un+1|/δ0→0) 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 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, 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. Figure8 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, for example, |un-un+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 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 fb≈2.35 kHz and fb≈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 fb≈1–5.2 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].Figure 8
The bifurcation diagram corresponding to the maximum velocity of bead 4 as a function of the actuation frequencyfb and for a value of the actuation amplitude of a=0.24816 μm is presented with smooth curves. The color indexing is the same as in Figures 3 and 5. The black squares with error bars represent the experimentally measured mean values of the velocities which were obtained using three experimental runs. Furthermore, the light gray areas enclosed by vertical dash-dotted black lines correspond to the frequency gaps according to Table 2 (see also Figure 2). Note that panel (b) is a zoom-in of panel (a).
(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.


## 4.1. Driving in the Second Spectral Gap

We first consider a fixed driving frequency offb=6.3087kHz, 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 a≈22.75 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 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 Figure 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 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 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. See also the discussion in Section 4.3 for more details on possible explanations for disagreement in high amplitude regimes.Figure 3
Comparison between experimental and theoretical results for a driving frequencyfb=6.3087 kHz, which lies in the second spectral gap. (a) Profile of solution for a driving amplitude of a=32.16 nm. Experimental values (green shaded areas) are shown as a superposition of the velocities measured over the 50 ms time window of recorded data. The numerically predicted extrema (over one period) of the exact solution obtained via Newton’s method are shown with blue (open) circles. (b) The quantity v4 as a function of the amplitude a of the drive. For the numerically exact solutions, v4 is the maximum velocity of the 4th bead over one period of motion. Note that the blue segments correspond to stable parametric regions while the red and green ones correspond to real and oscillatory unstable parametric regions, respectively. Black squares with error bars represent the experimentally measured mean values of the quantity v4 over three experimental runs, where v4 was calculated as the maximum measured over the 50 ms time window of recorded data. Finally, the dashed gray line is the maximum velocity (over one period) of bead 4 as predicted by the linear solution of (17). (c) Zoomed out version of panel (b). The gray circles correspond to numerical results that were obtained by dynamically running the equations of motion and monitoring the maximum velocity v4 over the experimental time window. The experimental data in this panel are shown again by black squares.
(a)
(b)
(c)Figure 4
Displacement profiles and corresponding Floquet multipliers of time-periodic solutions obtained at a driving frequency offb=6.3087 kHz for a driving amplitude of (a) a=0.2 nm, (b) a=30 nm, (c) a=40 nm, and (d) a=0.1 μm. These solutions correspond to the (A)–(D) labels of Figure 3(b). The lighter masses (i.e., the Steel beads) oscillate (nearly) out-of-phase and are shown as blue diamonds whereas the heavier masses (the Tungsten Carbide ones) are shown as orange squares.
(a)
(b)
(c)
(d)


## 4.2. Driving in the First Spectral Gap

We now consider a fixed driving frequency offb=3.7718kHz, 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 (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 [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 (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.Figure 5
Same as Figure3 but for a driving frequency of fb=3.7718 kHz, which lies in the first spectral gap. (a) Profile of solution for a driving amplitude of a=0.24816 μm. (b) The bifurcation diagram for fb=3.7718 kHz (shown here) is similar to the one presented in Figure 3(b), which was for a driving frequency in the second gap (fb=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 (see text), rather than to chaos, as in the case of the second gap breathers.
(a)
(b)
(c)Figure 6
Displacement profiles and corresponding Floquet multipliers of time-periodic solutions obtained at a driving frequency offb=3.7718 kHz for a driving amplitude of (a) a=50 nm, (b) a=0.94 μm, (c) a=1.2 μm, and (d) a=0.94 μm. These solutions correspond to the (A)–(D) labels of Figure 5(b). The lighter masses (i.e., the Steel beads) oscillate in-phase and are shown as blue diamonds whereas the heavier masses (the Tungsten Carbide ones) are shown as orange squares.
(a)
(b)
(c)
(d)Figure 7
(a) Displacement versus time over one period of motion of two adjacent steel massesu3(t) and u4(t) for a solution in the second spectral gap (with fb=6.3087 kHz and a=40 nm). In this case the motion is (nearly) out-of-phase. (b) Same as (a), but for a solution in the first spectral gap (with fb=3.7718 kHz and a=50 nm), in which case the motion is in-phase. (c) Velocities versus time for the motion of steel masses v2(t) and v15(t) for a solution with frequency fb=3.4429 kHz and amplitude a=50 nm. Solid lines correspond to numerical results while the dash-dotted lines correspond to experimental data (mean of the three runs). Notice that the tail oscillates with twice the frequency of the primary nodes near the left boundary, where 2fb lies in the second optical band.
(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 frequencyfb∈[3.3550,3.6850] will have a second harmonic that lies in the second optical band. Indeed, the breather solution depicted in Figure 7(c) has a frequency fb=3.4429, but the tail oscillates with twice that frequency (i.e., it is resonating with the linear mode at a frequency of 2fb). 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 Figures3 and 5. Note that the linear asymptotic equilibrium (shown by the dashed gray line and corresponding to |un-un+1|/δ0→0) 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 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, 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. Figure8 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, for example, |un-un+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 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 fb≈2.35 kHz and fb≈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 fb≈1–5.2 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].Figure 8
The bifurcation diagram corresponding to the maximum velocity of bead 4 as a function of the actuation frequencyfb and for a value of the actuation amplitude of a=0.24816 μm is presented with smooth curves. The color indexing is the same as in Figures 3 and 5. The black squares with error bars represent the experimentally measured mean values of the velocities which were obtained using three experimental runs. Furthermore, the light gray areas enclosed by vertical dash-dotted black lines correspond to the frequency gaps according to Table 2 (see also Figure 2). Note that panel (b) is a zoom-in of panel (a).
(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.



---

*Source: 830978-2015-07-09.xml*