Discrete breathers in a mechanical metamaterial

Mechanical metamaterials are engineered structures [4, 12, 15, 54, 19, 59] whose macroscopic properties are primarily controlled by their geometry and may differ considerably from those of their building blocks [16, 30, 33, 39, 55]. In recent years, there has been a lot of interest in nonlinear dynamic response of flexible mechanical metamaterials, a new class of engineered materials that exploit large deformation and mechanical instabilities of their components to yield a desired collective response [4, 22]. Examples include metamaterials consisting of rotating rigid elements that are connected by flexible hinges [23, 25], multistable kirigami sheets [31], chains of bistable shells [52] and beams [43], as well as origami-inspired [56, 58] and linkage-based [60] deployable structures. These metamaterials can be designed to enable potential applications that include morphing surfaces, soft robotics, reconfigurable devices, mechanical logic and controlled energy absorption [11, 21, 27, 37, 40, 42, 34]. Recent studies have demonstrated that metamaterials of this type can be designed to control propagation of a variety of nonlinear waves [25, 22, 43, 41, 56, 57, 60].

In this work we consider the flexible mechanical metamaterial that was recently studied experimentally and theoretically in [25, 24, 22]. The experimentally realized system, schematically shown in Fig. 1, consists of a chain of pairs of cross-type rigid units made of LEGO bricks and connected by thin flexible polyester or plastic hinges [25, 24]. Under certain assumptions, the system can be described by a discrete model that assigns two degrees of freedom to each pair of rigid units: horizontal displacement and rotation. This system, in turn, can be approximated at the continuum level by a Klein-Gordon equation with cubic nonlinearity, a nonlinear wave-bearing system that possesses both soliton and cnoidal wave solutions [22]. In [25], the authors use a combination of experiments, direct numerical simulations of the discrete system and analysis of the continuum model to investigate traveling waves in this system that correspond to elastic vector solitons on the continuum level. They demonstrate that the metamaterial lattice may be designed to exhibit amplitude gaps where soliton propagation is forbidden, which, in turn, enables the design of soliton splitters and diodes. In [24] the anomalous nature of the soliton collisions in this system is explored. These developments clearly illustrate the promise of this type of nonlinear lattice in regards to the wave dynamics and interactions.

In this work we demonstrate that in certain parameter regimes the discrete system derived in [25] also exhibits a plethora of spatially localized, time-periodic patterns in the form of discrete breathers. These structures arise in terms of the angle and strain (relative displacement) variables. Similar to discrete breathers observed in other settings, including Josephson junction arrays [5, 51], forced-damped arrays of coupled pendula [17], electrical lattices [38, 26, 44], micromechanical systems [48, 47, 46] and granular chains [6, 14, 61, 8], they emerge as a result of the interplay of (discrete) dispersion and nonlinearity [1, 2, 28] and appear to be generic in the gaps of the linear excitation spectrum, as we will show below.

To construct such solutions for the metamaterial system, we start by analyzing the dispersion relation, which features optical and acoustic branches. We show that when the angle $\phi_{0}$ measuring the vertical offset between the neighboring horizontal hinges, takes values in certain parameter-dependent intervals, there is a frequency gap between the optical and acoustic branches that enables existence of discrete breathers. We then use the iterations of Newton’s method with a suitable initial guess and (once converged to a member of a solution family) parameter continuation to compute branches of discrete breather solutions that have frequency inside the gap and either bifurcate from or exist near the edges of the optical and acoustic bands. Stability of the obtained solutions is investigated using Floquet analysis.

As our first example, we consider the system parameters from [25] and show that in this case a branch of discrete breather solutions bifurcates from the edge of the optical band provided that the offset angle $\phi_{0}$ is above a certain threshold. Floquet analysis reveals that this branch eventually undergoes period-doubling bifurcations, and we compute the corresponding double-period solutions and investigate their stability.

As a second example, we consider a different set of parameters that enables existence of breathers for small offset angles $\phi_{0}$ in a certain interval. Choosing two different values in this interval, we compute branches of solutions that exist in the gap between the optical and acoustic bands. Here, our computation reveals a complex bifurcation diagram in the energy-frequency plane involving branches of symmetric and asymmetric discrete breather solutions and emergence of instability modes associated with real and complex Floquet multipliers. In particular, we find that the onset of real instability can take place via collisions of complex multipliers, as well as symmetry-breaking and period-doubling bifurcations. Another mechanism involves critical points of the breather’s energy as a function of its frequency (effectively, a saddle-center bifurcation), in line with the stability criterion established in [32] for discrete breathers in Fermi-Pasta-Ulam and Klein-Gordon lattices. We investigate the fate of some of the unstable solutions by perturbing them along the corresponding eigenmodes and show that in each case the ensuing dynamic evolution leads to a discrete breather that is effectively stable if one neglects the presence of small-magnitude complex eigenvalues. The computed primary branches have a snake-like form with multiple turning points, and the solution profiles often evolve in a nontrivial way along a branch, e.g., via the emergence of additional peaks or troughs in the strain and angle variables describing a discrete breather with even symmetry. Some features of the obtained bifurcation diagrams are reminiscent of the “snake-and-ladder” patterns observed in other nonlinear systems [3, 13, 50], although a detailed exploration of such a phenomenology is outside the scope of the present work.

The rest of the paper is organized as follows. In Sec. 2 we introduce the discrete model and formulate the problem. Analysis of the dispersion relation for the linearized system is presented in Sec. 3. In Sec. 4 we discuss a solution branch bifurcating from the edge of the optical mode for the parameter values in [25] and exhibiting period-doubling bifurcations. In Sec. 5 we consider another set of parameters and describe the complex bifurcation diagrams involving branches that exist in the gap between the optical and acoustic bands. Concluding remarks can be found in Sec. 6.

Motivated by experimental and theoretical investigations in [25], we consider a chain that consists of $2\times N$ cross-type rigid units of mass $m$ and moment of inertia $J$ connected by thin flexible hinges, as shown in Fig. 1. The neighboring horizontal hinges are shifted in the vertical direction by $a\tan\phi_{0}$, where $a$ is the center-to-center horizontal distance between the neighboring units (see the bottom panel of Fig. 1). The hinges are modeled as a combination of three linear springs, with stiffness parameters $k_{l}$, $k_{s}$ and $k_{\theta}$ corresponding to longitudinal stretching, shearing and bending, respectively. Following [25], we describe the dynamics of the system by two degrees of freedom for $n$-th vertical pair of rigid units: the longitudinal displacement $u_{n}(t)$ and the rotation angle $\theta_{n}(t)$ at time $t$. Here it is assumed [25] that the two rigid units in each vertical pair have the same displacement and rotate by the same amount but in the opposite directions, with positive direction of rotation defined as shown in the bottom panel of Fig. 1. Introducing dimensionless variables

$$\tilde{u}_{n}=\dfrac{u_{n}}{a},\qquad\tilde{t}=t\sqrt{\dfrac{k_{l}}{m}}$$

and parameters

$$\alpha=\dfrac{a}{2\cos\phi_{0}}\sqrt{\dfrac{m}{J}},\qquad K_{s}=\dfrac{k_{s}}{k_{l}},\qquad K_{\theta}=\dfrac{4k_{\theta}\cos^{2}\phi_{0}}{k_{l}a^{2}},$$

one obtains [25]

$$\begin{split}&\ddot{u}_{n}=u_{n+1}-2u_{n}+u_{n-1}-\frac{\cos(\theta_{n+1}+\phi_{0})-\cos(\theta_{n-1}+\phi_{0})}{2\cos(\phi_{0})}\\ &\dfrac{1}{\alpha^{2}}\ddot{\theta}_{n}=-K_{\theta}(\theta_{n+1}+4\theta_{n}+\theta_{n-1})+K_{s}\cos(\theta_{n}+\phi_{0})\bigg{(}\sin(\theta_{n+1}+\phi_{0})+\sin(\theta_{n-1}+\phi_{0})\\ &-2\sin(\theta_{n}+\phi_{0})\bigg{)}-\sin(\theta_{n}+\phi_{0})\bigg{(}2\cos(\phi_{0})(u_{n+1}-u_{n-1})+4\cos(\phi_{0})-\cos(\theta_{n+1}+\phi_{0})\\ &-2\cos(\theta_{n}+\phi_{0})-\cos(\theta_{n-1}+\phi_{0})\bigg{)},\end{split}$$

(1)

where we dropped the tildes in the rescaled displacement and time variables, and double dot denotes second time derivative. The total energy of the system is [25]

$$H=\sum_{n=1}^{N}\left[(\Delta_{n}^{l})^{2}+K_{s}(\Delta_{n}^{s})^{2}+\frac{K_{\theta}}{8\cos^{2}(\phi_{0})}(2(\delta_{n}^{h})^{2}+(\delta_{n}^{v})^{2})+\dot{u}_{n}^{2}+\frac{1}{4\alpha^{2}\cos^{2}(\phi_{0})}\dot{\theta}_{n}^{2}\right],$$

(2)

where

$$\begin{split}\delta_{n}^{h}&=\theta_{n+1}+\theta_{n},\qquad\delta_{n}^{v}=2\theta_{n},\\ \Delta_{n}^{l}&=u_{n+1}-u_{n}+\frac{1}{2\cos(\phi_{0})}\left[2\cos(\phi_{0})-\cos(\phi_{0}+\theta_{n})-\cos(\phi_{0}+\theta_{n+1})\right],\\ \Delta_{n}^{s}&=\frac{1}{2\cos(\phi_{0})}\left[\sin(\phi_{0}+\theta_{n+1})-\sin(\phi_{0}+\theta_{n})\right]\end{split}$$

characterize the deformation associated with horizontal ($\delta_{n}^{h}$, $\Delta_{n}^{l}$, $\Delta_{n}^{s}$) and vertical ($\delta_{n}^{v}$) hinges.

We consider discrete breather (DB) solutions of (1). These are time-periodic nonlinear waves with frequency $\omega$ and the corresponding period $T=2\pi/\omega$,

$$u_{n}(t+T)=u_{n}(t),\qquad\theta_{n}(t+T)=\theta_{n}(t)$$

(3)

that are spatially localized in terms of strain

$$w_{n}(t)=u_{n+1}(t)-u_{n}(t)$$

(4)

and angle $\theta_{n}(t)$ variables.

To obtain conditions for existence of DB solutions bifurcating from the linear modes, we need to study the linear spectrum of the problem first. To that effect, we linearize (1) around the undeformed configuration. This yields

$$\begin{split}&\ddot{u}_{n}=u_{n+1}-2u_{n}+u_{n-1}+\dfrac{1}{2}\tan\phi_{0}(\theta_{n+1}-\theta_{n-1})\\ &\dfrac{1}{\alpha^{2}}\ddot{\theta}_{n}=(K_{s}\cos^{2}\phi_{0}-\sin^{2}\phi_{0}-K_{\theta})(\theta_{n+1}+\theta_{n-1})\\ &-2(K_{s}\cos^{2}\phi_{0}+\sin^{2}\phi_{0}+2K_{\theta})\theta_{n}-\sin(2\phi_{0})(u_{n+1}-u_{n-1}).\end{split}$$

(5)

Considering plane-wave solutions $u_{n}(t)=Ue^{i(kn-\omega t)}$, $\theta_{j}(t)=\Theta e^{i(kn-\omega t)}$ of (5) in the limit of an infinite chain ($N\to\infty$), we obtain the following solvability condition:

$$\begin{split}&\left(\omega^{2}-4\sin^{2}\frac{k}{2}\right)\left[\frac{\omega^{2}}{\alpha^{2}}-2(K_{\theta}-K_{s}\cos^{2}\phi_{0}+\sin^{2}\phi_{0})\cos(k)-2(2K_{\theta}+K_{s}\cos^{2}\phi_{0}+\sin^{2}\phi_{0})\right]\\ &-2\tan(\phi_{0})\sin(2\phi_{0})\sin^{2}(k)=0,\end{split}$$

which yields explicit (but cumbersome) expressions for the acoustic, $\omega_{-}(k)$, and optical, $\omega_{+}(k)$, branches of the dispersion relation between the wave number $k$ and the frequency $\omega$. The two branches satisfy

$$\omega_{-}(0)=0,\qquad\omega_{+}(0)=\alpha\sqrt{2(3K_{\theta}+2\sin^{2}\phi_{0})}>0.$$

(6)

We now examine the evolution of the dispersion relation when the parameters $\alpha$, $K_{s}$ and $K_{\theta}$ are fixed, while $\phi_{0}$ is varied. Due to $2\pi$-periodicity and even symmetry about $k=\pi$, it suffices to consider wave numbers $k$ in $[0,\pi]$. In what follows, we consider two sets of parameters $\alpha$, $K_{s}$ and $K_{\theta}$. In the first representative example, we set $\alpha=1.8$, $K_{s}=0.02$, and $K_{\theta}=1.5\times 10^{-4}$ from [25]. In the second case we keep the same value of $K_{s}$ and set $\alpha=5$ and $K_{\theta}=0.01$. In both cases

$$\alpha^{2}(K_{\theta}+2K_{s}\cos^{2}\phi_{0})<2$$

(7)

is satisfied for all $\phi_{0}$, and we thus have

$$\omega_{-}(\pi)=\alpha\sqrt{2(K_{\theta}+2K_{s}\cos^{2}\phi_{0})}<\omega_{+}(\pi)=2.$$

(8)

Furthermore, one can show that for these parameter values the acoustic branch $\omega_{-}(k)$ has the maximum value at $k=\pi$ given in (8) for all $\phi_{0}$. Meanwhile, as shown in Fig. 2(a), the optical branch $\omega_{+}(k)$ has a maximum at $k=\pi$ and a minimum at $k=0$ for $0\leq\phi_{0}<\phi^{\prime}_{0}$, where

$$\phi^{\prime}_{0}=\arccos\left(\sqrt{\frac{1+\frac{K_{\theta}}{2}-\frac{1}{\alpha^{2}}}{1-K_{s}}}\right)$$

(9)

is obtained by setting $\omega^{\prime\prime}_{+}(\pi)=0$ at $\phi_{0}=\phi^{\prime}_{0}$ and using (7). The corresponding inflection point at $k=\pi$ is shown in Fig. 2(b). For $\phi_{0}^{\prime}<\phi_{0}<\phi_{0}^{\prime\prime}$, where

$$\phi^{\prime\prime}_{0}=\arcsin\left(\sqrt{\frac{1}{\alpha^{2}}-\frac{3}{2}K_{\theta}}\right),$$

(10)

$k=\pi$ becomes a local minimum, and $\omega_{+}(k)$ reaches its maximum at $k=k_{max}$ in $(0,\pi)$ and a global minimum at $k=0$ (see Fig. 2(c)). At $\phi_{0}=\phi_{0}^{\prime\prime}$, the case shown in Fig. 2(d), we have $\omega_{+}(0)=\omega_{+}(\pi)=2$, which together with the second expression in (6) yields (10). For $\phi_{0}>\phi^{\prime\prime}_{0}$, the optical branch has a global minimum $\omega_{+}(\pi)=2$ at $k=\pi$. As shown in Fig. 2(e), it has a local minimum at $k=0$ and the maximum at $k=k_{max}$ in $(0,\pi)$ until $\phi_{0}$ reaches the value

$$\begin{split}&\phi_{0}^{\prime\prime\prime}=\arccos\bigg{(}\frac{1}{2\alpha}\bigg{[}\Big{(}-2+\alpha^{2}(K_{s}(3K_{\theta}+2)+5K_{\theta}+4)\\ &-\sqrt{\alpha^{4}(K_{s}(3K_{\theta}+2)+K_{\theta})^{2}-4\alpha^{2}(K_{s}(3K_{\theta}-2)+5K_{\theta})+4}\Big{)}/(K_{s}+1)\bigg{]}^{1/2}\bigg{)},\end{split}$$

(11)

where $k=0$ becomes an inflection point ($\omega_{+}^{\prime\prime}(0)=0$); see Fig. 2(f). For $\phi_{0}>\phi_{0}^{\prime\prime\prime}$ the optical branch is inverted and has the maximum value at $k=0$ and the minimum value at $k=\pi$, as shown in Fig. 2(g).

We obtain $\phi_{0}^{{}^{\prime}}=0.5736$, $\phi_{0}^{{}^{\prime\prime}}=0.5888$ and $\phi_{0}^{{}^{\prime\prime\prime}}=0.6032$ for the parameters $\alpha=1.8$, $K_{s}=0.02$, $K_{\theta}=1.5\times 10^{-4}$. In the case $\alpha=5$, $K_{s}=0.02$, $K_{\theta}=0.01$, the evolution of the optical branch is similar to Fig. 2 but the critical values are $\phi_{0}^{{}^{\prime}}=0.1240$, $\phi_{0}^{{}^{\prime\prime}}=0.1588$ and $\phi_{0}^{{}^{\prime\prime\prime}}=0.1959$.

Let $k_{min}$ denote the wave number where the optical branch $\omega_{+}(k)$ reaches its minimum value. From the above discussion it follows that $k_{min}=0$ for $0\leq\phi_{0}\leq\phi_{0}^{\prime\prime}$, with the minimum value $\omega_{+}(0)=\alpha(6K_{\theta}+4\sin^{2}\phi_{0})^{1/2}$, and $k_{min}=\pi$ for $\phi_{0}>\phi_{0}^{\prime\prime}$, with the minimum value $\omega_{+}(\pi)=2$. Recalling that the acoustic branch has a maximum at $k=\pi$, we find that when

$$G=\omega_{+}(k_{min})-\omega_{-}(\pi)>0,$$

(12)

there is a band gap between the two branches. See Fig. 3 for examples of such a gap.

A DB solution with frequency $\omega$ inside the gap, i.e., $\omega_{-}(\pi)<\omega<\omega_{+}(k_{\min})$, may exist provided that

$$S=\omega_{+}(k_{min})-\frac{1}{2}\omega_{+}(k_{max})>0$$

(13)

holds in addition to (12) and $\omega>\omega_{+}(k_{max})/2$. Here $k_{max}$ is the wavenumber where the optical branch $\omega_{+}(k)$ reaches its maximum value. The fact that $\omega$ does not coincide with either optical or acoustic values for any wave number means that the breather is not in resonance with any linear modes, while the condition (13) eliminates the second harmonic resonances by ensuring that $2\omega>\omega_{+}(k)$ for all wave numbers. This enables both the spatial localization (due to its presence in the bandgap) and the non-resonance of the breather, as discussed, e.g., in [35].

Fig. 4 shows $G$ and $S$ defined in (12) and (13), respectively, as functions of $\phi_{0}$ for the first parameter set.

Both functions have a corner at $\phi_{0}=\phi_{0}^{\prime\prime}$ where $k_{min}$ changes from $0$ to $\pi$. Noting that $G$ changes sign from negative to positive for $\phi_{0}<\phi_{0}^{\prime\prime}$, when $k_{min}=0$, we set

$$G=\omega_{+}(0)-\omega_{-}(\pi)=\alpha\left(\sqrt{6K_{\theta}+4\sin^{2}(\phi_{0})}-\sqrt{2K_{\theta}+4K_{s}\cos^{2}(\phi_{0})}\right)=0$$

to find the critical angle

$$\phi_{0}^{*}=\arccos\sqrt{\frac{1+K_{\theta}}{1+K_{s}}}$$

(14)

above which (12) holds. The function $S$ in Fig. 4(b) also changes sign for $\phi_{0}<\phi_{0}^{\prime}$, where $k_{min}=0$ and $k_{max}=\pi$, so that

$$S=\omega_{+}(0)-\frac{1}{2}\omega_{+}(\pi)=\alpha\sqrt{6K_{\theta}+4\sin^{2}(\phi_{0})}-1=0$$

at

$$\phi_{0}^{**}=\arcsin\bigg{(}\sqrt{\frac{1}{4\alpha^{2}}-\frac{3}{2}K_{\theta}}\bigg{)},$$

(15)

and hence (13) holds for $\phi_{0}>\phi_{0}^{**}$. We find that $\phi_{0}^{*}=0.1400$ and $\phi_{0}^{**}=0.2811$ in this case. Thus for $\phi_{0}>0.2811$, both (12) and (13) hold, and DB solutions may exist with frequencies $\omega$ in the interval $(\omega_{+}(k_{\max})/2,\omega_{+}(k_{\min}))$; otherwise, first or second resonances set in. The example at $\phi_{0}=26\pi/180\approx 0.4538$, where (12) and (13) hold for $1<\omega<1.57906$, is shown in Fig. 3(a). As shown in Fig. 4(b), the frequency gap increases until $\phi_{0}^{\prime\prime}=0.5888$ and then starts decreasing. Note that for $\phi_{0}<\phi_{0}^{\prime\prime}$, DB solutions bifurcating from the optical band emerge from $k=0$ mode, while for $\phi_{0}$ above this threshold the breathers bifurcate from the $k=\pi$ mode.

The functions $G(\phi_{0})$ and $S(\phi_{0})$ for the second parameter set are shown in Fig. 5. Recall that in this case (9), (10) and (11) yield $\phi_{0}^{\prime}=0.1240$, $\phi_{0}^{\prime\prime}=0.1588$ and $\phi_{0}^{\prime\prime\prime}=0.1959$. One can see that (12) holds ($F(\phi_{0})>0$) for $\phi_{0}>\phi_{0}^{*}$, where $\phi_{0}^{*}=0.0992$ is found from (14). Meanwhile, $S(\phi_{0})$ is positive for $0\leq\phi_{0}<\phi_{0}^{***}$. To find this value, we observe that it is above $\phi_{0}^{\prime\prime\prime}$, which means that $k_{\min}=\pi$ and $k_{\max}=0$ in (13). Thus,

$$S=2-\dfrac{1}{2}\alpha\sqrt{6K_{\theta}+4\sin^{2}(\phi_{0})}=0$$

must hold at $\phi_{0}=\phi_{0}^{***}$, which yields

$$\phi_{0}^{***}=\arcsin\bigg{(}\sqrt{\frac{4}{\alpha^{2}}-\frac{3}{2}K_{\theta}}\bigg{)}.$$

(16)

We obtain $\phi_{0}^{***}=0.3906$ for the second parameter set. Thus, in this case (12) and (13) both hold when $0.0992<\phi_{0}<0.3906$. Examples of dispersion relations with band gaps for this parameter regime are shown in panels (b) and (c) of Fig. 3. Note that in both cases the maximum of the acoustic branch lies above the half of the maximum of the optical one, and hence the frequency range where DB solutions may exist includes the entire gap between the two bands. This is in contrast to the example shown in Fig. 3(a) for the first parameter set, where the breather frequency must exceed $\omega_{+}(\pi)/2=1$.

We first discuss DB solutions bifurcating from the optical $k=0$ mode for $\phi_{0}<\phi_{0}^{\prime\prime}$ for the parameters considered in [25] and associated with the experimental implementation of the metamaterial in that work: $\alpha=1.8$, $K_{s}=0.02$, $K_{\theta}=1.5\times 10^{-4}$. We set $\phi_{0}=26\pi/180\approx 0.4538$, which enables existence of DB solutions with frequency $\omega$ in $(1,1.57906)$. The corresponding dispersion relation is shown in Fig. 3(a).

To obtain the breathers with frequency $\omega$ and corresponding period $T=2\pi/\omega$, we consider a chain of $N=200$ elements and solve iteratively using Newton’s method the following equations:

$$\begin{pmatrix}\mathbf{u}(T)-\mathbf{u}(0)\\ \dot{\mathbf{u}}(T)-\dot{\mathbf{u}}(0)\\ \boldsymbol{\theta}(T)-\boldsymbol{\theta}(0)\\ \dot{\boldsymbol{\theta}}(T)-\dot{\boldsymbol{\theta}}(0)\end{pmatrix}=\mathbf{0},$$

where the vector functions $\mathbf{u}(t)$, $\dot{\mathbf{u}}(t)$, $\boldsymbol{\theta}(t)$, and $\dot{\boldsymbol{\theta}}(t)$ have the components $u_{n}(t)$, $\dot{u}_{n}(t)$, $\theta_{n}(t)$, and $\dot{\theta}_{n}(t)$, $n=1,\dots,N$, respectively. We perform numerical continuation in the frequency $\omega$, starting with $\omega=1.57$, just below the edge of the optical band at $k=0$. The initial guess is of the form

$$u_{n}=\varepsilon_{u}\tanh[\delta(n-N/2)],\quad\theta_{n}=\varepsilon_{\theta}\operatorname{sech}[\delta(n-N/2)],$$

(17)

where $\varepsilon_{u}$, $\varepsilon_{\theta}$ and $\delta$ are small. The dynamical evolution of Eq. (1) (over the prescribed period $T$) is performed using the symplectic fourth-order Runge-Kutta-Nyström algorithm [9] with free-end boundary conditions.

To study the linear stability of the obtained solutions, we use Floquet analysis. Setting $u_{n}(t)=\hat{u}_{n}(t)+\epsilon v_{n}(t)$ and $\theta_{n}(t)=\hat{\theta}_{n}(t)+\epsilon\gamma_{n}(t)$ in (1), where $\hat{u}_{n}(t)$ and $\hat{\theta}_{n}(t)$ comprise the DB solutions, and considering $O(\epsilon)$ terms, we obtain the linearized system

$$\begin{split}&\ddot{v}_{n}=v_{n+1}+v_{n-1}-2v_{n}-\left[\frac{-\sin(\hat{\theta}_{n+1}+\phi_{0})\gamma_{n+1}+\sin(\hat{\theta}_{n-1}+\phi_{0})\gamma_{n-1}}{2\cos(\phi_{0})}\right]\\ &\frac{1}{\alpha^{2}}\ddot{\gamma}_{n}=-K_{\theta}(\gamma_{n+1}+4\gamma_{n}+\gamma_{n-1})+K_{s}[\cos(\hat{\theta}_{n}+\phi_{0})\cos(\hat{\theta}_{n+1}+\phi_{0})\gamma_{n+1}\\ &-\sin(\hat{\theta}_{n}+\phi_{0})\sin(\hat{\theta}_{n+1}+\phi_{0})\gamma_{n}+\cos(\hat{\theta}_{n}+\phi_{0})\cos(\hat{\theta}_{n-1}+\phi_{0})\gamma_{n-1}\\ &-\sin(\hat{\theta}_{n}+\phi_{0})\sin(\hat{\theta}_{n-1}+\phi_{0})\gamma_{n}-2(\cos^{2}(\hat{\theta}_{n}+\phi_{0})-\sin^{2}(\hat{\theta}_{n}+\phi_{0}))\gamma_{n}]\\ &-[2\sin(\hat{\theta}_{n}+\phi_{0})\cos(\phi_{0})(v_{n+1}-v_{n-1})+2\cos(\hat{\theta}_{n}+\phi_{0})\cos(\phi_{0})(\hat{u}_{n+1}-\hat{u}_{n-1})\gamma_{n}\\ &+4\cos(\hat{\theta}_{n}+\phi_{0})\cos(\phi_{0})\gamma_{n}-(\cos(\hat{\theta}_{n}+\phi_{0})\cos(\hat{\theta}_{n+1}+\phi_{0})\gamma_{n}\\ &-\sin(\hat{\theta}_{n}+\phi_{0})\sin(\hat{\theta}_{n+1}+\phi_{0})\gamma_{n+1})-2(\cos^{2}(\hat{\theta}_{n}+\phi_{0})-\sin^{2}(\hat{\theta}_{n}+\phi_{0}))\gamma_{n}\\ &-(\cos(\hat{\theta}_{n}+\phi_{0})\cos(\hat{\theta}_{n-1}+\phi_{0})\gamma_{n}-\sin(\hat{\theta}_{n}+\phi_{0})\sin(\hat{\theta}_{n-1}+\phi_{0})\gamma_{n-1})],\end{split}$$

which is used to compute the monodromy matrix $\mathcal{F}$ defined by

$$\begin{pmatrix}\mathbf{v}(T)\\ \dot{\mathbf{v}}(T)\\ \boldsymbol{\gamma}(T)\\ \dot{\boldsymbol{\gamma}}(T)\end{pmatrix}=\mathcal{F}\begin{pmatrix}\mathbf{v}(0)\\ \dot{\mathbf{v}}(0)\\ \boldsymbol{\gamma}(0)\\ \dot{\boldsymbol{\gamma}}(0)\end{pmatrix},$$

where the vector functions $\mathbf{v}(t)$, $\dot{\mathbf{v}}(t)$, $\boldsymbol{\gamma}(t)$, and $\dot{\boldsymbol{\gamma}}(t)$ have the components $v_{n}(t)$, $\dot{v}_{n}(t)$, $\gamma_{n}(t)$, and $\dot{\gamma}_{n}(t)$, $n=1,\dots,N$, respectively. The Floquet multipliers $\mu$ are the eigenvalues of the matrix $\mathcal{F}$. The existence of a Floquet multiplier $\mu$ satisfying $|\mu|>1$ indicates the presence of instability. When the relevant instability-inducing multiplier is real, we refer to the instability as exponential, given the exponential nature of the associated growth. When such real multipliers arise, they come in pairs $(\mu,1/\mu)$ (one of which is outside, while the other is inside the unit circle). In the case of a complex multiplier quartet $(\mu,1/\mu,\bar{\mu},1/\bar{\mu})$ with $|\mu|>1$, the instability is referred to as oscillatory, given that oscillations accompany the exponential growth due to the imaginary part of the associated multipliers. The fact that the multipliers come in real pairs or complex quartets is a generic by-product of the Hamiltonian nature of the underlying lattice dynamical system.

Fig. 6(a) shows the energy $H$ of the breathers bifurcating from the $k=0$ mode as a function of the frequency $\omega$. As we will see below, using this pair of independent and dependent variables to illustrate our bifurcation diagrams allows us to connect the change in monotonicity of the energy-frequency curve with a potential stability change [32]. As illustrated in the insets, the amplitude of both the strain (4) and angle variables increases as the frequency is decreased away from the edge of the optical band, i.e., as the strength of the nonlinear contribution increases. The maximum modulus of the Floquet multipliers computed for this solution branch is shown by the blue curve in Fig. 6(c). One can see that it exceeds unity and rapidly increases near the end of the continuation. As illustrated in the inset (see also panels (a) and (b) of Fig. 7, which show the Floquet multipliers at the beginning and the end of the continuation, respectively), this is due to the emergence of a pair $(\mu,1/\mu)$ of real Floquet multipliers from $\mu=-1$ at $\omega=1.05155$. One of these has modulus greater than one and hence leads to an exponential instability, at the point $d$ in Fig. 6(b), which corresponds to a period-doubling bifurcation [28]. A second pair of real Floquet multipliers emerges from $\mu=-1$ at $\omega=1.05006$, leading to another exponential instability mode (not further explored). As the frequency is decreased, these multipliers first move away from the unit circle along the real line and then start moving back toward it, eventually colliding at $\mu=-1$ at $\omega=1.0499$, which corresponds to the point $e$ in Fig. 6(b) and is associated with another period-doubling bifurcation.

To compute the double-period solutions that arise as a result of the bifurcations at the points $d$ and $e$ along the single-period solution branch, we used the same iterative procedure as discussed above with the initial guess consisting of a single-period solution with twice the frequency perturbed along the corresponding unstable mode. Solutions along the bifurcating branches were then obtained using parameter continuation in frequency or energy. The resulting energy as a function of frequency for the double-period solutions (red and green curves) is shown in Fig. 6(b) for each case together with the single-period solution branch (blue curve) discussed above. The double-period solution curves are plotted at twice their actual frequency in order to facilitate the comparison with the single-period solution curve. Insets in Fig. 6(b) show examples of the symmetric breather solutions along the different double-period solution curves. As the insets of Fig. 6(c) reveal, the Floquet spectra of the double-period and single-period solution branches are markedly different. While the single-period solutions, as noted above, are characterized by an exponential period-doubling instability associated with a Floquet multiplier $\mu<-1$ for frequencies below the value at the bifurcation point $d$, the double-period branches exhibit an exponential instability associated with a Floquet multiplier satisfying $\mu>1$. As the bifurcation points are approached, the corresponding pairs of real multipliers collide at $\mu=-1$ for the parent single-period branch and at $\mu=1$ for the bifurcating branches.

To examine the nature of these bifurcations further, we plot in the top panel of Fig. 6(d) the largest modulus of real Floquet multipliers $\mu$ as a function of $\omega$ along the green and blue curves near the bifurcation point $d$. One can see that at the period-doubling bifurcation point $d$ the single-period branch develops an exponential instability associated with a Floquet multiplier $\mu<-1$ via a subcritical pitchfork bifurcation of the double-period branch, which has a pair of real multipliers $(\mu,1/\mu)$ with $\mu>1$. In the bottom panel of Fig. 6(d), we show the second largest modulus of the real Floquet multipliers near the bifurcation point $e$, where the second pair of real multipliers emerges near $\mu=-1$ for the single-period branch and near $\mu=1$ for the bifurcating red branch. Due to the presence of the first pair of real multipliers, all solutions are unstable near the bifurcation point $e$, as indicated in Fig. 6(c).

Note that the upper branch of the multivalued energy-frequency function corresponding to the unstable green double-period solution curve bifurcating from the point $d$ has a local minimum and a local maximum, marked by the dashed vertical lines in Fig. 6(e). As illustrated in the first four panels in Fig. 6(f), these extrema are associated with a change of multiplicity of the Floquet multiplier at $\mu=1$ along this branch and subsequent emergence or collision of a second pair of real Floquet multipliers. The change in multiplicity of the unit Floquet multiplier when $H^{\prime}(\omega)$ changes sign is consistent with the energy-based stability criterion proved in [32] for discrete breathers in Fermi-Pasta-Ulam and Klein-Gordon lattices. Note, however, that in this case the change in multiplicity does not lead to a stability change due to the presence of an additional pair of non-unit real multipliers at these frequency values. As we trace the solution curve toward the point $d$, this pair collides at $\mu=1$ on the unit circle at a bifurcation point and subsequently briefly remains on it (see panels 4 and 5 in Fig. 6(f)), while the solutions are still unstable due to the presence of complex multipliers $\mu$ satisfying $|\mu|>1$ (not shown in panel 5). However, as illustrated in panels 7 and 8 in Fig. 6(f), two pairs of real multipliers subsequently emerge on the real axis via collisions of complex conjugate pairs of multipliers. One of the pairs eventually collides on the unit circle at another bifurcation point, leaving a single pair (panel 9), which in turn collides at $\mu=1$ at the point $d$.

The enlarged view of the Floquet multiplier curve for the single-period solution branch and the insets shown in Fig. 7(c) reveal that the onset of the period-doubling instability is preceded by small-magnitude oscillatory instabilities associated with pairs of multipliers colliding on the unit circle and then moving slightly off it in the form of a quartet as discussed above. Note also that the Floquet multipliers $\mu$ form an arc on or near the unit circle. Using the linearization (5) of Eq. (1) about the uniform equilibrium state for an infinite chain, one can show [10] that the background state of the breather with period $T$ contributes the Floquet multipliers

$$\mu=e^{\pm i\omega_{\pm}(k)T},$$

(18)

where we recall from Sec. 3 that $\omega_{+}(k)$ and $\omega_{-}(k)$ are the optical and acoustic branches of the dispersion relation. As we vary $k$ from $0$ to $\pi$, we obtain arcs of multipliers along the unit circle. Such arcs corresponding to the upper optical ($\omega_{+}(k)$, red arc) and the bottom acoustic ($-\omega_{-}(k)$, light blue arc) bands are depicted in panels (a), (b) and (c) of Fig. 8 for different values of $\omega$ (and hence different $T=2\pi/\omega$ in (18)) along with the numerically computed Floquet multipliers (dark blue crosses) for the obtained DB solutions. There are also symmetric arcs (not shown in the figure) corresponding to the bottom optical ($-\omega_{+}(k)$) and the upper acoustic ($\omega_{-}(k)$) bands.

Under the mapping given by (18), the left ends of the arcs corresponding to the top optical and bottom acoustic bands, respectively, seen in Fig. 8, are associated with $\omega_{+}(\pi)$ and $-\omega_{-}(\pi)$. As $\omega$ is decreased, the two ends approach each other along the unit circle and eventually coincide when

$$e^{i2\pi\omega_{+}(\pi)/\omega}=e^{-i2\pi\omega_{-}(\pi)/\omega},$$

which yields

$$\frac{\omega_{+}(\pi)+\omega_{-}(\pi)}{\omega}=n,$$

where $n$ is a positive integer. We find that the first such collision takes place when $n=2$, which together with (8) yields

$$\omega=\frac{2+\alpha\sqrt{2(K_{\theta}+2K_{s}\cos^{2}\phi_{0})}}{2}\approx 1.2293.$$

This predicted value of $\omega=1.2293$ is close to the first significant peak shown in Fig. 7(c), although there are also two smaller peaks to the right of it at $\omega=1.231$ and $\omega=1.239$. This discrepancy between predicted and actual collision frequency values may be attributed to numerical accuracy of computing the Floquet multipliers, as well as possible effects of weak nonlinearity.

The solution curve shown in Fig. 6(a) was continued until the frequency $\omega=0.9972$, and thus includes solutions with frequencies $\omega\leq 1$. As noted in Sec. 3, these frequencies are associated with second harmonic resonances of the DB solution with the linear waves that have frequencies in the optical band. As a result, the corresponding solutions are no longer localized and instead possess non-decaying oscillatory wings. Such solutions are known as phantom breathers [36] or nanoptera [7, 49]. The latter term stems from their non-vanishing tails given the resonance with the linear modes. An example of a phantom breather with frequency $\omega=0.9972$ (red curve) is shown in Fig. 9 along with the regular (localized) DB solution at $\omega=1.02$ (dashed blue).

We now consider the Fourier spectrum associated with the dynamic evolution of the obtained breathers with prescribed frequency $\tilde{\omega}$. Fig. 10 shows the Fast Fourier Transform (FFT) results involving the dynamics simulated over a course of $100$ oscillation periods for two different values of $\tilde{\omega}$, along with the acoustic and optical bands shaded in gray. In the case $\tilde{\omega}=1.1$ (panel (a)), there are only two peaks at nonzero frequencies for the displayed range, at $\tilde{\omega}$ and $2\tilde{\omega}$, and the latter is clearly above the top of the optical band (the right shaded strip) at $\omega=2$. When $\tilde{\omega}=1.02$ (panel (b)), one can see a third nonzero-frequency peak in addition to $\tilde{\omega}$ and $2\tilde{\omega}$. This peak is at $\tilde{\omega}/2$ and is associated with the period-doubling instability, which is present at this frequency. Note that $2\tilde{\omega}$ is above the optical band (the right shaded strip), and $\tilde{\omega}/2$ is above the acoustic band (the left shaded strip), so there are no resonances with either optical or acoustic linear waves. In contrast, in the case $\tilde{\omega}=0.9972$ (not shown), the peak at $2\tilde{\omega}$ is just inside the optical band, and the second-harmonic resonance results in the phantom breather structure shown in Fig. 9.

As we have seen, the existence of DB solutions with frequencies inside the band gap requires rather large angles $\phi_{0}$ (above $16^{\circ}$) for the set of model parameters used in the previous subsection. Since large offset angles may render the present description of the system with only two degrees of freedom somewhat less accurate [29], we consider in what follows the parameters $\alpha=5$, $K_{s}=0.02$, $K_{\theta}=0.01$, which allow breather existence at smaller values of $\phi_{0}$.