remarkRemark Generalized breathers in time-periodic nonlinear latticesC.
Chong, D. E. Pelinovsky, G. Schneider

# On the Existence of Generalized Breathers and Transition Fronts in Time-
Periodic Nonlinear Lattices

Christopher Chong Department of Mathematics, Bowdoin College, Brunswick, ME
04011, USA
e-mail<EMAIL_ADDRESS>Dmitry E. Pelinovsky Department of Mathematics
and Statistics, McMaster University, Hamilton, Ontario, Canada, L8S 4K1
e-mail<EMAIL_ADDRESS>Guido Schneider Institut für Analysis, Dynamik
und Modellierung, Universität Stuttgart, 70569 Stuttgart, Germany
e-mail<EMAIL_ADDRESS>
###### Abstract

We prove the existence of a class of time-localized and space-periodic
breathers (called $q$-gap breathers) in nonlinear lattices with time-periodic
coefficients. These $q$-gap breathers are the counterparts to the classical
space-localized and time-periodic breathers found in space-periodic systems.
Using normal form transformations, we establish rigorously the existence of
such solutions with oscillating tails (in the time domain) that can be made
arbitrarily small, but finite. Due to the presence of the oscillating tails,
these solutions are coined generalized $q$-gap breathers. Using a multiple-
scale analysis, we also derive a tractable amplitude equation that describes
the dynamics of breathers in the limit of small amplitude. In the presence of
damping, we demonstrate the existence of transition fronts that connect the
trivial state to the time-periodic ones. The analytical results are
corroborated by systematic numerical simulations.

##  1 Introduction

The classical discrete breather is a fundamental coherent structure of
nonlinear lattices. They can be found in many fields, ranging from photonics,
electrical circuits, condensed matter physics, molecular biology, chemistry,
and phononics [12]. Breathers are relevant for applications, such as
information storage and transfer in the context of photonic crystals [4], but
are also rich mathematically and have inspired countless numerical and
analytical studies [23, 9]. The discrete breather is localized in space and
periodic in time with temporal frequency lying within a frequency gap [12].
Spatially periodic media can have frequency gaps, and hence, discrete
breathers are possible in such systems [21].

If breathers can be found in the frequency gap of spatially periodic media,
what can be found in the wavenumber gap of temporally periodic media? While
this question is a natural one to ask, it has only been very recently
addressed. In the context of a photonic time crystal, it was formally shown in
[34] that structures that are localized in time and periodic in space can be
found in the wavenumber bandgap of temporally periodic media. The structure
reported on had same defining features as the classic breather, but with the
role of space and time switched. Such solutions are called $q$-gap breathers,
where $q$ stands for the wavenumber.

In the presence of damping, so-called transition fronts are possible in a
$q$-gap, which connect the trivial state to time-periodic ones. $q$-gap
breathers and transition fronts were studied numerically and experimentally in
the context of a nonlinear phononic lattice in [7]. The experimental platform
therein was based on the one developed in [25], where bifurcations of time-
periodic solutions were studied.

It is the purpose of this paper to establish rigorously the existence of
$q$-gap breathers and transition fronts and to provide a tractable analytical
approximation of their dynamics. $q$-gap breathers are a new type of
structure, and are distinct from $q$-breathers, which are localized in
wavenumber and periodic in time [13]. Temporal localization can also be
achieved via other mechanisms, including zero-wavenumber gain modulation
instability [29] and nonlinear resonances [5, 46]. Integrable equations admit
such solutions explicitly, e.g., the Akhmediev breathers of the nonlinear
Schrödinger (NLS) equation [1] and its discrete counterpart, the Ablowitz-
Ladik lattice [2]. A feature that distinguishes $q$-gap breathers from other
temporally localized structures, like the ones just described, is the fact
that the underlying wavenumber lies in a $q$-gap.

Wavenumber bandgaps for the (possible) existence of $q$-gap breathers can be
found in a wide class of temporally periodic lattices. Indeed, there have been
many recent advances in experimental platforms for time-varying systems,
including photonic [41, 44, 45, 32], electric [35, 27, 36], and phononic
examples [38, 42, 31, 33, 24]. Controllable temporal localization has
potential applications in the creation of phononic frequency combs [16] (see
also [5, 46]), energy harvesting [28, 39], or acoustic signal processing [19].
The alternate mechanism for temporal localization that $q$-gap breathers
afford and the wide availability of platforms in which they may be implemented
suggest the potential utility of $q$-gap breathers in photonic, phononic,
electrical, and even chemical or biological applications [7].

###  1.1 Model equations and physical motivation

The mathematical model for the present study is a time-periodic nonlinear
lattice,

(1)
$\underline{m}\ddot{u}_{n}+c\dot{u}_{n}+k(t)u_{n}=F(u_{n+1}-u_{n})-F(u_{n}-u_{n-1})$

with mass $\underline{m}$, damping parameter $c\geq 0$, the time-periodic
modulation of the spring parameter $k(t)=k(t+T)$ for a period $T>0$, and the
inter-particle force $F$. Assuming Dirichlet boundary conditions
$u_{0}(t)=u_{N+1}(t)=0$ for some integer $N$, we have a $2N$-dimensional
dynamical system obtained from (1) at $n=1,2,\ldots,N$. We use
$U:=(u_{1},u_{2},\dots,u_{N})$ for further references in the main results.

We will consider a polynomial form of the inter-particle force

(2) $F(w)=K_{2}w-K_{3}w^{2}+K_{4}w^{3},\qquad K_{2}>0$

in which Eq. (1) corresponds to the classical Fermi-Pasta-Ulam-Tsingou (FPUT)
lattice if $k(t)=0$ and $c=0$ [11, 15]. The FPUT lattice is a central equation
in the study of nonlinear waves [43], partly due to its relevance as a model
in phononic, electrical, and biological systems (among others), its
mathematical richness [3], and its place in history as the first test-bed for
numerical simulations [8].

One concrete motivation for studying system (1) with a time-periodic stiffness
term $k(t)=k(t+T)$ is that it describes an array of repelling magnets
surrounded by time modulated coils. It was in this setting that $q$-gap
breathers were observed experimentally [7]. In this case $F(w)$ models the
repulsive force of the magnets, and is given by

(3) $F(w)=-\frac{a_{1}}{(d+w)^{a_{2}}},\quad$

where $d,a_{1},a_{2}>0$ are material parameters. Using the Taylor expansion of
Eq. (3) at $w=0$ gives a correspondence to the FPUT model with $F$ given by
Eq. (2) with

(4) $K_{2}=\frac{a_{2}a_{1}}{d^{a_{2}+1}},\quad
K_{3}=\frac{a_{2}(a_{2}+1)a_{1}}{2d^{a_{2}+2}},\quad
K_{4}=\frac{a_{2}(a_{2}+1)(a_{2}+2)a_{1}}{6d^{a_{2}+3}}.$

For $k(t)=k(t+T)$, we will use a specific choice for illustrations that is
motivated by the experimental set-up of [7]. In particular, we consider a
piecewise constant time-periodic parameter function $k(t)$ in the form:

(5) $k(t)=\left\\{\begin{array}[]{cl}k_{a},&t\in[0,\tau_{d}T),\\\
k_{b},&t\in[\tau_{d}T,T),\end{array}\right.$

for a $\tau_{d}\in[0,1]$ and where $k_{a},k_{b}$ are the so-called modulation
amplitude parameters and $\tau_{d}$ is the duty-cycle. Using the rescaling

$u_{n}(t)\rightarrow\frac{K_{2}}{K_{3}}u_{n}\left(\sqrt{\frac{K_{2}}{\underline{m}}}t\right)$

leads to the normalized parameter values with
$\underline{m},K_{2},K_{3}\rightarrow 1$. We note that the results of this
paper are applicable for more general parameter choices and time-periodic
coefficients $k(t)=k(t+T)$ and to lattices with $K_{j}(t)=K_{j}(t+T)$ for
$j=2,3,\ldots$.

(a) (b)

(c) (d)

Figure 1:  (a) The real (red) and imaginary (blue) part of the Floquet
exponent $\gamma$ as a function of Fourier wavenumber $q\in[0,\pi]$ in the
infinite lattice. The gray shaded region corresponds to the wavenumber
bandgap. The gray markers correspond to the Floquet exponents with a lattice
size of $N=10$. The $m_{0}=3$ exponents lie in wavenumber bandgap (larger
black markers). The parameter values are $\underline{m}=K_{2}=1$, $c=0$,
$T=1/0.37$, $\tau_{d}=0.5$, $k_{a}=0.6$ and $k_{b}=0.8$. (b) Floquet
multipliers corresponding to panel (a). One of the $m_{0}=3$ multipliers
(larger black markers) has modulus exceeding unity. (c) Conceptualization of a
generalized $q$-gap breather. The first component of $U^{+}_{\rm
hom}=(u_{1},u_{2},\dots,u_{N})$ is shown. (d) The analytical approximation of
the solution shown in panel (c), namely $\mathcal{U}^{+}$. The first component
of $\mathcal{U}^{+}$ is shown.

###  1.2 Summary of main results

We will develop rigorous proofs of the existence of oscillating homoclinic
solutions for $c=0$ and heteroclinic solutions for small damping $c>0$ of Eq.
(1) with time-periodic stiffness $k(t)$. Since the tails of homoclinic
solutions have small oscillations that do not vanish at infinity, the
solutions can also be thought of as generalized $q$-gap breathers. Similar
nomenclature has been adopted of the description of classical breathers with
non-zero tails in space-time continuous systems [17] and with spatially
periodic coefficients [10]. See [14] for discussion of how the interchange of
time and space variables affects derivation and justification of the
homoclinic solutions.

Before stating the main theorems of the paper, let us describe intuitively the
generalized $q$-gap breathers for $c=0$. In the presence of time-periodic
stiffness $k(t)$, some wavenumbers may fit into the gap in the dispersion
relationship, as seen in Fig. 1(a). The corresponding Floquet multipliers are
shown in panel (b) (details on the Floquet theory follow in Sec. 2). Unlike
the situation for frequency gaps in (linear) spatially periodic media,
exponential growth of Fourier modes occurs if the wavenumber is inside the gap
of the dispersion relation since the Floquet exponent in the gap has positive
real part, or equivalently, the Floquet multiplier has modulus exceeding
unity. Due to this (parametric) instability, initializing Eq. (1) with such a
Fourier mode will initially lead to growth. However, as the amplitude
increases, the nonlinearity of the system comes into play, and, as we shall
prove later, has a localizing affect on the dynamics, see Fig. 1(c). This
solution, however, cannot decay to zero. This is due to the presence of
neutrally stable modes (i.e., the Floquet multipliers lying on the unit
circle). During the dynamic evolution, all of the Fourier modes will couple
(due to the nonlinearity). The presence of the neutrally stable modes causes
the small oscillations, as seen in the tails of Fig. 1(c). From a dynamical
systems point of view, the trivial state has one unstable direction, one
stable direction, and $2N-2$ neutral directions. While genuine homoclinic
solutions could be found in the intersection of the associated one-dimensional
stable and unstable manifolds, it cannot be expected that such an intersection
exists in an $2N$-dimensional phase space. Thus, only homoclinic solutions
with small oscillating ripples for $|t|\to\infty$ exist. These solutions lie
in the intersection of the $2N-1$-dimensional center-stable manifold with the
$2N-1$-dimensional center-unstable manifold of the origin, for which we use
the time-reversibility of the system (1) with $k(t)$ given by (5). The
distance the wavenumber of the unstable Fourier mode is to the edge of the gap
defines a small parameter $\varepsilon>0$ (how exactly is detailed in Sec. 2).
With normal form transformations for time-periodic systems (details in Secs.
3-6) it can be shown that the oscillating ripples can be made arbitrarily
small, i.e., of order $\mathcal{O}(\varepsilon^{M-1})$ at the time scale of
$\mathcal{O}(\varepsilon^{-M+2})$ with some $M\in{\mathbb{N}}$ arbitrarily
large but fixed, see Fig. 1(c). Using a multiple-scale analysis (details in
Sec. 7), one can derive an explicit approximation of a genuine homoclinic
orbit, see Fig. 1(d), which agrees with the leading order of the generalized
breather’s profile.

We are now ready to present the main theorem on the existence of two
homoclinic orbits with oscillating ripples, i.e., generalized $q$-gap
breathers. The assumptions (Spec) and (Coeff) are described in Secs. 2 and 5
respectively. The small parameter $\varepsilon$ is defined in (Spec) and
explicitly in the representation for $k(t)=k(t+T)$.

######  Theorem 1.1.

Assume the spectral condition (Spec) and the normal form coefficient condition
(Coeff) are satisfied. Then for every $M\in\mathbb{N}$ with $M\geq 3$ there
exists an $\varepsilon_{0}>0$ and $C_{0}>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0})$ the system (1) with the time-periodic
coefficient (5) and $c=0$ possesses two generalized homoclinic solutions
$U_{\rm hom}^{\pm}\in
C^{1}([-\varepsilon^{-M+2},\varepsilon^{-M+2}],{\mathbb{R}}^{N})$ satisfying

$\sup_{t\in[-\varepsilon^{-M+2},\varepsilon^{-M+2}]}\|U^{\pm}_{\rm
hom}(t)-\mathcal{U}^{\pm}(t)\|+\|(U^{\pm}_{\rm
hom})^{\prime}(t)-(\mathcal{U}^{\pm})^{\prime}(t)\|\leq
C_{0}\varepsilon^{M-1}$

where $\mathcal{U}^{\pm}(t):\mathbb{R}\to\mathbb{R}^{N}$ satisfy
$\lim\limits_{|t|\to\infty}\|\mathcal{U}^{\pm}(t)\|+\|(\mathcal{U}^{\pm})^{\prime}(t)\|=0$
and can be approximated as

$(\mathcal{U}^{\pm})_{n}(t)=\pm\varepsilon A(\varepsilon
t)g(t)\sin(q_{m_{0}}n)+\mathcal{O}(\varepsilon^{2}),$

where $g(t+T)=-g(t)$ and $A(\tau)=\alpha{\rm sech}(\beta\tau)$ are uniquely
defined, real-valued functions with some $\alpha,\beta>0$, see Eq. (69) below.

In case of small damping $c=\mathcal{O}(\varepsilon)$ each of the two
homoclinic orbits $U_{\rm hom}^{\pm}$ of Theorem 1.1 breaks up. There exist
two non-zero anti-periodic solutions $\mathcal{U}_{\rm
per}^{\pm}(t+T)=-\mathcal{U}_{\rm per}^{\pm}(t)$ with a $2N$-dimensional
stable manifold for $c>0$.

Therefore, as can be seen by counting the dimensions, the one-dimensional
unstable manifold from the zero equilibrium intersects the $2N$-dimensional
stable manifolds of one the two non-zero anti-periodic solutions
$\mathcal{U}_{\rm per}^{\pm}$ transversally. In contrast to the oscillating
homoclinic orbits, the heteroclinic orbits have no oscillating ripples as
$t\to-\infty$ and converge to the orbits of the anti-periodic solutions
$\mathcal{U}_{\rm per}^{\pm}$ as $t\to+\infty$, see Figure 2.

Figure 2:  Example of the transition front $U^{+}_{\rm
het}=(u_{1},u_{2},\dots,u_{N})$, where the first component is shown. The first
component of the anti-periodic solution $\mathcal{U}^{+}_{\rm per}$ is also
shown as the light gray line.

Existence of the anti-periodic solutions $\mathcal{U}_{\rm per}^{\pm}$ is
guaranteed by the following theorem.

######  Theorem 1.2.

Assume the spectral condition (Spec) and the normal form coefficient condition
(Coeff). Fix $\widetilde{c}\geq 0$. Then there exists an $\varepsilon_{0}>0$
and $C_{0}>0$ such that for all $\varepsilon\in(0,\varepsilon_{0})$ the system
(1) with the time-periodic coefficient (5) and $c=\widetilde{c}\varepsilon>0$
possesses anti-periodic solutions $\mathcal{U}_{\rm per}^{\pm}$ such that
$\mathcal{U}_{\rm per}^{\pm}(t+T)=-\mathcal{U}_{\rm per}(t)$ and

$\sup_{t\in\mathbb{R}}\|\mathcal{U}^{\pm}_{\rm
per}(t)\|+\|(\mathcal{U}^{\pm}_{\rm per})^{\prime}(t)\|\leq C_{0}\varepsilon.$

The following theorem presents the main result on the existence of the
heteroclinic orbits (transition fronts) between the trivial solution $0$ and
the anti-periodic solutions $\mathcal{U}_{\rm per}^{\pm}$.

######  Theorem 1.3.

Assume the spectral condition (Spec) and the normal form coefficient condition
(Coeff). Fix $\widetilde{c}>0$. Then there exists an $\varepsilon_{0}>0$ such
that for all $\varepsilon\in(0,\varepsilon_{0})$ the system (1) with the time-
periodic coefficient (5) and $c=\widetilde{c}\varepsilon>0$ possesses two
heteroclinic solutions $U_{\rm het}^{\pm}\in
C^{1}({\mathbb{R}},{\mathbb{R}}^{N})$ such that

$\lim_{t\to-\infty}U^{\pm}_{\rm
het}(t)=0,\qquad\lim_{t\to-\infty}(U^{\pm}_{\rm het})^{\prime}(t)=0$

and

$\lim_{t\to+\infty}\inf_{t_{0}\in[0,T]}\|U^{\pm}_{\rm het}(t)-\mathcal{U}_{\rm
per}^{\pm}(t+t_{0})\|+\|(U^{\pm}_{\rm het})^{\prime}(t)-(\mathcal{U}_{\rm
per}^{\pm})^{\prime}(t+t_{0})\|=0,$

where $\mathcal{U}_{\rm per}^{\pm}$ are the anti-periodic solutions of (1)
from Theorem 1.2.

######  Remark 1.4.

Because of the quadratic nonlinearity in the FPUT system (1), the homoclinic,
anti-periodic, and heteroclinic orbits in Theorems 1.1-1.3 are not related by
the sign reflection even though the leading orders obtained from the cubic
normal form are related by the sign reflection, see Eq. (66). However, if
$K_{3}=0$, these orbits are related by the sign reflection up to any orders
due to the symmetry of the FPUT system (1).

The article is organized as follows. Section 2 presents the Floquet and
spectral analysis of the linearized FPUT system and introduce assumption
(Spec). Preparations for the normal form transformations are described in
Section 3. Normal form transformations are described in Section 4. The proof
of Theorem 1.1 is given in Section 5, where assumption (Coeff) is introduced.
Section 6 contains the proof of Theorems 1.2 and 1.3. A multiple-scale
analysis is carried out in Section 7, which provides tractable approximations
for both breathers and fronts. It also allows verification of (Coeff) through
direct computation. Numerical illustrations of the main results are described
in Section 8. Section 9 concludes the paper with a summary and brief
discussions.

Funding: This work was partially supported by the National Science Foundation
under grant number DMS-2107945 (C. Chong). D. E. Pelinovsky is partially
supported by AvHumboldt Foundation as Humboldt Research Award. G. Schneider is
partially supported by the Deutsche Forschungsgemeinschaft DFG through the
cluster of excellence ’SimTech’ under EXC 2075-390740016.

##  2 The linearized system

The linearized FPUT system at the trivial (zero) equilibrium is given by

(6)
$\underline{m}\ddot{u}_{n}+c\dot{u}_{n}+k(t)u_{n}=K_{2}(u_{n+1}-2u_{n}+u_{n-1})$

for $n=1,2,\ldots,N$ with Dirichlet boundary conditions
$u_{0}(t)=u_{N+1}(t)=0$. The linearized system (6) is solved by a linear
superposition of the discrete Fourier sine modes:

$u_{n}(t)=\sum_{m=1}^{N}\widehat{u}_{m}(t)\sin(q_{m}n),\quad q_{m}:=\frac{\pi
m}{N+1},\quad 1\leq m\leq N.$

The $m$-th Fourier mode has the amplitude $\widehat{u}_{m}(t)$ for which the
linear FPUT equation (6) transforms to the linear Schrödinger equation

(7) $\mathcal{L}\widehat{u}_{m}=K_{2}\omega^{2}(q_{m})\widehat{u}_{m},$

where

$\displaystyle\mathcal{L}$
$\displaystyle:=-\underline{m}\partial_{t}^{2}-c\partial_{t}-k(t),\quad
k(t+T)=k(t)$

and

$\displaystyle\omega^{2}(q):=4\sin^{2}\left(\frac{q}{2}\right),\quad
q\in[0,\pi].$

We will review solutions of the spectral problem (7) by using the Floquet
theory and the spectral theory of the Schrödinger operators.

###  2.1 Floquet theory

To obtain the monodromy matrix associated to (7) for general time-periodic
coefficients $k(t)$, one may resort to numerical computation or perturbation
analysis [26]. However, in the case of piecewise constant $k(t)$ as in (5),
this can be done explicitly [6, 37] (see also [7]). For the convenience of the
readers, we summarize the relevant results.

Let $\lambda_{m}:=K_{2}\omega^{2}(q_{m})$ for each $1\leq m\leq N$, and define
$s_{a},s_{b}>0$ by

(8)
$s_{a,b}:=\sqrt{\frac{\lambda_{m}+k_{a,b}}{\underline{m}}-\frac{c^{2}}{4\underline{m}^{2}}}.$

Note that $s_{a},s_{b}$ also depend on $m=1,2,\dots,N$ but the index $m$ is
dropped from the notation for simplicity. We obtain the exact solution of (7):

(9)
$\widehat{u}_{m}(t)=\left\\{\begin{array}[]{cl}e^{-\frac{ct}{2\underline{m}}}\left[A_{0}\cos(s_{a}t)+B_{0}\sin(s_{a}t)\right],&t\in[0,\tau_{d}T),\\\
e^{-\frac{ct}{2\underline{m}}}\left[C_{0}\cos(s_{b}(t-\tau_{d}T))+D_{0}\sin(s_{b}(t-\tau_{d}T))\right],&t\in[\tau_{d}T,T),\end{array}\right.$

with some constants $A_{0},B_{0},C_{0},D_{0}$. By $C^{1}$-continuity across
$t=\tau_{d}T$, we obtain

(10) $\left[\begin{array}[]{c}C_{0}\\\
D_{0}\end{array}\right]=\left[\begin{array}[]{cc}\cos(s_{a}\tau_{d}T)&\sin(s_{a}\tau_{d}T)\\\
-\frac{s_{a}}{s_{b}}\sin(s_{a}\tau_{d}T)&\frac{s_{a}}{s_{b}}\cos(s_{a}\tau_{d}T)\end{array}\right]\left[\begin{array}[]{c}A_{0}\\\
B_{0}\end{array}\right].$

The monodromy matrix $J$ is obtained as a mapping

(11) $\left\\{\begin{array}[]{l}\widehat{u}_{m}(0)=A_{0},\\\
\widehat{u}_{m}^{\prime}(0)=s_{a}B_{0}-\frac{c}{2\underline{m}}A_{0}\end{array}\right.\qquad\Rightarrow\qquad\left\\{\begin{array}[]{l}\widehat{u}_{m}(T)=e^{-\frac{cT}{2\underline{m}}}A_{1},\\\
\widehat{u}_{m}^{\prime}(T)=e^{-\frac{cT}{2\underline{m}}}\left[s_{a}B_{1}-\frac{c}{2\underline{m}}A_{1}\right],\end{array}\right.$

with

(20) $\displaystyle\left[\begin{array}[]{c}A_{1}\\\ B_{1}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{cc}\cos(s_{b}(1-\tau_{d})T)&\sin(s_{b}(1-\tau_{d})T)\\\
-\frac{s_{b}}{s_{a}}\sin(s_{b}(1-\tau_{d})T)&\frac{s_{b}}{s_{a}}\cos(s_{b}(1-\tau_{d})T)\end{array}\right]\left[\begin{array}[]{c}C_{0}\\\
D_{0}\end{array}\right]=:J\left[\begin{array}[]{c}A_{0}\\\
B_{0}\end{array}\right].$

Since $\det(J)=1$ and

(21) $\displaystyle{\rm trace}(J)$
$\displaystyle=2\cos(s_{a}\tau_{d}T)\cos(s_{b}(1-\tau_{d})T)-\frac{s_{a}^{2}+s_{b}^{2}}{s_{a}s_{b}}\sin(s_{a}\tau_{d}T)\sin(s_{b}(1-\tau_{d})T),$

the eigenvalues $\rho_{1}$ and $\rho_{2}$ of $J$ satisfies

$\rho_{1}\rho_{2}=1,\quad\rho_{1}+\rho_{2}={\rm trace}(J),$

with only three possibilities:

  * •

${\rm trace}(J)>2$ implies $0<\rho_{1}<1<\rho_{2}=\rho_{1}^{-1}$,

  * •

$-2\leq{\rm trace}(J)\leq 2$ implies
$\rho_{1}=\overline{\rho}_{2}\in\mathbb{C}$ with $|\rho_{1,2}|=1$,

  * •

${\rm trace}(J)<-2$ implies $\rho_{2}=\rho_{1}^{-1}<-1<\rho_{1}<0$.

The Floquet exponents of the mapping (11) are given by
$\gamma_{1,2}=\upsilon_{1,2}-\frac{c}{2\underline{m}}$, where

  * •

${\rm trace}(J)>2$ implies $\upsilon_{1,2}=\pm\log(\rho_{2})/T$,

  * •

$-2\leq{\rm trace}(J)\leq 2$ implies $\upsilon_{1,2}=\pm i\arg(\rho_{1})/T$,

  * •

${\rm trace}(J)<-2$ implies $\upsilon_{1,2}=i\pi/T\pm\log(|\rho_{2}|)/T$.

If $c=0$, the trivial solution $U=0$ is spectrally stable if all Floquet
exponents are purely imaginary. This corresponds to the case with $-2\leq{\rm
trace}(J)\leq 2$. Figure 3(a) shows the Floquet multipliers $\rho$ in the
critical case where ${\rm trace}(J)=-2$ with $m=3$. The corresponding Floquet
exponents $\gamma=\upsilon$ are shown in Fig. 3(b).

(a) (b) (c)

Figure 3: Bifurcation scenario for the parameter set $\underline{m}=K_{2}=1$,
$c=0$, $T=1/0.37$ and $\tau_{d}=0.5$. With the critical modulation amplitude
parameters $k_{a}^{0}=0.5$ and $k_{b}^{0}=0.79$, the critical Fourier mode is
$m_{0}=3$. (a) Plot of the Floquet multipliers, $\rho$, in the complex plane
(the unit circle is shown for visual aid). The $m_{0}=3$ multiplier lies
exactly at -1 on the unit circle (larger black marker). (b) The real (red) and
imaginary (blue) part of the Floquet exponent $\gamma$ as a function of
Fourier wavenumber $q\in[0,\pi]$ in the infinite lattice. The gray shaded
region corresponds to the wavenumber bandgap. The gray markers correspond to
the Floquet exponents with a lattice size of $N=10$. The $m_{0}=3$ exponent
lies exactly at the left edge of the wavenumber bandgap (larger black marker).
(c) Spectral bands (blue curves) of the Schrödinger operator as a function of
$\ell=\mathrm{imag}(\gamma)$, the imaginary part of the Floquet exponent. The
gray dots show the corresponding values in the finite lattice with $N=10$. The
eigenvalues $\mu_{1}$ and $\mu_{2}$ define the edges of the band gap, shown as
the gray shaded region. The $m_{0}=3$ mode lies exactly at the top of the
first spectral band (larger black marker).

###  2.2 Spectral theory

Let us review the spectral properties of the Schrödinger operator

$\mathcal{L}:H^{2}({\mathbb{R}})\subset L^{2}({\mathbb{R}})\to
L^{2}({\mathbb{R}})$

with a $T$-periodic coefficient $k(t)=k(t+T)$ in the particular case of $c=0$.
The spectrum of $\mathcal{L}$ is purely continuous and consists of bands
disjoint from each other by some gaps:

(22)
$\sigma(\mathcal{L})=[\nu_{0},\mu_{1}]\cup[\mu_{2},\nu_{1}]\cup[\nu_{2},\mu_{3}]\cup[\mu_{4},\nu_{3}]\cup\cdots$

where $\\{\nu_{j}\\}_{j=0}^{\infty}$ are eigenvalues of $\mathcal{L}f=\nu f$
with periodic boundary conditions $f(t+T)=f(t)$ and
$\\{\mu_{j}\\}_{j=1}^{\infty}$ are eigenvalues of $\mathcal{L}g=\mu g$ with
anti-periodic boundary conditions $g(t+T)=-g(t)$. Eigenvalues
$\\{\nu_{j}\\}_{j=0}^{\infty}$ correspond to ${\rm trace}(J)=2$ of the Floquet
theory, whereas eigenvalues $\\{\mu_{j}\\}_{j=0}^{\infty}$ correspond to ${\rm
trace}(J)=-2$ and the spectral bands correspond to $-2\leq{\rm trace}(J)\leq
2$.

######  Remark 2.1.

Note that $\lambda$ is parameterized by the wavenumber $q$ in
$\lambda=K_{2}\omega^{2}(q)$, which in turn determines the Floquet exponent
$\gamma$. Within the spectral bands, the corresponding Floquet exponents are
purely imaginary and are of the form $\gamma=i\ell$. An example plot showing
the dependence of $\lambda$ on $\ell$ is shown in Fig. 3(c). For this example
there is one band gap, which is the shaded region in the figure. The band gap
edges are given by $\mu_{1}$ (the bottom of the gap) and $\mu_{2}$ (the top of
the gap). This representation of the spectrum will be useful later when we
derive an amplitude equation for the description of the envelope of the
breather in Sec. 7. In particular, the concavity of $\lambda(\ell)$ will play
a key role.

We define the bifurcation for a particular stiffness $k(t)\equiv k_{0}(t)$,
for which $\mathcal{L}\equiv\mathcal{L}_{0}$, if there exists an integer
$1\leq m_{0}\leq N$ such that $\lambda_{m_{0}}=K_{2}\omega^{2}(q_{m_{0}})$
coincides with the end point of $\sigma(\mathcal{L}_{0})$ and
$\lambda_{m}=K_{2}\omega^{2}(q_{m})$ for $m\neq m_{0}$ are located inside
$\sigma(\mathcal{L}_{0})$. For the bifurcation shown in Figure 3,
$\lambda_{m_{0}}$ coincides with $\mu_{1}$, for which the two bifurcating
Floquet exponents are given by $\gamma_{0}=i\frac{\pi}{T}$ with
$\ell_{0}:=\frac{\pi}{T}$. This bifurcation corresponds to
$\lambda_{1}^{\prime\prime}(\ell_{0})<0$ and

(23) $[0,4K]\subset[\nu_{0},\nu_{1}],\quad
K_{2}\omega^{2}(q_{m_{0}})=\mu_{1}.$

Equivalently, we can obtain a bifurcation when $\lambda_{m_{0}}$ coincides
with $\mu_{2}$ with $\lambda_{2}^{\prime\prime}(\ell_{0})>0$ and

(24) $[0,4K]\subset[\nu_{0},\nu_{1}],\quad
K_{2}\omega^{2}(q_{m_{0}})=\mu_{2}.$

The two bifurcating Floquet exponents $\gamma_{0}$ correspond to the Floquet
multiplier $\rho$ at $-1$. All other Floquet exponents $\gamma$ are assumed to
be on the imaginary axis bounded away from $0$ and $\gamma_{0}$. The
corresponding Floquet multipliers $\rho$ are on the unit circle bounded away
from $+1$ and $-1$.

The bifurcation in terms of Floquet multipliers is shown in Fig. 3(a), whereas
Fig. 3(b) shows the bifurcation in terms of Floquet exponents. The spectral
bands of the Schrödinger operator $\mathcal{L}_{0}$ are shown in Fig. 3(c).
Notice that the discrete mode $m_{0}=3$ lies exactly at
$(\ell,\lambda)=(\ell_{0},\mu_{1})$, where $\ell_{0}=\frac{\pi}{T}$.

###  2.3 Spectral assumption and defining the small parameter $\varepsilon$

With the linear theory in hand, we can now specify the spectral assumption as
follows.

(Spec) There exists a periodic coefficient $k_{0}(t+T)=k_{0}(t)$, for which
all Floquet exponents lie on the imaginary axis. With the exception of two
exponents at $\frac{i\pi}{T}$, they are assumed to be simple and non-zero. For
small $\varepsilon>0$ we assume that the two Floquet exponents at
$\frac{i\pi}{T}$ split symmetrically from the imaginary axis along the real
axis to the order of $\mathcal{O}(\varepsilon)$.

We can define $\varepsilon$ more explicitly by using the decomposition

(25) $k(t)=k_{0}(t)+\delta\varepsilon^{2},$

where $\delta$ is a proper sign factor. For $c=0$, the small parameter
$\varepsilon$ is related to the distance of the critical Floquet exponents
$\gamma=\upsilon$ from the imaginary axis in the following way. The real part
of the Floquet exponent $\gamma$ which depends on $\varepsilon$ is given by

$\mathrm{Re}(\gamma)=\frac{1}{T}\mathrm{cosh}^{-1}\left(-\frac{1}{2}\mathrm{trace}(J)\right).$

Since we know from (Spec) that $\mathrm{trace}(J)=-2$ with $\varepsilon=0$
(for $k(t)=k_{0}(t)$), a series expansion of the real part of the Floquet
exponent about $\varepsilon=0$ yields
$\mathrm{Re}(\gamma)=\mathcal{O}(\varepsilon)$, where we used the fact that
$\mathrm{cosh}^{-1}(1+w)\approx\sqrt{2w}$ and
$\mathrm{trace}(J)=-2+\mathcal{O}(\varepsilon^{2})$. In Sec. 7.2 we will show
that

(26)
$\mathrm{Re}(\gamma)=\varepsilon\frac{\sqrt{2}}{\sqrt{|\lambda^{\prime\prime}(\ell_{0})|}}+\mathcal{O}(\varepsilon^{2}),$

where $\lambda(\ell)$ is the corresponding band of $\mathcal{L}_{0}$ at
$\ell_{0}=\frac{\pi}{T}$ and the sign of $\delta$ is selected to be the
opposite of the sign of $\lambda^{\prime\prime}(\ell_{0})$.

######  Remark 2.2.

One can relate the small parameter $\varepsilon$ to the distance of the
bifurcating wavenumber $q_{m_{0}}$ to the band edge in the following way. For
a fixed $\varepsilon$, suppose the wavenumber bandgap is $[q_{\ell},q_{r}]$,
where the left edge $q_{\ell}$ and right edge $q_{r}$ depend on $\varepsilon$
and can be found by solving trace($J$)$=-2$ with $\lambda=K_{2}\omega^{2}(q)$.
Suppose that the critical wavenumber $q_{m_{0}}$ coincides with the left band
edge at the bifurcation point, (i.e., $q_{m_{0}}=q_{\ell}$ when
$\varepsilon=0$). Then, for $\varepsilon>0$, the distance to the band edge is
$\Delta q=q_{\ell}-q_{m_{0}}$. By inspection of (8) and (21), if one knows the
critical values of $k_{a}^{0}$ and $k_{b}^{0}$ then $\Delta q$ can be
determined by solving $\lambda(\ell(q_{m_{0}}))=\lambda(\ell(q_{m_{0}}+\Delta
q))+\delta\varepsilon^{2}$ which yields

(27) $\Delta
q=q_{m_{0}}-2\sin^{-1}\left(\sqrt{\sin^{2}(q_{m_{0}}/2)-\frac{\delta\varepsilon^{2}}{4K_{2}}}\right)=\frac{\delta\varepsilon^{2}}{\partial_{q}\lambda(\ell(q_{m_{0}}))}+\mathcal{O}(\varepsilon^{4}).$

Thus $\Delta q=\mathcal{O}(\varepsilon^{2})$.

##  3 Normal form transformations

The FPUT system (1) consists of $N$ oscillators with Dirichlet boundary
conditions. By augmenting the vector $U(t)\in\mathbb{R}^{N}$ with
$U^{\prime}(t)\in\mathbb{R}^{N}$ as the vector $V(t)\in\mathbb{R}^{2N}$, we
rewrite the $2N$-dimensional time-periodic system in the abstract form:

(28) $\dot{V}(t)=Q(t)V(t)+N(V(t)),$

with the time-periodic coefficient matrix
$Q(t)=Q(t+T)\in{\mathbb{R}}^{2N\times 2N}$ being piecewise continuous on
$[0,T]$ for a period $T>0$ and the nonlinear function
$N(V):\mathbb{R}^{2N}\to\mathbb{R}^{2N}$ being smooth at $V=0$ with $N(0)=0$
and $D_{V}N(0)=0$. The solutions of the linear system

(29) $\dot{V}(t)=Q(t)V(t)$

are, according to Floquet’s theorem, of the form

(30) $V(t)=P(t)e^{\Lambda t}V(0)$

with a $T$-time-periodic matrix function $P(t)=P(t+T)\in{\mathbb{R}}^{2N\times
2N}$ and a time-independent matrix $\Lambda\in{\mathbb{R}}^{2N\times 2N}$,
eigenvalues of which coincide with Floquet exponents in Section 2.

######  Remark 3.1.

Eigenvalues $\gamma$ of the matrix $\Lambda$ are uniquely defined in the
strip:

(31) $-\frac{\pi}{T}<{\rm Imag}(\gamma)\leq\frac{\pi}{T}.$

Eigenvalues of $\Lambda$ are generally complex-valued but we use the
presentation (30) with real $P(t)$ and real $\Lambda$ for convenience of the
normal form transformations. For example, if $\gamma=\alpha\pm i\beta$ are two
complex-conjugate eigenvalues of $\Lambda$, then the canonical form for the
corresponding block of $\Lambda\in{\mathbb{R}}^{2N\times 2N}$ is

$\left[\begin{array}[]{cc}\alpha&\beta\\\ -\beta&\alpha\end{array}\right].$

######  Remark 3.2.

As preparations for the normal form transformation, we can consider the time-
periodic system (28) on the double period $2T$. The advantage of this approach
is that the bifurcating Floquet exponents in the stripe (31) correspond to
zero Floquet exponents in the $2T$-periodic system. The solution (30) of the
linear system (29) can be rewritten in the form:

(32) $V(t)=\widetilde{P}(t)e^{\widetilde{\Lambda}t}V(0),$

with the time-periodic matrix function
$\widetilde{P}(t)=\widetilde{P}(t+2T)\in{\mathbb{R}}^{2N\times 2N}$ and the
time-independent matrix $\widetilde{\Lambda}\in{\mathbb{R}}^{2N\times 2N}$.
According to the assumption (Spec) at the bifurcation point,
$\widetilde{\Lambda}$ has a double zero eigenvalue and all other (purely
imaginary) eigenvalues of $\widetilde{\Lambda}$ are simple and bounded away
from $0$.

######  Remark 3.3.

For the normal form transformation in Section 4 we need the following property
of $\widetilde{P}(t)$. Comparing the two representations of the fundamental
matrix solution

$\Phi(t)=P(t)e^{\Lambda t}=\widetilde{P}(t)e^{\widetilde{\Lambda}t},$

we obtain

$\widetilde{P}(t)=P(t)e^{\pi it/T}=\sum_{m\in{\mathbb{Z}}}P_{m}e^{2\pi
imt/T}e^{\pi it/T},$

with $P_{m}$ being constant $2N\times 2N$-matrices.

We now transform the system (28) on the double period to a convenient form for
which the linear part is autonomous in $t$. Let $V(t)=\widetilde{P}(t)W(t)$,
then $W(t)\in\mathbb{R}^{2N}$ satisfies the time-periodic system:

(33)
$\dot{W}(t)=\widetilde{\Lambda}W(t)+\widetilde{P}(t)^{-1}N(\widetilde{P}(t)W(t)).$

We define the projection $\Pi_{0}$ on the subspace associated with the double
zero eigenvalue of $\widetilde{\Lambda}$ by

(34) $\Pi_{0}=\frac{1}{2\pi i}\int_{\Gamma_{0}}(\lambda
I-\widetilde{\Lambda})^{-1}d\lambda,$

where $\Gamma_{0}$ is a closed curve surrounding the origin in the $\lambda$
plane counter-clockwise. The projection on the two other $(2N-2)$ eigenvalues
of $\widetilde{\Lambda}$ on the imaginary axis is defined by
$\Pi_{h}=I-\Pi_{0}$. The range of $\Pi_{0}$ is two-dimensional and the range
of $\Pi_{h}$ is $(2N-2)$-dimensional.

We apply these projections on system (33) and find for $W_{0}=\Pi_{0}W$ and
$W_{h}=\Pi_{h}W$ that

(35) $\displaystyle\dot{W}_{0}(t)$ $\displaystyle=$
$\displaystyle\Lambda_{0}W_{0}(t)+N_{0}(W_{0},W_{h}),$ (36)
$\displaystyle\dot{W}_{h}(t)$ $\displaystyle=$
$\displaystyle\Lambda_{h}W_{h}(t)+N_{h}(W_{0},W_{h})+H(W_{0}),$

where we have introduced $\Pi_{0}\widetilde{\Lambda}=\Lambda_{0}\Pi_{0}$,
$\Pi_{h}\widetilde{\Lambda}=\Lambda_{h}\Pi_{h}$,

$\displaystyle N_{0}(W_{0},W_{h})$
$\displaystyle:=\Pi_{0}\widetilde{P}(t)^{-1}N(\widetilde{P}(t)W(t)),$
$\displaystyle N_{h}(W_{0},W_{h})+H(W_{0})$
$\displaystyle:=\Pi_{h}\widetilde{P}(t)^{-1}N(\widetilde{P}(t)W(t)).$

The splitting into $N_{h}(W_{0},W_{h})+H(W_{0})$ is justified with
$N_{h}(W_{0},0)=0$. System (35)-(36) is extended by the additional equation
$\dot{\varepsilon}=0$, where $\varepsilon$ is the bifurcation parameter in
(Spec).

######  Remark 3.4.

In the context of the time-periodic system (33) on the double period, we
recall that $W_{0}$ represents the modes associated to the two Floquet
exponents which split from the double zero and leave the imaginary axis and
that $W_{h}$ represents the modes associated to the other $(2N-2)$ Floquet
exponents which stay on the imaginary axis for small bifurcation parameter
$\varepsilon$.

We use the normal form transformations to reduce the order of $H(W_{0})$ in
terms of powers of $\|W_{0}\|$.

######  Lemma 3.5.

For every $M\geq 2$ there exists an $\varepsilon_{0}>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0})$ there exists a change of coordinates

$W_{0,M}=W_{0},\qquad W_{h,M}=W_{h}+G(W_{0})$

such that system (35)-(36) transforms into

(37) $\displaystyle\dot{W}_{0}$ $\displaystyle=$
$\displaystyle\Lambda_{0}W_{0}+N_{0,M}(W_{0},W_{h,M}),$ (38)
$\displaystyle\dot{W}_{h,M}$ $\displaystyle=$
$\displaystyle\Lambda_{h}W_{h,M}+N_{h,M}(W_{0},W_{h,M})+H_{M}(W_{0}),$

with $N_{h,M}(W_{0},0)=0$ and $H_{M}(W_{0})=\mathcal{O}(\|W_{0}\|^{M})$.

######  Proof 3.6.

We set $W_{h,2}=W_{h}$ and then inductively

$W_{h,n+1}=W_{h,n}+G_{n}(W_{0}),$

with $G_{n}(W_{0})$ being a $n$-linear mapping in $W_{0}$. After the
transformations we have a system of the form

(39) $\displaystyle\dot{W}_{0}$ $\displaystyle=$
$\displaystyle\Lambda_{0}W_{0}+N_{0,n+1}(W_{0},W_{h,n+1}),$ (40)
$\displaystyle\dot{W}_{h,n+1}$ $\displaystyle=$
$\displaystyle\Lambda_{h}W_{h,n+1}+N_{h,n+1}(W_{0},W_{h,n+1})+H_{n+1}(W_{0}),$

with

$\displaystyle N_{0,n+1}(W_{0},W_{h,n+1})$ $\displaystyle=$ $\displaystyle
N_{0,n}(W_{0},W_{h,n}-G_{n}(W_{0})),$ $\displaystyle
N_{h,n+1}(W_{0},W_{h,n+1})$ $\displaystyle=$ $\displaystyle
N_{h,n}(W_{0},W_{h,n}-G_{n}(W_{0}))$ $\displaystyle $
$\displaystyle-(D_{W_{0}}G_{n}(W_{0}))(N_{0,n}(W_{0},W_{h,n})-N_{0,n}(W_{0},0)),$
$\displaystyle H_{n+1}(W_{0})$ $\displaystyle=$ $\displaystyle
H_{n}(W_{0})+\Lambda_{h}G_{n}(W_{0})$ $\displaystyle $
$\displaystyle-(D_{W_{0}}G_{n}(W_{0}))(\Lambda_{0}W_{0}+N_{0,n}(W_{0},0)),$

and so $N_{h,n+1}(W_{0},0)=0$. In order to have
$H_{n+1}(W_{0})=\mathcal{O}(\|W_{0}\|^{n+1})$ if
$H_{n}(W_{0})=\mathcal{O}(\|W_{0}\|^{n})$ we have to choose $G_{n}$ such that

(41)
$H_{n,n}(W_{0})+\Lambda_{h}G_{n}(W_{0})-(D_{W_{0}}G_{n}(W_{0}))\Lambda_{0}W_{0}=0,$

where $H_{n,n}(W_{0})$ is the $n$-linear part of $H_{n}(W_{0})$, i.e.
$H_{n}(W_{0})-H_{n,n}(W_{0})=\mathcal{O}(\|W_{0}\|^{n+1})$. Since

$H_{n,n}(W_{0})(t)=H_{n,n}(W_{0})(t+2T)$

we also have to choose $G_{n}(W_{0})(t)=G_{n}(W_{0})(t+2T)$. Using Fourier
series

$G_{n}(W_{0})(t)=\sum_{m\in{\mathbb{Z}}}G_{n}(W_{0})[m]e^{\pi imt/T}$

it is easy to see that (41) can be solved w.r.t. $G_{n}$ since according to
the assumption (Spec), none of the eigenvalues $\lambda_{j}$ of $\Lambda_{h}$
are located at $0$ and $\Lambda_{0}$ has a double zero eigenvalue. As a
result, the term $H(W_{0})$ can be made arbitrarily small in terms of powers
of $\|W_{0}\|$.

######  Remark 3.7.

The sequence of normal form transformations is not convergent and so we stop
after $M-1$ transformations with some large but fixed $M\in\mathbb{N}$, for
which we have $H_{M}(W_{0})=\mathcal{O}(\|W_{0}\|^{M})$. Note that the minimum
of $H_{M}(W_{0})$ is attained for $M=\mathcal{O}(1/\|W_{0}\|)$-many
transformations, after which $H_{M}$ is exponentially small in terms of
$\|W_{0}\|$, cf. [30].

##  4 Normal form transformations for the reduced system

If we ignore the terms $H_{M}(W_{0})$ in (37)-(38), then $\\{W_{h,M}=0\\} $ is
an invariant subspace for (37)-(38). In this two-dimensional subspace the
reduced system is obtained by setting $W_{h,M}=0$ in the first equation. So we
consider the two-dimensional ODE

(42) $\dot{W}_{0}=\Lambda_{0}W_{0}+N_{0,M}(W_{0},0).$

At the bifurcation point $\Lambda_{0}$ posseses a double eigenvalue
$\lambda=0$ with algebraic multiplicity two and geometric multiplicity one.
Thus, we have a Jordan-block of size two. The eigenvector of
$\widetilde{\Lambda}$ is denoted with $\varphi_{1}$ and the generalized
eigenvector with $\varphi_{2}$, i.e. $\widetilde{\Lambda}\varphi_{1}=0$ and
$\widetilde{\Lambda}\varphi_{2}=\varphi_{1}$. If we introduce coordinates $A$,
$B$ by

(43) $W_{0}=A\varphi_{1}+B\varphi_{2}$

we can rewrite (42) as the following two-dimensional system

(44a) $\displaystyle\dot{A}$ $\displaystyle=$ $\displaystyle B+f_{A}(t,A,B),$
(44b) $\displaystyle\dot{B}$ $\displaystyle=$
$\displaystyle\widetilde{\varepsilon}^{2}A+f_{B}(t,A,B),$

where $f_{A}$ and $f_{B}$ stand for real-valued nonlinear terms which are of
the form

$f_{A}(t,A,B)=\sum_{n=2}^{\infty}\sum_{j=0}^{n}\sum_{m\in{\mathbb{Z}}}f_{A,n,j,m}A^{j}B^{n-j}e^{2im\pi
t/T}e^{i(n-1)\pi t/T}$

and similarly for $f_{B}$, where $f_{A,n,j,m}$ and $f_{B,n,j,m}$ are
independent of time. To simplify notations, we introduce here the normalized
small parameter $\widetilde{\varepsilon}$ for the distance of the two Floquet
exponents from the imaginary axis. According to the expansion (26), we have

(45)
$\widetilde{\varepsilon}=\varepsilon\frac{\sqrt{2}}{\sqrt{|\lambda^{\prime\prime}(\ell_{0})|}}+\mathcal{O}(\varepsilon^{2}).$

For analyzing (44) we simplify this system by eliminating a number of terms by
another normal form transformation

######  Lemma 4.1.

There exists an $\tilde{\varepsilon}_{0}>0$ such that for all
$\tilde{\varepsilon}\in(0,\tilde{\varepsilon}_{0})$ there exists a change of
coordinates

$A=A_{3}+F_{A}(A_{3},B_{3}),\qquad B=B_{3}+F_{B}(A_{3},B_{3}),$

with $F_{A}$ and $F_{B}$ polynomials not containing linear terms, such that
system (44a)-(44b) transforms into

$\displaystyle\dot{A}_{3}$
$\displaystyle=B_{3}+f_{A,3,3,-1,2}A_{3}^{3}+f_{A,3,2,-1,2}A_{3}^{2}B_{3}+f_{A,3,1,-1,2}A_{3}B_{3}^{2}+f_{A,3,0,-1,2}B_{3}^{3}+\mathcal{O}(|A_{3}|^{4}+|B_{3}|^{4}),$
$\displaystyle\dot{B}_{3}$
$\displaystyle=\widetilde{\varepsilon}^{2}A_{3}+f_{B,3,3,-1,2}A_{3}^{3}+f_{B,3,2,-1,2}A_{3}^{2}B_{3}+f_{B,3,1,-1,2}A_{3}B_{3}^{2}+f_{B,3,0,-1,2}B_{3}^{3}+\mathcal{O}(|A_{3}|^{4}+|B_{3}|^{4}),$

with real-valued coefficients $f_{A,n,j,m,2}$ and $f_{B,n,j,m,2}$.

######  Proof 4.2.

It is well known that all terms which have a pre-factor which is oscillating
in time can be eliminated by a normal form transform or equivalently by
averaging, cf. [18, 40]. The technique is elaborated in the Normal Form
Theorem III of [22, Theorem III.13]. For the quadratic terms there is no term
which has a pre-factor which is constant in time and so all quadratic terms
can be eliminated by a transformation

$\displaystyle A$
$\displaystyle=A_{2}+\sum_{j=0}^{2}\sum_{m\in{\mathbb{Z}}}g_{A,2,j,m,1}A_{2}^{j}B_{2}^{2-j}e^{2im\pi
t/T}e^{i\pi t/T},$ $\displaystyle B$
$\displaystyle=B_{2}+\sum_{j=0}^{2}\sum_{m\in{\mathbb{Z}}}g_{B,2,j,m,1}A_{2}^{j}B_{2}^{2-j}e^{2im\pi
t/T}e^{i\pi t/T}.$

By suitably choosing the coefficients $g_{A,2,j,m,1}$ and $g_{B,2,j,m,1}$ we
find

$\displaystyle\dot{A}_{2}$
$\displaystyle=B_{2}+\sum_{n=3}^{\infty}\sum_{j=0}^{n}\sum_{m\in{\mathbb{Z}}}f_{A,n,j,m,1}A_{2}^{j}B_{2}^{n-j}e^{2im\pi
t/T}e^{i(n-1)\pi t/T},$ $\displaystyle\dot{B}_{2}$
$\displaystyle=\widetilde{\varepsilon}^{2}A_{2}+\sum_{n=3}^{\infty}\sum_{j=0}^{n}\sum_{m\in{\mathbb{Z}}}f_{B,n,j,m,1}A_{2}^{j}B_{2}^{n-j}e^{2im\pi
t/T}e^{i(n-1)\pi t/T},$

with new time-independent coefficients $f_{A,n,j,m,1}$ and $f_{B,n,j,m,1}$.
For simplifying the cubic terms we make a near identity transformation

$\displaystyle A_{2}$
$\displaystyle=A_{3}+\sum_{j=0}^{3}\sum_{m\in{\mathbb{Z}}}g_{A,3,j,m,2}A_{3}^{j}B_{3}^{3-j}e^{2im\pi
t/T},$ $\displaystyle B_{2}$
$\displaystyle=B_{3}+\sum_{j=0}^{3}\sum_{m\in{\mathbb{Z}}}g_{B,3,j,m,2}A_{3}^{j}B_{3}^{3-j}e^{2im\pi
t/T}.$

Again by suitably choosing the coefficients $g_{A,3,j,m,2}$ and
$g_{B,3,j,m,2}$ we find

$\displaystyle\dot{A}_{3}$
$\displaystyle=B_{3}+f_{A,3,3,-1,2}A_{3}^{3}+f_{A,3,2,-1,2}A_{3}^{2}B_{3}+f_{A,3,1,-1,2}A_{3}B_{3}^{2}+f_{A,3,0,-1,2}B_{3}^{3}$
$\displaystyle\quad+\sum_{n=4}^{\infty}\sum_{j=0}^{n}\sum_{m\in{\mathbb{Z}}}f_{A,n,j,m,2}A_{3}^{j}B_{3}^{n-j}e^{2im\pi
t/T}e^{i(n-1)\pi t/T},$ $\displaystyle\dot{B}_{3}$
$\displaystyle=\widetilde{\varepsilon}^{2}A_{3}+f_{B,3,3,-1,2}A_{3}^{3}+f_{B,3,2,-1,2}A_{3}^{2}B_{3}+f_{B,3,1,-1,2}A_{3}B_{3}^{2}+f_{B,3,0,-1,2}B_{3}^{3}$
$\displaystyle\quad+\sum_{n=4}^{\infty}\sum_{j=0}^{n}\sum_{m\in{\mathbb{Z}}}f_{B,n,j,m,2}A_{3}^{j}B_{3}^{n-j}e^{2im\pi
t/T}e^{i(n-1)\pi t/T},$

with new coefficients $f_{A,n,j,m,2}$ and $f_{B,n,j,m,2}$.

##  5 Proof of Theorem 1.1

Here we obtain the oscillating homoclinic solutions with small tails for
$c=0$. The bifurcating solutions scale as
$A_{3}(t)=\widetilde{\varepsilon}\widetilde{A}(\tau)$ and
$B_{3}(t)=\widetilde{\varepsilon}^{2}\widetilde{B}(\tau)$, with
$\tau=\widetilde{\varepsilon}t$. For the rescaled variables we find

(46a) $\displaystyle\partial_{\tau}\widetilde{A}$ $\displaystyle=$
$\displaystyle\widetilde{B}+\mathcal{O}(\widetilde{\varepsilon}),$ (46b)
$\displaystyle\partial_{\tau}\widetilde{B}$ $\displaystyle=$
$\displaystyle\widetilde{A}+f_{B,3,3,-1,2}\widetilde{A}^{3}+\mathcal{O}(\widetilde{\varepsilon}).$

Ignoring the terms of order $\mathcal{O}(\widetilde{\varepsilon})$ we find two
homoclinic solutions to the origin, see the left panel of Figure 4, if the
following sign condition holds.

(Coeff) Assume that

$f_{B,3,3,-1,2}<0.$

######  Remark 5.1.

It is shown in Section 7 that (Coeff) can generally be satisfied either at the
bifurcation (23) or (24).

######  Remark 5.2.

The truncated system (46), with $\mathcal{O}(\widetilde{\varepsilon})$ terms
neglected, admits an explicit solution

(47) $\left\\{\begin{array}[]{l}\widetilde{A}_{\rm
hom,0}(\tau)=\sqrt{2|f_{B,3,3,-1,2}|^{-1}}{\rm sech}(\tau-\tau_{0}),\\\
\widetilde{B}_{\rm hom,0}(\tau)=-\sqrt{2|f_{B,3,3,-1,2}|^{-1}}{\rm
tanh}(\tau-\tau_{0}){\rm sech}(\tau-\tau_{0}),\end{array}\right.$

where $f_{B,3,3,-1,2}<0$ and $\tau_{0}\in\mathbb{R}$ is to our disposal.

In the invariant subspace $\\{W_{h,M}=0\\}$ the homoclinic orbits persist in
the reduced system (42) if the following reversibility condition holds.

(Rev) Assume that there exists $t_{0}\in[0,T]$ such that
$k(t-t_{0})=k(t_{0}-t)$.

######  Remark 5.3.

Condition (Rev) is satisfied for $k(t)$ defined in (5) with either
$t_{0}=\frac{1}{2}\tau_{d}T$ or $t_{0}=\frac{1}{2}(\tau_{d}+1)T$. There exist
many other possible time-periodic coefficients $k(t)$ satisfying the
reversibility condition (5.3) which lead to the same kind of solutions as in
Theorems 1.1.

As a consequence of the reversibility, if
$t\mapsto(\widetilde{A}(t),\widetilde{B}(t))$ is a solution, so is

$t\mapsto(\widetilde{A}(2t_{0}-t),-\widetilde{B}(2t_{0}-t)).$

Hence, the one-dimensional unstable manifold of the origin of the time
$2T$-mapping transversely intersects the fixed space of reversibility
$\\{(\widetilde{A},\widetilde{B}):\widetilde{B}=0\\}$ and continue from the
upper half to the lower half of the phase plane, see the right panel of Figure
4. Hence, by extending $(\widetilde{A}(t),\widetilde{B}(t))_{t\leq t_{0}}$ by
its mirror picture $(\widetilde{A}(2t_{0}-t),-\widetilde{B}(2t_{0}-t))_{t\geq
t_{0}}$ at the fixed space of reversibility we constructed a homoclinic orbit

$(A,B)(t)=(\widetilde{\varepsilon}\widetilde{A}_{\rm hom}(\varepsilon
t),\widetilde{\varepsilon}^{2}\widetilde{B}_{\rm
hom}(\widetilde{\varepsilon}t))$

for the reduced system (44), where $\widetilde{A}_{\rm hom}=\widetilde{A}_{\rm
hom,0}+\mathcal{O}(\widetilde{\varepsilon})$ and $\widetilde{B}_{\rm
hom}=\widetilde{B}_{\rm hom,0}+\mathcal{O}(\widetilde{\varepsilon})$ with the
leading-order solution given by (47). In the original variables the homoclinic
orbit is denoted with $\mathcal{W}$ and corresponds to the truncation
$W_{h,M}=0$ for a given $M\in\mathbb{N}$.

(a) (b)

Figure 4: (a) Two homoclinic solutions of Eq. (46) with
$\mathcal{O}(\widetilde{\varepsilon})$ terms neglected. (b) A sketch of the
transversal intersection of the unstable manifold with the fixed space of
reversibility.

The rest of this section contains the proof of Theorem 1.1, which we rewrite
in notations of Sections 3 and 4 as follows.

######  Theorem 5.4.

Assume the validity of (Spec), (Coeff), and (Rev). Then there exists an
$\widetilde{\varepsilon}_{0}>0$ and $C_{0}>0$ such that for all
$\widetilde{\varepsilon}\in(0,\widetilde{\varepsilon}_{0})$ and every
$M\in\mathbb{N}$ with $M\geq 3$ the system (35)-(36) possesses a generalized
homoclinic solution $W_{\rm
hom}:[-\widetilde{\varepsilon}^{-M+2},\widetilde{\varepsilon}^{-M+2}]\to{\mathbb{R}}^{2N}$
with

$\sup_{t\in[-\widetilde{\varepsilon}^{-M+2},\widetilde{\varepsilon}^{-M+2}]}\|W_{\rm
hom}(t)-\mathcal{W}(t)\|\leq C_{0}\widetilde{\varepsilon}^{M-1}$

with $\lim\limits_{|t|\to\infty}\mathcal{W}(t)=0$. Moreover, for
$\mathcal{W}(t)=(\mathcal{W}_{0}(t),\mathcal{W}_{h}(t))$ we have

$\sup_{t\in[-\widetilde{\varepsilon}^{-M+2},\widetilde{\varepsilon}^{-M+2}]}\|\mathcal{W}_{h}(t)\|\leq
C_{0}\widetilde{\varepsilon}^{2}$

and

$\sup_{t\in[-\widetilde{\varepsilon}^{-M+2},\widetilde{\varepsilon}^{-M+2}]}\|\mathcal{W}_{0}(t)-\widetilde{\varepsilon}A_{\rm
hom,0}(\widetilde{\varepsilon}^{2}t)\varphi_{1}\|\leq
C_{0}\widetilde{\varepsilon}^{2},$

with $A_{\rm hom,0}$ given by (47).

######  Proof 5.5.

To prove persistence of $\mathcal{W}$ if the terms
$H_{M}(W_{0})=\mathcal{O}(\|W_{0}\|^{M})$ are taken into account, we use again
the reversibility. Obviously, it is impossible that the one-dimensional
unstable manifold transversally intersects the $N$-dimensional fixed space of
reversibility for the full system (37)-(38). Therefore, it can only be
expected that the homoclinic solutions persist as solutions with small tails
for $|t|\to\infty$. This can rigorously be shown by intersecting the
$2N-1$-dimensional center-unstable manifold with the fixed space of
reversibility. Obviously this intersection is a transversal intersection.

On the center-unstable manifold the solutions converge towards the center
manifold for $t\to-\infty$ with some exponential rate. However, the solutions
on the center manifold can grow slowly, hence it remains to obtain bounds for
such solutions. In a first step we apply another normal form transformation

$W_{0,M}=\widetilde{W}_{0,M}+Q_{0}(W_{0},W_{h,M}),\qquad
W_{h,M}=\widetilde{W}_{h,M}+Q_{h}(W_{0},W_{h,M})$

to eliminate the bilinear terms which are linear in $W_{0}$ and linear in
$W_{h,M}$ from the full system (37)-(38), where the

$Q_{0,h}(W_{0},W_{h,M})=\sum_{m\in{\mathbb{Z}}}Q_{0,h}(W_{0},W_{h,M})[m]e^{2im\pi
t/T}e^{i\pi t/T}$

are bilinear mappings in their arguments. With some slight abuse of notation
we skip the tildes and reconsider the system (37)-(38) but now with $N_{0,M}$
and $N_{h,M}$ additionally satisfying

$\displaystyle\|N_{0,M}(W_{0},W_{h,M})\|$ $\displaystyle\leq$ $\displaystyle
C(\|W_{0}\|^{3}+\|W_{0}\|^{2}\|W_{h,M}\|+\|W_{h,M}\|^{2}),$
$\displaystyle\|N_{h,M}(W_{0},W_{h,M})\|$ $\displaystyle\leq$ $\displaystyle
C(\|W_{0}\|^{2}\|W_{h,M}\|+\|W_{h,M}\|^{2}).$

The transformations are possible due to the spectral assumption (Spec).

In a second step we introduce the deviation $\widetilde{W}_{0}$ from the
homoclinic orbit $\mathcal{W}_{0}$ by
$W_{0}=\mathcal{W}_{0}+\widetilde{W}_{0}$. The subsequent estimates on the
deviation $\widetilde{W}_{0}$ have already been carried out in a number of
papers, cf. [17, 10]. We use the cutoff functions to estimate the solutions on
$[-\xi_{0},\xi_{0}]$ with a suitable chosen large $\xi_{0}$ as
$\tilde{\varepsilon}\to 0$.

The homoclinic orbit can be estimated pointwise by
$\|\mathcal{W}_{0}(t)\|\leq\widetilde{\varepsilon}q(\widetilde{\varepsilon}t)$
with a smooth $q\in L^{1}({\mathbb{R}})$. We denote the stable part of
$\widetilde{W}_{0}$ by $\widetilde{W}_{0,s}$ and the projection on the stable
part by $\Pi_{s}$. We find for a large $\xi_{0}$ determined below that

$\displaystyle\|\widetilde{W}_{0,s}(t)\|$
$\displaystyle=\left\|\int_{-\xi_{0}}^{t}e^{\Lambda_{0}(t-\tau)}\Pi_{0}\left[N_{0,M}(\mathcal{W}_{0}+\widetilde{W}_{0},W_{h,M})-N_{0,M}(\mathcal{W}_{0},0)\right](\tau)d\tau\right\|$
$\displaystyle\leq
CY(t)\int_{-\infty}^{t}e^{-\widetilde{\varepsilon}(t-\tau)}\widetilde{\varepsilon}^{2}q^{2}(\widetilde{\varepsilon}\tau)d\tau+CY(t)^{2}\int_{-\infty}^{t}e^{-\widetilde{\varepsilon}(t-\tau)}d\tau$
$\displaystyle\leq
C\widetilde{\varepsilon}Y(t)+C\widetilde{\varepsilon}^{-1}Y(t)^{2},$

where

$Y(t):=\sup_{\tau\in[-\xi_{0},t]}(\|\widetilde{W}_{0}(\tau)\|+\|W_{h,M}(\tau)\|).$

If the solution is in the fixed space of reversibility at $t=0$ we find

$\displaystyle\|W_{h,M}(t)\|$
$\displaystyle=\left\|\int_{0}^{t}e^{\Lambda_{h}(t-\tau)}(N_{h,M}(W_{0},W_{h,M})(\tau)+H_{M}(W_{0})(\tau))d\tau\right\|$
$\displaystyle\leq
CY(t)\int_{0}^{t}\widetilde{\varepsilon}^{2}q^{2}(\widetilde{\varepsilon}\tau)d\tau+C|t|Y(t)^{2}+C\int_{0}^{t}\widetilde{\varepsilon}^{M}q^{M}(\widetilde{\varepsilon}\tau)d\tau$
$\displaystyle\leq
C\widetilde{\varepsilon}Y(t)+C|t|Y(t)^{2}+C\widetilde{\varepsilon}^{M-1}.$

Summarizing the estimates yields

$Y(\xi_{0})\leq
C\widetilde{\varepsilon}Y(\xi_{0})+C\widetilde{\varepsilon}^{-1}Y(\xi_{0})^{2}+C|\xi_{0}|Y(\xi_{0})^{2}+C\widetilde{\varepsilon}^{M-1}$

and so $Y(\xi_{0})\leq C\varepsilon^{M-1}$ for
$\xi_{0}\leq\widetilde{C}\varepsilon^{-M+2}$ and $M\geq 3$. This completes the
proof of the theorem in view of the normal form transformations.

##  6 Proof of Theorem 1.2 and Theorem 1.3

Here we consider the case of small damping $c>0$ where reversibility no longer
holds. For consistency of our analysis, we assume that
$c=\widetilde{c}\widetilde{\varepsilon}$ with $\widetilde{c}=\mathcal{O}(1)>0$
fixed. With the same analysis as for (44) and (46) we end up in
$\\{W_{h,M}=0\\}$ at the reduced system for the rescaled variables which is
now given by

(48a) $\displaystyle\partial_{\tau}\widetilde{A}$ $\displaystyle=$
$\displaystyle\widetilde{B}+\mathcal{O}(\widetilde{\varepsilon}),$ (48b)
$\displaystyle\partial_{\tau}\widetilde{B}$ $\displaystyle=$
$\displaystyle\widetilde{A}-\frac{\widetilde{c}}{\underline{m}}\widetilde{B}+f_{B,3,3,-1,2}\widetilde{A}^{3}+\mathcal{O}(\widetilde{\varepsilon}),$

where we used the expression $\gamma=\upsilon-\frac{c}{2\underline{m}}$ for
the Floquet exponents of the mapping (11) in Section 2.1. Ignoring the terms
of order $\mathcal{O}(\widetilde{\varepsilon})$ we first find two fixed points
$(\widetilde{A},\widetilde{B})_{\pm}=(\pm 1/\sqrt{|f_{B,3,3,-1,2}}|,0)$ if
$f_{B,3,3,-1,2}<0$. Since the fixed points are hyperbolic they persist as
fixed points of the time-$2T$-mapping even if the terms of order
$\mathcal{O}(\widetilde{\varepsilon})$ are included. This proves Theorem 1.2.

(a) (b)

Figure 5: (a) Two heteroclinic solutions for small $c>0$ of Eq. (48) with
$\mathcal{O}(\widetilde{\varepsilon})$ terms neglected. (b) The heteroclinic
solution of Eq. (48) in the intersection of the one-dimensional unstable
manifold from the zero equilibrium point (in black) and the stable manifold
from the positive fixed point (in gray).

Next we find two heteroclinic solutions connecting the origin with the two
fixed points $(\widetilde{A},\widetilde{B})_{\pm}$. In order to prove the
persistence of these heteroclinic solutions we use that the one-dimensional
unstable manifold of the origin transversely intersects the two-dimensional
stable manifold of the two other fixed points of the reduced system. These
heteroclinic solutions are shown on Figure 5.

To prove the persistence of heteroclinic solutions if the terms
$H_{M}(W_{0})=\mathcal{O}(|W_{0}|^{M})$ are taken into account we use that the
stable manifold of the two other fixed points is $2N$-dimensional and so the
one-dimensional unstable manifold of the origin transversally intersects the
$2N$-dimensional stable manifold of the two anti-periodic solutions from
Theorem 1.2. Theorem 1.3 is proven by using the transformations of Sections 3
and 4, with the decomposition $W=(W_{0},W_{h})$.

##  7 Multiple-scale analysis and checking the assumption (Coeff)

For the system (1) with (2) we already checked the spectral condition (Spec)
used in Theorems 1.1-1.3. Here we establish the validity of the normal form
coefficient condition (Coeff).

We do so by formally deriving the reduced systems (46) and (48) via a
multiple-scale analysis, which will yield an explicit and convenient formula
for $f_{B,3,3,-1,2}$ after adjusting the notations. Truncating (48) with
$\mathcal{O}(\widetilde{\varepsilon})$ terms neglected yields the scalar
equation:

(49)
$\partial_{\tau}^{2}\widetilde{A}+\frac{\widetilde{c}}{\underline{m}}\partial_{\tau}\widetilde{A}=\widetilde{A}+f_{B,3,3,-1,2}\widetilde{A}^{3}$

######  Remark 7.1.

The scalar equation (49) is recovered in Eq. (66) below with the
correspondence between $\varepsilon$ and $\widetilde{\varepsilon}$ given by
(45). Note that the definition
$A=\tilde{\varepsilon}\widetilde{A}(\tilde{\varepsilon}t)$, see Eq. (43), is
different from the definition $A(\varepsilon t)$, see Eq. (61) below. For Eq.
(49), the amplitude $\widetilde{A}$ is introduced based on the Floquet
theorem, the decomposition into subspaces, and the normal form
transformations. For Eq. (66), the amplitude $A$ is introduced directly by the
perturbation expansions in powers of $\varepsilon$, see Eq. (60).

We define as in (25),

(50) $k(t)=k_{0}(t)+\delta\varepsilon^{2},$

where $k_{0}(t)$ is the potential for the bifurcation (23) or (24),
$\varepsilon$ is a small parameter for the asymptotic expansion, and $\delta$
is a proper sign factor.

###  7.1 Bloch modes of the unperturbed linear problem

We start with the linear problem for $c=0$. In order to construct the
asymptotic expansion, we define the Bloch function $v(t)=e^{i\ell
t}f_{j}(\ell,t)$ for the $j$-th spectral band $\lambda_{j}(\ell)$ of
$\mathcal{L}_{0}v=\lambda v$, where $f_{j}(\ell,t+T)=f_{j}(\ell,t)$. The
parameter $\ell$ is defined in the Brillouin zone $[0,\frac{2\pi}{T})$ with
$\ell=0$ and $\ell=\frac{\pi}{T}$ being at the ends of each spectral band
corresponding to periodic and anti-periodic solutions respectively.

Let $L^{2}_{\rm per}$ be the Hilbert space of $2\pi$-periodic functions with
the inner product $\langle\cdot,\cdot\rangle$ and the induced norm
$\|\cdot\|_{L^{2}_{\rm per}}$. The $L^{2}_{\rm per}$-normalized Bloch function
$f_{j}(\ell,t)$ is a $T$-periodic solution of the spectral problem

(51)
$\left[-\underline{m}(\partial_{t}+i\ell)^{2}-k_{0}(t)\right]f_{j}(\ell,t)=\lambda_{j}(\ell)f_{j}(\ell,t),$

which can be differentiated in $\ell$ as

(52)
$\displaystyle\left[-\underline{m}(\partial_{t}+i\ell)^{2}-k_{0}(t)-\lambda_{j}(\ell)\right]\partial_{\ell}f_{j}(\ell,t)=2i\underline{m}(\partial_{t}+i\ell)f_{j}(\ell,t)+\lambda_{j}^{\prime}(\ell)f_{j}(\ell,t)$

and

$\displaystyle\left[-\underline{m}(\partial_{t}+i\ell)^{2}-k_{0}(t)-\lambda_{j}(\ell)\right]\partial^{2}_{\ell}f_{j}(\ell,t)=4i\underline{m}(\partial_{t}+i\ell)\partial_{\ell}f_{j}(\ell,t)$
(53)
$\displaystyle\qquad+2\lambda_{j}^{\prime}(\ell)\partial_{\ell}f_{j}(\ell,t)+\left[\lambda_{j}^{\prime\prime}(\ell)-2\underline{m}\right]f_{j}(\ell,t),$

Projecting to $f(\ell,t)$ in $L^{2}_{\rm per}$ yields from (52) and (53):

(54) $\displaystyle\lambda_{j}^{\prime}(\ell)$
$\displaystyle={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-}2\underline{m}\langle
f_{j}(\ell,\cdot),i\partial_{t}f_{j}(\ell,\cdot)\rangle{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}+}2\underline{m}\ell,$

and

$\displaystyle\lambda_{j}^{\prime\prime}(\ell)$
$\displaystyle+2\lambda_{j}^{\prime}(\ell)\langle
f_{j}(\ell,\cdot),\partial_{\ell}f_{j}(\ell,\cdot)\rangle=-4\underline{m}\langle
f_{j}(\ell,\cdot),i\partial_{t}\partial_{\ell}f_{j}(\ell,\cdot)\rangle$ (55)
$\displaystyle\qquad+4\underline{m}\langle
f_{j}(\ell,\cdot),\partial_{\ell}f_{j}(\ell,\cdot)\rangle\langle
f_{j}(\ell,\cdot),i\partial_{t}f_{j}(\ell,\cdot)\rangle+2\underline{m},$

where the normalization condition $\|f_{j}(\ell,\cdot)\|_{L^{2}_{\rm per}}=1$
has been used.

######  Remark 7.2.

If $\ell_{0}=\frac{\pi}{T}$ for the bifurcating mode and the band gap
$(\mu_{1},\mu_{2})$ has a non-zero width, then necessarily,
$\lambda_{j}^{\prime}(\ell_{0})=0$. If the bifurcation (23) occurs at
$\mu_{1}$, then $j=1$ and $\lambda_{1}^{\prime\prime}(\ell_{0})<0$. See Fig.
3(c) for an example. If the bifurcation (24) occurs at $\mu_{2}$, then $j=2$
and $\lambda_{2}^{\prime\prime}(\ell_{0})>0$. See Fig. 8(a).

###  7.2 Perturbation of the linear problem

Let us consider asymptotic solutions of the linear equation

(56)
$(\mathcal{L}_{0}-K_{2}\omega^{2}(q_{m_{0}})-\delta\varepsilon^{2})\widehat{u}(t)=0,$

which follows from (7) and (50) at $m=m_{0}$. As in (23), we take
$\mu_{1}=K_{2}\omega^{2}(q_{m_{0}})$ for which $\ell_{0}=\frac{\pi}{T}$ is
selected in the first spectral band
$\\{\lambda_{1}(\ell)\\}_{\ell\in[0,\frac{2\pi}{T})}$. The set-up for the
second spectral band $\\{\lambda_{2}(\ell)\\}_{\ell\in[0,\frac{2\pi}{T})}$
when $\mu_{2}=K_{2}\omega^{2}(q_{m_{0}})$ as in (24) is essentially identical
(see remark 7.3). Expanding

$\widehat{u}(t)=A(\varepsilon t)e^{i\ell_{0}t}f_{1}(\ell_{0},t)+\varepsilon
B(\varepsilon
t)e^{i\ell_{0}t}g_{1}(\ell_{0},t)+\varepsilon^{2}e^{i\ell_{0}t}C(\varepsilon
t)h_{1}(\ell_{0},t)+\mathcal{O}(\varepsilon^{3}),$

with $A,B,C$ and $g_{1},h_{1}$ to be determined, we obtain for
$\mu_{1}=K_{2}\omega^{2}(q_{m_{0}})$ at the order of
$\mathcal{O}(\varepsilon)$ that

$\displaystyle
B\left[-\underline{m}(\partial_{t}+i\ell_{0})^{2}-k_{0}(t)-\mu_{1}\right]g_{1}(\ell_{0},t)=2mA^{\prime}(\partial_{t}+i\ell_{0})f_{1}(\ell_{0},t),$

Since $\lambda_{1}(\ell_{0})=\mu_{1}$ and $\lambda_{1}^{\prime}(\ell_{0})=0$,
comparing with (52) yields

$g_{1}(t)=\partial_{\ell}f_{1}(\ell_{0},t)\quad\mbox{\rm and}\quad
B(\varepsilon t)=-iA^{\prime}(\varepsilon t).$

At the order of $\mathcal{O}(\varepsilon^{2})$, we obtain the linear
inhomogeneous equation,

$\displaystyle
C\left[-\underline{m}(\partial_{t}+i\ell_{0})^{2}-k_{0}(t)-\mu_{1}\right]h_{1}(t)$
$\displaystyle=-2\underline{m}iA^{\prime\prime}(\partial_{t}+i\ell_{0})\partial_{\ell}f_{1}(\ell_{0},t)$
(57) $\displaystyle\quad+\underline{m}A^{\prime\prime}f_{1}(\ell_{0},t)+\delta
Af_{1}(\ell_{0},t).$

Comparing (57) with (53) yields
$h_{1}(t)=\partial_{\ell}^{2}f_{1}(\ell_{0},t)$ and $C(\varepsilon
t)=-\frac{1}{2}A^{\prime\prime}(\varepsilon t)$ if and only if $A(\tau)$,
$\tau:=\varepsilon t$ satisfies the amplitude equation

(58)
$\frac{1}{2}\lambda_{1}^{\prime\prime}(\ell_{0})A^{\prime\prime}(\tau)+\delta
A(\tau)=0.$

Alternatively, the amplitude equation (58) can be obtained by projecting the
linear inhomogeneous equation (57) to $f_{1}(\ell_{0},t)$ in $L^{2}_{\rm per}$
and using equations (54) and (55) with $\lambda_{1}^{\prime}(\ell_{0})=0$.

Since $\delta$ and $\lambda_{1}^{\prime\prime}(\ell_{0})$ have opposite signs,
$A(\tau)$ of Eq. (58) will experience exponential growth with rate

(59)
$\varepsilon\frac{\sqrt{2}}{\sqrt{|\lambda_{1}^{\prime\prime}(\ell_{0})|}},$

which is the leading order for Eq. (45). This is an approximation of the real
part of the Floquet exponent associated to the $m_{0}$-th mode of the linear
problem with $k(t)=k_{0}(t)+\delta\varepsilon^{2}$, see Eq. (26) and Fig. 6(c)
for example.

######  Remark 7.3.

If $\lambda_{1}^{\prime\prime}(\ell_{0})<0$ for the bifurcating mode (23),
then $\delta=+1$ is selected from the condition that
$\mu_{1}=K_{2}\omega^{2}(q_{m_{0}})$ is inside the band gap of the perturbed
operator $\mathcal{L}=\mathcal{L}_{0}-\delta\varepsilon^{2}$. On the other
hand, if $\lambda_{2}^{\prime\prime}(\ell_{0})>0$ for the bifurcating mode
(24), then $\delta=-1$ is selected from the condition that
$\mu_{2}=K_{2}\omega^{2}(q_{m_{0}})$ is inside the band gap of $\mathcal{L}$.

###  7.3 Perturbation of the nonlinear problem

Let us now consider asymptotic solutions of the nonlinear equation (1) with
(50) and $c=\varepsilon\tilde{c}$, where $\tilde{c}\geq 0$ is fixed and
$\delta=-{\rm sign}(\lambda_{1}^{\prime\prime}(\ell_{0}))$. As in (23), we
take $\mu_{1}=K_{2}\omega^{2}(q_{m_{0}})$, for which $\ell_{0}=\frac{\pi}{T}$
is selected in the first spectral band
$\\{\lambda_{1}(\ell)\\}_{\ell\in[0,\frac{2\pi}{T})}$. To simplify
computations, we will use expansions in terms of real-valued functions only.
Expanding

(60) $u_{n}(t)=\varepsilon
U_{n}^{(1)}(t)+\varepsilon^{2}U_{n}^{(2)}(t)+\varepsilon^{3}U_{n}^{(3)}(t)+\mathcal{O}(\varepsilon^{4}),$

we select the leading order in the form

(61) $U_{n}^{(1)}(t)=A(\varepsilon t)g_{1}(t)\sin(q_{m_{0}}n),$

where the amplitude $A(\varepsilon t)$ is real according to (49) and the
eigenfunction

(62)
$g_{1}(t):=\left[e^{i\ell_{0}t}f_{1}(\ell_{0},t)+e^{-i\ell_{0}t}\bar{f}_{1}(\ell_{0},t)\right]$

is a real-valued, $T$-anti-periodic solution of
$\mathcal{L}_{0}g_{1}=\mu_{1}g_{1}$ satisfying $g_{1}(t+T)=-g_{1}(t)$. At the
order of $\mathcal{O}(\varepsilon^{2})$ we obtain

$\displaystyle\mathcal{L}_{0}U_{n}^{(2)}+K_{2}(U_{n+1}^{(2)}-2U_{n}^{(2)}+U_{n-1}^{(2)})=H_{n}^{(2)},$

with

$\displaystyle H_{n}^{(2)}$
$\displaystyle=2\underline{m}\partial_{\tau}\partial_{t}U_{n}^{(1)}+\tilde{c}\partial_{t}U^{(1)}_{n}+K_{3}\left[(U_{n+1}^{(1)}-U_{n}^{(1)})^{2}-(U_{n}^{(1)}-U_{n-1}^{(1)})^{2}\right],$

where $\tau=\varepsilon t$ and $K_{3}$ is the coefficient of the quadratic
term in (1). In the explicit form, we obtain

$\displaystyle H_{n}^{(2)}$
$\displaystyle=[2\underline{m}A^{\prime}(\tau)+\tilde{c}A(\tau)]g_{1}^{\prime}(t)\sin(q_{m_{0}}n)+K_{3}A(\tau)^{2}g_{1}(t)^{2}F_{n}^{(2)},$

with

$\displaystyle F_{n}^{(2)}$
$\displaystyle=\left[\sin(q_{m_{0}}(n+1))-\sin(q_{m_{0}}n)\right]^{2}-\left[\sin(q_{m_{0}}n)-\sin(q_{m_{0}}(n-1))\right]^{2}$
$\displaystyle=-2\sin(q_{m_{0}})(1-\cos(q_{m_{0}}))\sin(2q_{m_{0}}n).$

The solution for $U_{n}^{(2)}(t)$ can be written in the form

$\displaystyle U_{n}^{(2)}(t)$
$\displaystyle=\left[A^{\prime}(\tau)+\frac{\tilde{c}}{2\underline{m}}\right]h_{1}(t)\sin(q_{m_{0}}n)+K_{3}A(\tau)^{2}h_{2}(t)\sin(2q_{m_{0}}n),$

where $h_{1}$ and $h_{2}$ are solutions of the linear inhomogeneous equations:

(63) $\displaystyle(\mathcal{L}_{0}-K_{2}\omega^{2}(q_{m_{0}}))h_{1}$
$\displaystyle=2\underline{m}g_{1}^{\prime}(t),$ (64)
$\displaystyle(\mathcal{L}_{0}-K_{2}\omega^{2}(2q_{m_{0}}))h_{2}$
$\displaystyle=-2\sin(q_{m_{0}})(1-\cos(q_{m_{0}}))g_{1}(t)^{2}.$

It follows from the linear theory that the real, $T$-anti-periodic solution
for $h_{1}(t)$ exists in the form:

$h_{1}(t)=-ie^{i\ell_{0}t}\partial_{\ell}f_{1}(\ell_{0},t)+ie^{-i\ell_{0}t}\partial_{\ell}\bar{f}_{1}(\ell_{0},t).$

There exists a unique $T$-periodic solutions for $h_{2}(t)$ if and only if the
non-resonance condition is met:

(65) $K_{2}\omega^{2}(2q_{m_{0}})\notin\cup_{j=1}^{\infty}\lambda_{j}(0).$

This non-resonance condition is satisfied if the spectral assumption (Spec) is
satisfied.

At the order of $\mathcal{O}(\varepsilon^{3})$, we obtain

$\displaystyle\mathcal{L}_{0}U_{n}^{(3)}+K_{2}(U_{n+1}^{(3)}-2U_{n}^{(3)}+U_{n-1}^{(3)})=H_{n}^{(3)},$

with

$\displaystyle H_{n}^{(3)}$ $\displaystyle=\delta
U_{n}^{(1)}+2m\partial_{\tau}\partial_{t}U_{n}^{(2)}+m\partial_{\tau}^{2}U_{n}^{(1)}+\tilde{c}\partial_{t}U_{n}^{(2)}+\tilde{c}\partial_{\tau}U_{n}^{(1)}$
$\displaystyle\quad+2K_{3}\left[(U_{n+1}^{(1)}-U_{n}^{(1)})(U_{n+1}^{(2)}-U_{n}^{(2)})-(U_{n}^{(1)}-U_{n-1}^{(1)})(U_{n}^{(2)}-U_{n-1}^{(2)})\right]$
$\displaystyle\quad-
K_{4}\left[(U_{n+1}^{(1)}-U_{n}^{(1)})^{3}-(U_{n}^{(1)}-U_{n-1}^{(1)})^{3}\right],$

where $K_{4}$ is the coefficient of the cubic term in (1).

By using Euler’s formulas, we obtain

$\displaystyle\left[\sin(q_{m_{0}}(n+1))-\sin(q_{m_{0}}n)\right]^{3}-\left[\sin(q_{m_{0}}n)-\sin(q_{m_{0}}(n-1))\right]^{3}$
$\displaystyle\qquad=-3(1-\cos(q_{m_{0}}))^{2}\sin(q_{m_{0}}n)$
$\displaystyle\qquad\qquad+\frac{1}{2}\left[1-3\cos(q_{m_{0}})+3\cos(2q_{m_{0}})-\cos(3q_{m_{0}})\right]\sin(3q_{m_{0}}n)$

and

$\displaystyle\left[\sin(q_{m_{0}}(n+1))-\sin(q_{m_{0}}n)\right]\left[\sin(2q_{m_{0}}(n+1))-\sin(2q_{m_{0}}n)\right]$
$\displaystyle\qquad\qquad-\left[\sin(q_{m_{0}}n)-\sin(q_{m_{0}}(n-1))\right]\left[\sin(2q_{m_{0}}n)-\sin(2q_{m_{0}}(n-1))\right]$
$\displaystyle\qquad=-\left[2\sin(q_{m_{0}})-\sin(2q_{m_{0}})\right]\sin(q_{m_{0}}n)$
$\displaystyle\qquad\qquad-\left[\sin(q_{m_{0}})+\sin(2q_{m_{0}})-\sin(3q_{m_{0}})\right]\sin(3q_{m_{0}}n)$

Hence, we obtain

$\displaystyle
H_{n}^{(3)}=D_{1}(t)\sin(q_{m_{0}}n)+D_{2}(t)\sin(2q_{m_{0}}n)+D_{3}(t)\sin(3q_{m_{0}}n),$

where we are only interested to write explicitly the coefficient for the
bifurcating mode:

$\displaystyle D_{1}(t)$ $\displaystyle=\delta
A(\tau)g_{1}(t)+mA^{\prime\prime}(\tau)g_{1}(t)+2mA^{\prime\prime}(\tau)h_{1}^{\prime}(t)$
$\displaystyle\quad{+\tilde{c}\left[A^{\prime}(\tau)g_{1}(t)+2A^{\prime}(\tau)h_{1}^{\prime}(t)\right]+\frac{\tilde{c}^{2}}{2\underline{m}}A(\tau)h_{1}^{\prime}(t)}$
$\displaystyle\quad+3K_{4}(1-\cos(q_{m_{0}}))^{2}A(\tau)^{3}g_{1}(t)^{3}-2K_{3}^{2}\left[2\sin(q_{m_{0}})-\sin(2q_{m_{0}})\right]A(\tau)^{3}g_{1}(t)h_{2}(t).$

Projecting $D_{1}(t)$ to $g_{1}(t)$ gives the amplitude equation for
$A(\tau)$:

(66)
$\frac{1}{2}\lambda_{1}^{\prime\prime}(\ell_{0})\left[A^{\prime\prime}(\tau)+\frac{\tilde{c}}{\underline{m}}A^{\prime}(\tau)+\frac{\tilde{c}^{2}}{4\underline{m}^{2}}A(\tau)\right]+\left[\delta-\frac{\tilde{c}^{2}}{4\underline{m}}\right]A(\tau)+\chi
A(\tau)^{3}=0,$

where

(67) $\lambda_{1}^{\prime\prime}(\ell_{0})=2m+4m\frac{\langle
g_{1},h_{1}^{\prime}\rangle}{\|g_{1}\|^{2}}$

and

(68) $\displaystyle\chi$
$\displaystyle=3K_{4}(1-\cos(q_{m_{0}}))^{2}\frac{\langle
g_{1}^{2},g_{1}^{2}\rangle}{\|g_{1}\|^{2}}-2K_{3}^{2}\left[2\sin(q_{m_{0}})-\sin(2q_{m_{0}})\right]\frac{\langle
g_{1}^{2},h_{2}\rangle}{\|g_{1}\|^{2}}.$

Since the linear part of (66) should be identical to the linear amplitude
equation (58), the new formula for $\lambda_{1}^{\prime\prime}(\ell_{0})$ must
be identical to the previous equation (55) with
$\lambda_{1}^{\prime}(\ell_{0})=0$. The expression for $\chi$ is defined in
terms of real quantities only.

######  Remark 7.4.

Equation (66) for $\tilde{c}=0$ is analogous to the stationary NLS equation
that can be derived in the context of spatially periodic media for the
description of breathers [21]. A similar equation was derived in [34] for a
(space-time continuous) photonic time crystal. Eq. (66) is equivalent to (49)
with the correspondence (45) and the appropriate definitions of amplitudes $A$
and $\widetilde{A}$. The coefficient $f_{B,3,3,-1,2}$ is constant proportional
to $\chi/\lambda_{1}^{\prime\prime}(\ell_{0})$.

######  Remark 7.5.

The coefficient in front of the second derivative in the amplitude equation
(66) comes from an expansion of the imaginary parts of the spectral curves at
the spectral gaps, see Fig. 3(c). The coefficient in front of the second
derivative changes sign at every spectral boundary. Since the coefficient in
front of the cubic coefficient in (66) does not change sign, the homoclinic
and heteroclinic solutions exist as bifurcating solutions at every spectral
gap associated with the anti-periodic eigenfunctions. In other words, if
$\chi<0$, then we pick the bifurcating mode at (23) with
$\lambda_{1}^{\prime\prime}(\ell_{0})<0$ and $\delta=+1$. If $\chi>0$, then we
pick the bifurcating mode at (24) with
$\lambda_{2}^{\prime\prime}(\ell_{0})>0$ and $\delta=-1$.

##  8 Comparison with numerical simulations

We now conduct a number of numerical simulations to illustrate the main
results of the paper. We start with the simplest case, and work up in
complexity.

Equation (66) with $\tilde{c}=0$ and $\delta=-{\rm
sign}(\lambda_{1}^{\prime\prime}(\ell_{0}))=-{\rm sign}(\chi)$ has the
following homoclinic solution

(69)
$A(\tau)=\sqrt{\frac{2}{|\chi|}}\operatorname{sech}\left(\sqrt{\frac{2}{|\lambda^{\prime\prime}_{1}(\ell_{0})|}}\tau\right)$

See Fig. 4(a) for an example plot of this solution in the $(A,A^{\prime})$
phase plane. Returning to ansatz Eq. (60), we have the following leading order
approximation in terms of the original lattice variables,

(70)
$u_{n}(t)=\varepsilon\sqrt{\frac{2}{|\chi|}}\operatorname{sech}\left(\sqrt{\frac{2}{|\lambda^{\prime\prime}_{1}(\ell_{0})|}}\varepsilon
t\right)g_{1}(t)\sin(q_{m_{0}}n).$

To make practical use of this approximation, the first step is to identify the
bifurcation value $k_{0}$ (namely the critical modulation amplitude values
$k_{a}^{0}$ and $k_{b}^{0}$) and critical mode number $m_{0}$ such that
$\mathrm{trace}(J)=-2$ where $J$ is the monodromy matrix defined in Eq. (20).
This corresponds to the bifurcation scenario shown in Fig. 3. In this case,
the Floquet exponent corresponding to $m_{0}$ will be purely imaginary and
will be of the form $i\ell_{0}=\frac{i\pi}{T}$. The corresponding Bloch mode
$e^{i\ell_{0}t}f_{1}(\ell_{0},t)$ is obtained by solving Eq. (51). This can be
done explicitly, see Sec. 2.1 or the appendix of [7] for details. Next, we
compute $\lambda^{\prime\prime}_{1}(\ell_{0})$ using Eq. (55) with $j=1$ which
depends only on $f_{1}$ and its derivatives, which can also be computed
explicitly. Equivalently, one can determine
$\lambda^{\prime\prime}_{1}(\ell_{0})$ using Eq. (67).

###  8.1 Examples with $c=0$ and $K_{3}=0$

The nonlinear coefficient $\chi$ can be computed from $g_{1}$ if $K_{3}=0$,
i.e., if there is no quadratic nonlinearity (the case of $K_{3}\neq 0$ is
discussed below). To compute $\chi$ in this case, we substitute
$g_{1}(t)=2{\rm Re}\left(e^{i\ell_{0}t}f_{1}(\ell_{0},t)\right)$ into Eq. (68)
and evaluate.

As our first example, we chose parameters that correspond to the spectral
picture in Fig. 3 and $K_{3}=0$ and $K_{4}=-0.8$. In this case the critical
mode $m_{0}=3$ lies at the top of the first spectral band, namely
$(\ell,\lambda)=(\pi/T,\mu_{1})$. This demonstrates that the spectral
condition (Spec) is satisfied and we choose $\delta=+1$. It can be seen from
Fig. 3(c), or via direct calculation, that
$\lambda^{\prime\prime}(\ell_{0})<0$. By choosing $K_{4}<0$, we have that
$\chi<0$, and thus (Coeff) is satisfied.

(a) (b) (c)

Figure 6: A generalized $q$-gap breather bifurcating from $\mu_{1}$. The
parameter values are $c=0$, $T=1/0.37$, $K_{3}=0$, $K_{4}=-0.8$ and
$K_{2}=\underline{m}=1$. The critical modulation amplitude parameters are
$k_{a}^{0}=0.5$ and $k_{b}^{0}=0.79$. The critical mode number is $m_{0}=3$,
which lies at the top of the first spectral band, namely
$\lambda(q_{3})=\mu_{1}$, (see Fig. 3(a)). (a) Numerical simulation with
initial value $u_{n}(0)=10^{-4}\sin(q_{3}n)$ and $\varepsilon=0.1$. The
displacement of the 6th particle is shown as a function of time. The
modulation amplitude is $k_{a}=k_{a}^{0}+\varepsilon^{2}=0.06$,
$k_{b}=k_{b}^{0}+\varepsilon^{2}=0.8$. The dashed line shows an approximation
of the envelope given by Eq. (70) with $\chi=-0.6530$,
$\lambda^{\prime\prime}(\ell_{0})=-21.0222$ and $\varepsilon=0.1$. For this
value of $\varepsilon$, the distance of the wavenumber to the edge of the gap
is $\Delta q=q_{3}-q_{\ell}=0.007$. (b) Intensity plot of the solution shown
in panel (a). Color intensity corresponds to $|u_{n}|$. (c) Plot of the
amplitude of the breather ($\mathrm{max}_{t}u_{6}(t)$) for the numerical
simulation (blue dots) and prediction based on Eq. (70) (blue line) as a
function of $\Delta q$. The real part of the Floquet exponent corresponding to
$q_{3}$ (solid red squares) and asymptotic approximation (59) (red line) are
also shown, which indicate the growth (decay) rate of the breather.

To generate the generalized $q$-gap breather, we keep all parameters fixed,
but select $\varepsilon=0.1$ and $k_{a}=k_{a}^{0}+\varepsilon^{2}=0.06$,
$k_{b}=k_{b}^{0}+\varepsilon^{2}=0.8$. With these parameter values, the
$m_{0}=3$ mode lies in the spectral gap. The corresponding Floquet multipliers
and exponents and are shown in Fig. 1(a,b). We initialize the numerical
simulation with

$u_{n}(0)=10^{-4}\sin(q_{3}n)\quad\mbox{\rm and}\quad\dot{u}_{n}=0.$

For initial data with such small amplitude, the dynamics will initially be
nearly linear, and hence the solution will grow exponentially with rate given
by Re($\gamma$), which is the real part the Floquet exponent associated to
mode $m=3$ (see the larger black dot in Fig. 1(a)). According to Eq. (58), an
approximation of this growth rate is given by (59) with $\varepsilon=0.1$. As
the amplitude increases in the dynamic evolution, the affect of the
nonlinearity comes into play, which will cause the solution to experience
decay, such that the resulting waveform is localized in time. The time series
of the $u_{6}(t)$ node is shown in Fig. 6(a). The temporal localization occurs
uniformly throughout the lattice, as seen in the intensity plot of Fig. 6(b).
By construction, the wavenumber $q_{3}$ lies in a wavenumber bandgap. Thus,
the solution shown in Fig. 6(a,b) is a generalized $q$-gap breather. The
analytical prediction based on Eq. (70) is shown as the dashed-line in panel
(a). For the sake of clarity, only the envelope of the approximation is shown,
which is simply a plot of the local maximums (and minimums) of Eq. (70).

(a) (b)

Figure 7: Numerical simulation with all parameters as in Fig. 6 but larger
values of $\varepsilon$. (a) With $\Delta q=0.032$ ($\varepsilon=0.22)$ a
temporally localized structure is initially observed, but over a longer time
interval a recurrence of breathing patterns emerges. The dashed line shows an
approximation of the envelope given by Eq. (70). (b) With $\Delta q=0.046$
($\varepsilon=0.26)$ the second breather like structure emerges before the
temporal localization of the first pulse can be achieved. The envelope
approximation given by Eq. (70) is still quite good until the breakdown of the
temporal localization.

For $\varepsilon=0.1$, the distance of the wavenumber to the edge of the gap
is $\Delta q=q_{\ell}-q_{3}=0.007$, where $q_{\ell}$ is the wavenumber at the
(left) edge of the gap. Recall from Eq. (27) that $\Delta
q=\mathcal{O}(\varepsilon^{2})$. Since the amplitude of the breather is
$\mathcal{O}(\varepsilon)$, see Eq. (70), the amplitude grows like
$\mathcal{O}(\sqrt{\Delta q})$. This observation was made numerically and
experimentally in [7], which we have now proved. This amplitude trend is
consistent with the trend found for discrete breathers in space-periodic
systems where it is well known that the breather amplitude grows like
$\mathcal{O}\left(\sqrt{\Delta\omega}\right)$, where $\Delta\omega$ is the
difference between the breather frequency and the edge of the frequency
spectrum [12]. Using (27) and Eq. (70) allows us to obtain an analytical
prediction of the breather amplitude dependence of the distance to the band
edge, see Fig. 6(c). For small values of $\Delta q$ (and hence $\varepsilon$)
the agreement is very good. The growth parameter gives an indication of how
wide or narrow the breather will be, with larger growth parameters
corresponding to narrower solutions. The prediction from the linear theory is
given by the real part of the Floquet exponent corresponding to mode $m_{0}$,
and the approximation from the perturbation analysis is (59). A comparison of
these two quantities are shown as the red markers and lines, respectively of
Fig. 6(c). The trends in Fig. 6(c) demonstrate that the $q$-gap breathers
becomes larger in amplitude and more narrow as the wavenumber goes deeper into
the gap.

In Fig. 6(c) the $q$-gap breathers are generated up until $\Delta q=0.042$
(which corresponds to $\varepsilon=0.25$). For $\Delta q=0$, the width of the
wavenumber bandgap is $q_{r}-q_{\ell}\approx 0.28$ (where $q_{r}$, $q_{\ell}$
are the right and left edges of the bandgap, respectively). Thus, the branch
of solutions shown Fig. 6(c) extends to roughly 15% of the width of the
bandgap. For $\Delta q>0.042$ we did not observe a coherent temporal
localization. Indeed, for all the breathers observed numerically, the
localization is obtained for a finite interval of time. For longer time
simulations the amplitude of the breather can grow again (leading to a
repeated appearance of breathers), see Fig. 7(a). Similar observations have
been made for k-gap solitons in photonic systems [34]. While Theorem 1.1
guarantees that a temporally localized structures exists over a finite
temporal interval, there is no statement about the dynamics beyond this
interval. The numerical simulations suggest the tail of the breather can
experience repeated growth. We observed as $\Delta q$ becomes larger, the time
between consecutive peaks of the pulses becomes smaller. In other words, the
emergence of the “second” breather occurs faster as $\Delta q$ becomes larger.
Thus, for sufficiently large $\Delta q$ the structure is not temporally
localized since the second breather emerges “too soon”, see Fig. 7(b). This is
the reason why we only show $\Delta q\leq 0.042$ in Fig. 6(c).

(a) (b) (c)

Figure 8: (a) Spectral bands (blue curves) of the Schrödinger operator as a
function of $\ell$, the imaginary part of the Floquet exponent. The gray dots
shown the corresponding values in the finite lattice with $N=10$. The
parameter values are $c=0$, $T=1/0.37$ and $K_{2}=\underline{m}=1$. The
critical modulation amplitude parameters are $k_{a}^{0}=0.5$ and
$k_{b}^{0}=0.7487$. The critical mode number is $m_{0}=4$, which lies at the
bottom of the second spectral band, namely $\lambda(\ell(q_{4}))=\mu_{2}$ (see
larger black marker). The eigenvalues $\mu_{1}$ and $\mu_{2}$ define the edges
of the band gap, shown as the gray shaded region. (b) Example $q$-gap breather
bifurcating from $\mu_{2}$ with nonlinear coefficients $K_{3}=0$ and
$K_{4}=0.8$. The initial value is $u_{n}(0)=10^{-4}\sin(q_{4}n)$ and
$\varepsilon=0.1$. The displacement of the 6th particle is shown as a function
of time. The modulation amplitude is $k_{a}=k_{a}^{0}+\varepsilon^{2}=0.06$,
$k_{b}=k_{b}^{0}+\varepsilon^{2}=0.7587$. The dashed line shows an
approximation of the envelope given by Eq. (70) with $\chi=1.7707$,
$\lambda^{\prime\prime}(\ell_{0})=26.3821$ and $\varepsilon=0.1$. For this
value of $\varepsilon$, the distance of the wavenumber to the edge of the gap
is $\Delta q=q_{4}-q_{r}=-0.005$. (c) Plot of the amplitude of the breather
($\mathrm{max}_{t}u_{6}(t)$) for the numerical simulation (blue dots) and
prediction based on Eq. (70) (blue line) as a function of $\Delta q$. The real
part of the Floquet exponent corresponding to $q_{4}$ (solid red squares) and
asymptotic approximation
$\varepsilon\sqrt{\delta\lambda^{\prime\prime}(\ell_{0})/2}$ (red line) are
also shown, which indicate the growth (decay) rate of the q-gap breather.

Next, we consider an example where the breathers bifurcate from $\mu_{2}$. The
spectral bands corresponding to $c=0$, $T=1/0.37$, $K_{2}=\underline{m}=1$,
$k_{a}^{0}=0.5$ and $k_{b}^{0}=0.7487$ are shown in Fig. 8(a). For these
parameter values the mode $m_{0}=4$ lies at the bottom of the second spectral
band, namely $\lambda(q_{4})=\mu_{2}$. Notice that the concavity of the second
spectral band is opposite of the first band, namely
$\lambda_{2}^{\prime\prime}(\ell_{0})>0$. Thus, in order to satisfy the
(Coeff) condition, we require $\chi>0$. If $K_{3}=0$ this implies that
$K_{4}>0$ and the sign parameter is now $\delta=-1.$ Thus, for the next
numerical simulations, we fix $K_{3}=0$ and $K_{4}=0.8.$ The approximation
given in Eq. (70) is identical in this case, but we replace
$\lambda_{1}^{\prime\prime}(\ell)$ with $\lambda_{2}^{\prime\prime}(\ell)$,
and likewise for the underlying Bloch modes (where $f_{1}$ should be replaced
by $f_{2}$, etc.).

Figure 8(b) shows an example of the generalized $q$-gap breather with
$\varepsilon=0.1$, with corresponding envelope prediction given by Eq. (70).
Qualitatively, the results are similar to the example shown in Fig. 6(a).
Figure 8(c) shows the dependence of the breather amplitude and growth rate on
the parameter $\Delta q=q_{r}-q_{m_{0}}$. Note, since the breather is
bifurcating from the right edge of the wavenumber bandgap, the quantity
$\Delta q$ will be negative. The breather amplitude grows like
$\mathcal{O}(\sqrt{|\Delta q|})$.

###  8.2 Examples with $c=0$ and $K_{3}\neq 0$

Here we will consider $K_{3}\neq 0$. In particular, we will chose values of
the nonlinear coefficients in (4) that correspond to the modulated magnetic
lattice described in Sec. 1.1 so that the results obtained here are directly
relevant for the experimental set-up described in [7]. In the re-scaled
variables the interaction coefficients are $K_{2}=K_{3}=1$ and $K_{4}=0.8$.
Since $K_{3}\neq 0$ the sign of $\chi$ must be computed directly to see if the
relevant eigenvalue to bifurcate from is $\mu_{1}$ or $\mu_{2}$. $\chi$ will
depend on the function $h_{2}(t)$, which we can obtain by solving Eq. (64). It
will be convenient to estimate $h_{2}(t)$ numerically under the constraint
that $h_{2}(t)$ is $T$ periodic, which we achieve using a shooting method. In
particular, we apply Newton iterations on the map
$F(h^{0})=h_{2}(0;h^{0})-h_{2}(T;h^{0})$ where $h_{2}(t;h^{0})$ is the
solution of Eq. (64) with initial condition
$h^{0}=(h_{2}(0),\dot{h}_{2}(0))^{T}$. The Jacobian of the map $F$ is simply
$I-V(T)$, where $I$ is the 2x2 identity matrix and $V(T)$ is the solution to
the variational equation $\dot{V}=\frac{df}{dh}V$ with initial value $V(0)=I$
where $\frac{df}{dh}$ is the Jacobian corresponding to Eq. (64) [20].

(a) (b) (c)

Figure 9: (a) A generalized $q$-gap breather bifurcating from $\mu_{2}$ with
nonlinear coefficients corresponding to a magnetic lattice, namely $K_{3}=1$
and $K_{4}=0.8$. The other parameter values and spectral picture is identical
to Fig. 8(a). The initial value is $u_{n}(0)=10^{-4}\sin(q_{4}n)$ and
$\varepsilon=0.1$. The displacement of the 6th particle is shown as a function
of time. The modulation amplitude is $k_{a}=k_{a}^{0}+\varepsilon^{2}=0.06$,
$k_{b}=k_{b}^{0}+\varepsilon^{2}=0.7587$. The dashed line shows an
approximation of the envelope given by Eq. (70) with $\chi=1.997$,
$\lambda^{\prime\prime}(\ell_{0})=26.3821$ and $\varepsilon=0.1$. For this
value of $\varepsilon$, the distance of the wavenumber to the edge of the gap
is $\Delta q=q_{4}-q_{r}=-0.005$. (b) Same as panel (a) with
$\varepsilon=0.22$. (c) Plot of the amplitude of the breather
($\mathrm{max}_{t}u_{6}(t)$) for the numerical simulation (blue dots) and
prediction based on Eq. (70) (blue line) as a function of $\Delta q$. The real
part of the Floquet exponent corresponding to $q_{4}$ (solid red squares) and
asymptotic approximation (59) (red line) are also shown, which indicate the
growth (decay) rate of the breather.

In this example, we consider the spectral situation as shown in Fig. 8(a),
such that $m_{0}=4$ is the critical mode bifurcating from $\mu_{2}$. In this
case, $\lambda^{\prime\prime}(\ell_{0})>0$, so we chose $\delta=-1$ and we
must have $\chi>0$ to satisfy (Coeff). Upon computing $h_{2}(t)$ with the
shooting method and substituting into Eq. (68) with $K_{3}=1$ and $K_{4}=0.8$
we find $\chi=1.997>0$, as desired.

Figure 9(a) shows a numerical simulation of the lattice with initial
displacement $u_{n}(0)=10^{-4}\sin(q_{4}n)$ and $\varepsilon=0.1$, and panel
(b) shows a simulation with $\varepsilon=0.22$. By comparing panels (a) and
(b), we see once again that the $q$-gap breather becomes more narrow and
larger in amplitude as $\Delta q$ (and thus $\varepsilon$) increases in
magnitude. What is also apparent, especially in panel (b), is that the
numerical simulation is asymmetric, namely, the maximum is not simply the
minimum reflected about the $u=0$ line. Evidently, the asymmetric nature of
the FPUT potential with $K_{3}\neq 0$ is manifested through a lack of
reflection symmetry in the $q$-gap breather profile. Asymmetric breathing
profiles are well known in space-periodic FPUT systems with quadratic
nonlinearities [21]. The approximation given by Eq. (70) remains symmetric,
however, and thus one would expect the approximation not to do as well as in
the $K_{3}=0$ case. Indeed, inspection of Fig. 8(c) confirms this, where the
difference in amplitude between simulation and theory is larger than in the
$K_{3}=0$ case (see e.g., Fig. 8(c)). Nonetheless, the asymptotic behavior as
$\Delta q\rightarrow 0$ is correct, and in particular the breather amplitude
grows like $\mathcal{O}(\sqrt{|\Delta q|})$.

###  8.3 Examples with $c\neq 0$ and $K_{3}\neq 0$

In our final example, we include the effect of damping and select
$\tilde{c}=0.1$. By definition, $c=\varepsilon\tilde{c}$, so the critical
parameter set (when $\varepsilon=0$) will have $c=0$, like before. Thus, we
consider once again the parameter set that corresponds to Fig. 9(a). However,
the numerical solutions and asymptotic approximations will have non-zero
damping affect for $\varepsilon>0$. Since the bifurcation scenario is the same
as in Fig. 9(a) the initial condition for simulations will be of the same
form, namely $u_{n}(0)=10^{-4}\sin(q_{4}n)$. An example lattice simulation
with $\varepsilon=0.1$ is shown in Fig. 10(a) where the underlying damping
constant is $c=\varepsilon\tilde{c}=0.01$. The solution experiences an initial
growth, with growth rate given by the real part of the $m_{0}=4$ Floquet
exponent, but rather then decaying to a near zero amplitude, like in all the
previous examples with $c=0$, the solution approaches steady periodic motion
with period $2T$. A longer time evolution of the same solution is shown in
Fig. 11. In particular, panel (b) shows the dynamics are essentially periodic
for $t$ sufficiently larger.

In terms of the Poincaré map $F_{j}=U(2Tj)$, the trivial solution $U(t)=0$ is
clearly a fixed point. The $2T$-periodic solution that is approached in the
dynamic simulation is another fixed point. Thus, the solution shown in Fig.
10(a) is a transition front, since it connects two different fixed points.

Another example of the transition front for a larger value of $\varepsilon$ is
shown in Fig. 10(b). Despite the fact that the structure is not temporally
localized, the initial dynamics still resemble the “left” side of the $q$-gap
breather. Indeed, the homoclinic approximation from Eq. (70) is quite close to
the initial front dynamics (see the solid gray line of Fig. 10(a)). For this
reason, it is still reasonable to measure the amplitude of the front in the
same way we measured the amplitude for the $q$-gap breathers. A plot of the
front amplitude and real part of the $m_{0}=4$ Floquet exponent is shown in
10(c). The amplitude trend is similar to the non-damped case, but the
magnitude of the amplitude is smaller, as expected (compare panel (c) of Figs.
9 and 10).

(a) (b) (c)

Figure 10: (a) A transition front bifurcating from $\mu_{2}$ with nonlinear
coefficients corresponding to a damped magnetic lattice, namely $K_{3}=1$,
$K_{4}=0.8$ and $\tilde{c}=0.1$. The other parameter values and spectral
picture is identical to Fig. 8(a). The initial value is
$u_{n}(0)=10^{-4}\sin(q_{4}n)$ and $\varepsilon=0.1$. The displacement of the
6th particle is shown as a function of time. The modulation amplitude is
$k_{a}=k_{a}^{0}+\varepsilon^{2}=0.06$,
$k_{b}=k_{b}^{0}+\varepsilon^{2}=0.7587$ and the damping constant is
$c=\varepsilon\tilde{c}=0.01$. The dashed line shows an approximation of the
envelope given by Eq. (60) with $\chi=1.997$,
$\lambda^{\prime\prime}(\ell_{0})=26.3821$ and $\varepsilon=0.1$. For this
value of $\varepsilon$, the distance of the wavenumber to the edge of the gap
is $\Delta q=q_{4}-q_{r}=-0.005$. (b) Same as panel (a) with
$\varepsilon=0.22$. (c) Plot of the amplitude of the front
($\mathrm{max}_{t}u_{6}(t)$) for the numerical simulation (blue dots) and
prediction based on Eq. (60) (blue line) as a function of $\Delta q$. The real
part of the Floquet exponent corresponding to $q_{4}$ (solid red squares) and
asymptotic approximation $r_{0}^{+}\varepsilon$ (red line) are also shown,
which indicates the initial growth rate.

(a) (b)

Figure 11:  (a) Same as Fig. 10(a) but over a longer time interval. The dashed
line shows an approximation of the envelope given by Eq. (60), which accounts
for damping. The solid gray line shows the homoclinic approximation (i.e.,
with no damping) of the envelope given by Eq. (70). (b) A zoom of panel (a)
for large values of $t$. Here it can be seen the dynamics are very close to
periodic.

The amplitude equation Eq. (66) can be used to approximate the front dynamics.
However, the equation does not yield an explicit solution in the presence of
damping and so we will employ a qualitative and numerical analysis of Eq.
(66). A straightforward phase plane analysis shows that the trivial fixed
point $(A,A^{\prime})=(0,0)$ is a saddle with corresponding eigenvalues
$r_{0}^{\pm}=(-s_{3}\pm\sqrt{s_{3}^{2}+4s_{1}})/2$ and the fixed points
$(A,A^{\prime})=(\pm\sqrt{s_{1}/s_{2}},0)$ are spiral-sinks with corresponding
eigenvalues $r_{1}^{\pm}=(-s_{3}\pm i\sqrt{8s_{1}-s_{3}^{2}})/2$ where the
$s_{j}$ are the coefficients of Eq. (66), namely,

$s_{1}=\left(\frac{\tilde{c}^{2}}{4\underline{m}}-\delta\right)\frac{2}{\lambda^{\prime\prime}(\ell_{0})}+\frac{\tilde{c}^{2}}{4\underline{m}^{2}}>0,\qquad
s_{2}=\frac{2\chi}{\lambda^{\prime\prime}(\ell_{0})}>0,\qquad
s_{3}=\frac{\tilde{c}}{\underline{m}}>0.$

In the phase plane, there is a heteroclinic orbit that leaves the trivial
fixed point along the unstable eigenvector $(1,r_{+}^{0})^{T}$ and approaches
the $(\sqrt{s_{1}/s_{2}},0)$ fixed point. There is another heteroclinic orbit
that leaves the trivial fixed point along the unstable eigenvector
$(-1,-r_{+}^{0})^{T}$ and approaches the $(-\sqrt{s_{1}/s_{2}},0)$ fixed
point, see Fig. 5. To approximate the heteroclinic orbit, we numerically solve
Eq. (66) with initial condition
$A(0)=10^{-4},A^{\prime}(0)=10^{-4}\,r_{+}^{0}$. The resulting solution
$A(\tau)$ is then used in Eq. (60) to generate the approximation of the
lattice dynamics.

Examples are shown in Fig. 10(a,b), where the envelopes are shown for
$\varepsilon=0.1$ and $\varepsilon=0.22$, respectively. Once again, the
envelope dynamics are captured well by Eq. (60), especially for small
$\varepsilon$. The periodic oscillation of the envelope can be approximated by
the imaginary part of the eigenvalue associated to the non-trivial fixed
point, namely $\tau_{\rm env}=2\pi/(\sqrt{8s_{1}-s_{3}^{2}}/2)$. In terms of
the original lattice variables this translates to $\tau_{\rm
env}/\varepsilon$. For the example shown in Fig. 11(a) with $\varepsilon=0.1$,
the average peak-to-peak time of the envelope is $165.3$ time units, where as
$\tau_{\rm env}/\varepsilon=159.3$, which is quite close.

The front amplitude as a function of $\Delta q$ is shown as the solid blue
line in panel (c), and an approximation of the initial growth rate
$r_{+}^{0}\varepsilon$ is shown as the red line. Once again the asymptotic
behavior as $\Delta q\rightarrow 0$ is correct. Despite the presence of
damping, the front amplitude grows like $\mathcal{O}(\sqrt{|\Delta q|})$.

##  9 Conclusions

Generalized $q$-gap breathers are coherent structures that are localized in
time, periodic in space, and have wavenumber in a $q$-gap. They are the
natural counterpart of the discrete breathers of spatially periodic lattices,
which themselves are of fundamental importance in a diverse range of fields.

In the absence of damping, we proved rigorously the existence of generalized
$q$-gap breathers in a time-periodic FPUT lattice using normal form theory. In
particular, we proved the existence of oscillating homoclinic solutions over a
finite time interval with tails that can be made arbitrarily small, but
finite. These solutions bifurcate from one edge of the $q$-gap. Which of the
edges is determined by the nonlinear coefficients $K_{3},K_{4}$ and the
concavity of the spectral band from which the solution bifurcates. The
amplitude of the $q$-gap breather grows like $\mathcal{O}(\sqrt{\Delta q})$,
where $\Delta q$ is the distance of the underlying wavenumber to the band
edge. This result makes rigorous the numerical and experimental observations
of such $q$-gap breathers in [7]. We also provided a tractable analytical
approximation of such solutions using a multiple-scale analysis and
corroborated results with direct numerical simulations.

In the presence of damping we proved the existence of solutions that connect
the zero state to a time-periodic one, which we called the transition fronts.
The multiple-scale analysis also provided an accurate description of the front
solutions, although the underlying amplitude equation needed to be solved
numerically. The initial stage of the front dynamics were well described by
the undamped $q$-gap breather approximations.

Generalized $q$-gap breathers and transition fronts represent new types of
nonlinear wave structures. This work provided the first rigorous results in
their study, complementing earlier experimental and numerical work.
Nonetheless, there are still many open questions regarding $q$-gap breathers
and transition fronts. This includes the possible existence of genuine $q$-gap
breathers (i.e., with both tails decaying to zero), the numerically exact
computation of $q$-gap breathers (i.e., numerical roots of the appropriate map
up to a user-prescribed tolerance) and the exploration of such structures in
higher spatial dimensions or in settings beyond the FPUT realm. Indeed, any
system that is already described by a nonlinear wave equation that could be
adapted to be time-varying (in order to induce a $q$-gap) would be a candidate
for the implementation of $q$-gap breathers. This suggests that $q$-gap
breathers’ relevance, and hence the results of this work, could extend to a
wide range of fields.

## References

  * [1] N. Akhmediev, A. Ankiewicz, and J. M. Soto-Crespo, Rogue waves and rational solutions of the nonlinear Schrödinger equation, Phys. Rev. E, 80 (2009), p. 026601.
  * [2] A. Ankiewicz, N. Akhmediev, and J. M. Soto-Crespo, Discrete rogue waves of the Ablowitz-Ladik and Hirota equations, Phys. Rev. E, 82 (2010), p. 026602.
  * [3] G. P. Berman and F. M. Izrailev, The fermi-pasta-ulam problem: Fifty years of progress, Chaos, 15 (2005).
  * [4] M. Butt, S. Khonina, and N. Kazanskiy, Recent advances in photonic crystal optical devices: A review, Optics and Laser Technology, 142 (2021), p. 107265.
  * [5] L. S. Cao, D. X. Qi, R. W. Peng, M. Wang, and P. Schmelcher, Phononic frequency combs through nonlinear resonances, Phys. Rev. Lett., 112 (2014), p. 075505.
  * [6] M. Centurion, M. A. Porter, Y. Pu, P. G. Kevrekidis, D. J. Frantzeskakis, and D. Psaltis, Modulational instability in nonlinearity-managed optical media, Phys. Rev. A, 75 (2007), p. 063804.
  * [7] C. Chong, B. Kim, E. Wallace, and C. Daraio, Modulation instability and wavenumber bandgap breathers in a time layered phononic lattice, Phys. Rev. Research, 6 (2024), p. 023045.
  * [8] T. Dauxois, Fermi, Pasta, Ulam, and a mysterious lady, Phys. Today, 61(1) (2008), p. 55.
  * [9] S. V. Dmitriev, E. A. Korznikova, Y. A. Baimova, and M. G. Velarde, Discrete breathers in crystals, Physics-Uspekhi, 59 (2016), pp. 446–461.
  * [10] T. Dohnal, D. E. Pelinovsky, and G. Schneider, Travelling modulating pulse solutions with small tails for a nonlinear wave equation in periodic media, Nonlinearity, 37 (2024), p. 055005.
  * [11] E. Fermi, J. Pasta, S. Ulam, and M. Tsingou, Studies of nonlinear problems i, 1955.
  * [12] S. Flach and A. Gorbach, Discrete breathers: advances in theory and applications, Phys. Rep., 467 (2008), p. 1.
  * [13] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, $q$-Breathers and the Fermi-Pasta-Ulam Problem, Phys. Rev. Lett., 95 (2005), p. 064102.
  * [14] R. Fukuizumi and G. Schneider, Interchanging space and time in nonlinear optics modeling and dispersion management models, J. Nonlinear Sci., 32 (2022), pp. 29, 39.
  * [15] G. Gallavotti, The Fermi–Pasta–Ulam Problem: A Status Report, Springer-Verlag, Berlin, Germany, 2008.
  * [16] A. Ganesan, C. Do, and A. Seshia, Phononic frequency comb via intrinsic three-wave mixing, Phys. Rev. Lett., 118 (2017), p. 033903.
  * [17] M. D. Groves and G. Schneider, Modulating pulse solutions for quasilinear wave equations, Journal of Differential Equations, 219 (2005), pp. 221–258.
  * [18] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, vol. 42 of Appl. Math. Sci., Springer, Cham, 1983.
  * [19] W. Hartmann, Acoustic Signal Processing, Springer New York, New York, NY, 2007, pp. 503–530.
  * [20] M. Hirsch, S. Smale, and R. Devaney, Differential Equations, Dynamical Systems, and an Introduction to Chaos, Pure and Applied Mathematics - Academic Press, Elsevier Science, 2004.
  * [21] G. Huang and B. Hu, Asymmetric gap soliton modes in diatomic lattices with cubic and quartic nonlinearity, Phys. Rev. B, 57 (1998), p. 5746.
  * [22] G. Iooss and M. Adelmeyer, Topics in bifurcation theory and applications., vol. 3 of Adv. Ser. Nonlinear Dyn., Singapore: World Scientific, 2nd ed. ed., 1998.
  * [23] P. G. Kevrekidis, Non-linear waves in lattices: Past, present, future, IMA J. Appl. Math., 76 (2011), pp. 389–423.
  * [24] B. L. Kim, C. Chong, S. Hajarolasvadi, Y. Wang, and C. Daraio, Dynamics of time-modulated, nonlinear phononic lattices, Phys. Rev. E, 107 (2023), p. 034211.
  * [25] B. L. Kim, C. Daraio, C. Chong, S. Hajarolasvadi, and Y. Wang, Dynamics of time-modulated, nonlinear phononic lattices, arxiv:, (2022), p. 2209.06511.
  * [26] I. Kovacic, R. Rand, and S. M. Sah, Mathieu’s equation and its generalizations: Overview of stability charts and their features, Applied Mechanics Reviews, 70 (2018).
  * [27] A. B. Kozyrev, H. Kim, and D. W. van der Weide, Parametric amplification in left-handed transmission line media, Applied Physics Letters, 88 (2006), p. 264101.
  * [28] G. Lee, D. Lee, J. Park, Y. Jang, M. Kim, and J. Rho, Piezoelectric energy harvesting using mechanical metamaterials and phononic crystals, Commun Phys, 5 (2022), p. 94.
  * [29] L. Liu, W.-R. Sun, B. A. Malomed, and P. G. Kevrekidis, Time-localized dark modes generated by zero-wave-number-gain modulational instability, Phys. Rev. A, 108 (2023), p. 033504.
  * [30] E. Lombardi, Orbits homoclinic to exponentially small periodic orbits for a class of reversible systems. Application to water waves, Arch. Ration. Mech. Anal., 137 (1997), pp. 227–304.
  * [31] J. Marconi, E. Riva, M. Di Ronco, G. Cazzulani, F. Braghin, and M. Ruzzene, Experimental Observation of Nonreciprocal Band Gaps in a Space-Time-Modulated Beam Using a Shunted Piezoelectric Array, Physical Review Applied, 13 (2020), p. 031001.
  * [32] S. Mukherjee and M. C. Rechtsman, Observation of unidirectional solitonlike edge states in nonlinear floquet topological insulators, Phys. Rev. X, 11 (2021), p. 041057.
  * [33] H. Nassar, B. Yousefzadeh, R. Fleury, M. Ruzzene, A. Alù, C. Daraio, A. N. Norris, G. Huang, and M. R. Haberman, Nonreciprocity in acoustic and elastic materials, Nature Reviews Materials, 5 (2020), p. 667–685.
  * [34] Y. Pan, M.-I. Cohen, and M. Segev, Superluminal $k$-gap solitons in nonlinear photonic time crystals, Phys. Rev. Lett., 130 (2023), p. 233801.
  * [35] D. A. Powell, I. V. Shadrivov, and Y. S. Kivshar, Multistability in nonlinear left-handed transmission lines, Applied Physics Letters, 92 (2008), p. 264104.
  * [36] , Asymmetric parametric amplification in nonlinear left-handed transmission lines, Applied Physics Letters, 94 (2009), p. 084105.
  * [37] Z. Rapti, G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and B. A. Malomed, Modulational instability in bose–einstein condensates under feshbach resonance management, Physica Scripta, 2004 (2004), p. 27.
  * [38] J. R. Reyes-Ayona and P. Halevi, Observation of genuine wave vector ($k$ or $\beta$) gap in a dynamic transmission line and temporal photonic crystals, Applied Physics Letters, 107 (2015), p. 074101.
  * [39] M. Safaei, H. A. Sodano, and S. R. Anton, A review of energy harvesting using piezoelectric materials: state-of-the-art a decade later (2008–2018), Smart Materials and Structures, 28 (2019), p. 113001.
  * [40] J. A. Sanders, F. Verhulst, and J. Murdock, Averaging methods in nonlinear dynamical systems, vol. 59 of Appl. Math. Sci., New York, NY: Springer, 2nd ed. ed., 2007.
  * [41] M. Soljačić, M. Ibanescu, S. G. Johnson, Y. Fink, and J. D. Joannopoulos, Optimal bistable switching in nonlinear photonic crystals, Physical Review E, 66 (2002), p. 055601.
  * [42] G. Trainiti, Y. Xia, J. Marconi, G. Cazzulani, A. Erturk, and M. Ruzzene, Time-Periodic Stiffness Modulation in Elastic Metamaterials for Selective Wave Filtering: Theory and Experiment, Physical Review Letters, 122 (2019), p. 124301.
  * [43] A. Vainchtein, Solitary waves in fpu-type lattices, Physica D, 434 (2022), p. 133252.
  * [44] F. Y. Wang, G. X. Li, H. L. Tam, K. W. Cheah, and S. N. Zhu, Optical bistability and multistability in one-dimensional periodic metal-dielectric photonic crystal, Applied Physics Letters, 92 (2008), p. 211109.
  * [45] C.-P. Wen, W. Liu, and J.-W. Wu, Tunable terahertz optical bistability and multistability in photonic metamaterial multilayers containing nonlinear dielectric slab and graphene sheet, Applied Physics A, 126 (2020), p. 426.
  * [46] M. Yu, J. K. Jang, Y. Okawachi, A. G. Griffith, K. Luke, S. A. Miller, X. Ji, M. Lipson, and A. L. Gaeta, Breather soliton dynamics in microresonators, Nature Communication, 8 (2017), p. 14569.