# Travelling waves in nonlinear magneto-inductive lattices

M. Agaoglou, M. Fečkan, M. Pospíšil, V.M. Rothos, H. Susanto

► To cite this version:

M. Agaoglou, M. Fečkan, M. Pospíšil, V.M. Rothos, H. Susanto. Travelling waves in nonlinear magneto-inductive lattices. Journal of Differential Equations, 2016, 260 (2), pp.1717-1746. 10.1016/j.jde.2015.09.043. hal-01343794

HAL Id: hal-01343794

https://hal.science/hal-01343794

Submitted on 12 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Travelling waves in nonlinear magneto-inductive lattices

M. Agaoglou<sup>a</sup>, M. Fečkan<sup>b,c</sup>, M. Pospíšil<sup>c</sup>, V.M. Rothos<sup>a,*</sup>, H. Susanto<sup>d</sup>

<sup>a</sup> Lab of Nonlinear Mathematics and Department of Mechanical Engineering, Faculty of Engineering, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece

<sup>b</sup> Department of Mathematical Analysis and Numerical Mathematics, Comenius University in Bratislava, Mlynská dolina, 842 48 Bratislava, Slovakia

<sup>c</sup> Mathematical Institute of Slovak Academy of Sciences, Štefánikova 49, 814 73 Bratislava, Slovakia

<sup>d</sup> Department of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, United Kingdom

* Corresponding author. Michal.Feckan@fmph.uniba.sk, Michal.Pospisil@mat.savba.sk, rothos@auth.gr, husanto@essex.ac.uk.

We consider a lattice equation modelling one-dimensional metamaterials formed by a discrete array of nonlinear resonators. We focus on periodic travelling waves due to the presence of a periodic force. The existence and uniqueness results of periodic travelling waves of the system are presented. Our analytical results are found to be in good agreement with direct numerical computations.

Keywords: Travelling wave; Lattice wave; Forced magneto-inductive lattice

## 1. Introduction

In this work, we consider [26]

$$ \frac{d^2}{dt^2}(u_n - \lambda u_{n-1} - \lambda u_{n+1}) + \gamma \frac{d}{dt}u_n + u_n + \frac{d^2}{dt^2}(u_n^2 - \lambda u_{n-1}^2 - \lambda u_{n+1}^2) + \gamma \frac{d}{dt}u_n^2 - h(\omega t + pn) = 0 \quad (1) $$
where $\gamma \ge 0$, $\lambda, \omega > 0$, $p \ne 0$ are parameters and $h \in C(\mathbb{R}, \mathbb{R})$ is such that

$$h(z) = \sum_{k \in \mathbb{Z}} h_k e^{k_1 z}, \quad \sum_{k \in \mathbb{Z}} |h_k| < \infty.$$

The equation models the dynamics of electromagnetic waves in the so-called magneto-inductive metamaterials.

Metamaterials are artificial materials that are engineered to have properties that may not be found in nature [5]. The (structural rather than chemical) engineering is achieved by composing periodic inhomogeneities to create desirable effective behaviour. The invention ignited a new paradigm in electromagnetism, including cloaking devices [20] (see also [17] for a recent comprehensive review of electromagnetic manipulation enabled by metamaterials). Magnetic metamaterials are non-magnetic materials exhibiting magnetic properties in the Terahertz and optical frequencies. They were predicted theoretically in [12,28] and demonstrated experimentally in, e.g., [7,11,14,27]. When the materials composed of magnetic metamaterials, e.g., split-ring resonators (SRRs), are magnetically weakly coupled through their mutual inductances, one obtains magneto-inductive metamaterials [9,21–23]. Magneto-inductive metamaterials consisting of periodic arrays of SRRs have been made in one-, two- and three dimensions [19]. Both the dimensions of SRRs and their inter-space distance are small relative to the free space wavelength and, thus, the dynamics of electromagnetic fields in such settings are governed by the quasi magnetostatic approximation [13].

Model equations for the propagation of nonlinear electromagnetic waves in metamaterials can be grouped into two classes. The first approach assumes that effective metamaterials are homogeneous media with specific physical properties resulting in partial differential equations, such as coupled short-pulse equation [24] and higher-order nonlinear Schrödinger equations [25]. In the second class, metamaterials are modelled by arrays of coupled oscillators, i.e. lattice equations, such as a nonlinear Klein–Gordon equation [18] and coupled Klein–Gordon equations [15,16]. The governing equation (1) falls into the second approach with $\gamma$ representing the loss coefficient, $\lambda$ is the coupling parameter, $h$ is an external forcing that is periodic in time and varying in the spatial direction.

In the lattice equation (1) the nonlinearity appears in the coupling terms and in the damping. This is due to the assumption of the nonlinearity of the capacitance of the split-ring resonators that compose the magneto-inductive materials [26]. Note that the nonlinearity is different from our previous work [1,6], where it is in the onsite potential. The present work is to study the effect of such nonlinearity. Here, we investigate the existence of travelling periodic solutions of system (1) due to the periodic forcing, and the bifurcation of such solutions with small $\gamma$, $\lambda$ and $h$. We also study the modulational stability of the periodic solutions by computing Floquet multipliers of the linearized systems.

Note that in our governing equation (1) the nonlinearity in the coupling terms between the sites is akin to that in the Fermi–Pasta–Ulam lattices [10]. However, there is a significant difference in fact that the coupling in our governing equation is also the derivative term. The bifurcation structures of periodic solutions in general FPU lattices forced by periodic drive were studied in [8]. It is imperative to study the effect of the present nonlinear couplings to periodic solutions caused by the same periodic drive. Using Lyapunov–Schmidt reduction, we derive the asymptotic expressions of the bifurcating solutions.

The present paper is organized as follows. In Section 2 we are looking for a periodic travelling wave when its amplitude is limited by a function of the magnitude of the forcing. In the section,
we prove the existence of such waves rigorously. The resonance condition for the parameter values is also derived. In Section 3, we study the resonance region. Using the Lyapunov–Schmidt reduction, we show that the periodic solutions persist. Several explicit examples of the general results obtained in the previous sections are presented in Section 4. In Section 5 we solve the governing equations numerically. Comparisons with the analytical results are presented where we obtain good agreement.

## 2. Existence of small periodic solutions

Putting $u_n(t) = U(\omega t + pn)$, $z = \omega t + pn$ in (1), we obtain

$$ \omega^2(U(z) - \lambda U(z-p) - \lambda U(z+p))'' + \gamma \omega U'(z) + U(z) + \omega^2(U^2(z) - \lambda U^2(z-p) - \lambda U^2(z+p))'' + \gamma \omega (U^2(z))' - h(z) = 0. \quad (2) $$

We take Banach spaces

$$ X := \left\{ U \in C(\mathbb{R}, \mathbb{R}) \mid U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k \iota z}, \sum_{k \in \mathbb{Z}} |c_k| < \infty \right\}, $$

$$ Y := \left\{ U \in C^1(\mathbb{R}, \mathbb{R}) \mid U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k \iota z}, \sum_{k \in \mathbb{Z} \setminus \{0\}} |k| |c_k| < \infty \right\}, $$

$$ Z := \left\{ U \in C^2(\mathbb{R}, \mathbb{R}) \mid U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k \iota z}, \sum_{k \in \mathbb{Z} \setminus \{0\}} k^2 |c_k| < \infty \right\} $$

with the norms

$$ \|U\| := \sum_{k \in \mathbb{Z}} |c_k|, \quad \|U\|_1 := |c_0| + \sum_{k \in \mathbb{Z} \setminus \{0\}} |k| |c_k|, \quad \|U\|_2 := |c_0| + \sum_{k \in \mathbb{Z} \setminus \{0\}} k^2 |c_k|, $$

respectively. It is easy to verify that $Z \circled{\circ} Y \circled{\circ} X$ are compact embeddings and $\|U\| \le \|U\|_1 \forall U \in Y$, $\|U\|_1 \le \|U\|_2 \forall U \in Z$.

The following lemma is clear.

**Lemma 1.** If $U_1, U_2 \in X$ then $U_1 U_2 \in X$ and $\|U_1 U_2\| \le \|U_1\| \|U_2\|$.

By setting

$$ \mathcal{K}U := \omega^2 (U''(z) - \lambda U''(z-p) - \lambda U''(z+p)) + \gamma \omega U'(z) + U(z), $$

$$ \mathcal{F}(U, h) := -\omega^2 \left( (U^2(z))'' - \lambda (U^2(z-p))'' - \lambda (U^2(z+p))'' \right) - \gamma \omega (U^2(z))' + h(z), $$

equation (2) has the form

$$ \mathcal{K}U = \mathcal{F}(U, h). $$
We have the next result.

**Lemma 2.** Function $\mathcal{F}: Z \times X \to X$ fulfils

$$ \|\mathcal{F}(U, h)\| \le 2\omega(2\omega + 4\lambda\omega + \gamma)\|U\|_2^2 + \|h\|, \quad (3) $$

$$ \|\mathcal{F}(U_1, h) - \mathcal{F}(U_2, h)\| \le 2\omega(2\omega + 4\lambda\omega + \gamma)\|U_1 - U_2\|_2(\|U_1\|_2 + \|U_2\|_2), \quad (4) $$

$$ \|\mathcal{F}(U, h_1) - \mathcal{F}(U, h_2)\| \le \|h_1 - h_2\| \quad (5) $$

for any $U, U_1, U_2 \in Z$ and $h, h_1, h_2 \in X$.

**Proof.** Property (5) is obvious. Next, we use $\|U'\| = \|U\|_1 \le \|U\|_2$, $\|U''\| = \|U\|_2$ with Lemma 1 to derive

$$ \|(U^2(z+c))''\| = 2\|(U'(z+c))^2 + U(z+c)U''(z+c)\| \le 4\|U(z+c)\|_2^2 $$

for $c = -p, 0, p$. Moreover,

$$ \|U(z \pm p)\| = \sum_{k \in \mathbb{Z}} |c_k| \|e^{k_1 p}\| = \sum_{k \in \mathbb{Z}} |c_k| = \|U\| $$

for $U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k_1 z}$. Similarly we derive

$$ \|(U^2(z))'\| = 2\|U(z)U'(z)\| \le 2\|U\|_2^2. $$

Hence, property (3) is obtained. Analogically, we get

$$ \begin{aligned} \|(U_1^2(z+c) - U_2^2(z+c))''\| &\le \|(U_1(z+c) - U_2(z+c))''(U_1(z+c) + U_2(z+c))\| \\ &\quad + 2\|(U_1(z+c) - U_2(z+c))'(U_1(z+c) + U_2(z+c))'\| \\ &\quad + \|(U_1(z+c) - U_2(z+c))(U_1(z+c) + U_2(z+c))''\| \\ &\le 4\|U_1 - U_2\|_2\|U_1 + U_2\|_2 \end{aligned} $$

for $c = -p, 0, p$, and

$$ \begin{aligned} \|(U_1^2(z) - U_2^2(z))'\| &\le \|(U_1(z) - U_2(z))'(U_1(z) + U_2(z))\| \\ &\quad + \|(U_1(z) - U_2(z))(U_1(z) + U_2(z))'\| \le 2\|U_1 - U_2\|_2\|U_1 + U_2\|_2 \end{aligned} $$

since $a^2 - b^2 = (a-b)(a+b)$ for any $a, b \in \mathbb{C}$. Now, (4) follows easily. $\square$

Next, if $U \in Z$ with $U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k_1 z}$ then

$$ \mathcal{K}U(z) = \sum_{k \in \mathbb{Z}} (1 - \omega^2 k^2 + 2\lambda\omega^2 k^2 \cos kp + i\gamma\omega k) c_k e^{k_1 z} \quad (6) $$

and so $\mathcal{K} \in L(Z, X)$ with
$$ \|\mathcal{K}\|_{L(Z,X)} \leq 1 + \omega^2(1 + 2\lambda) + \gamma\omega. $$

If

$$ \Theta := \inf_{k \in \mathbb{Z} \setminus \{0\}} \left\{ 1, \sqrt{\left(\frac{1}{k^2} + \omega^2(2\lambda \cos kp - 1)\right)^2 + \frac{\gamma^2 \omega^2}{k^2}} \right\} > 0 \quad (7) $$

is a constant depending on $\gamma$, $\lambda$, $\omega$ and $p$, then we also have $\mathcal{K}^{-1} \in L(X, Z) \subset L(X, Y) \subset L(X)$. So $\mathcal{K}^{-1}: X \to Z$ is continuous such that

$$ \|\mathcal{K}^{-1}\|_{L(X,Z)} \leq \frac{1}{\Theta}. \quad (8) $$

Now we can prove the following existence result on (2) when all parameters except $h$ are fixed.

**Theorem 1.** Assume (7) along with

$$ 0 < \|h\| < \frac{\Theta^2}{8\omega(2\omega + 4\lambda\omega + \gamma)}. \quad (9) $$

Then equation (2) has a unique solution $U(h) \in B(\rho_h)$ in a closed ball

$$ B(\rho_h) := \{U \in Z \mid \|U\|_2 \leq \rho_h\}, $$

where

$$ \rho_h = \frac{\Theta - \sqrt{\Theta^2 - 8\omega(2\omega + 4\lambda\omega + \gamma)\|h\|}}{4\omega(2\omega + 4\lambda\omega + \gamma)}. \quad (10) $$

Moreover, $U(h)$ can be approximated by an iteration process. Finally, it holds

$$ \|U(h_1) - U(h_2)\|_2 \leq \frac{\|h_1 - h_2\|}{\sqrt{\Theta^2 - 8\omega(2\omega + 4\lambda\omega + \gamma)\max\{\|h_1\|, \|h_2\|\}}} \quad (11) $$

for any $h_1, h_2 \in X$ satisfying (9).

**Proof.** We rewrite (2) as a parameterized fixed point problem

$$ U = \mathcal{R}(U, h) := \mathcal{K}^{-1}\mathcal{F}(U, h) $$

in $Z$. We already know that $\mathcal{R}: Z \times X \to Z$ is continuous and by (3), (8) such that

$$ \|\mathcal{R}(U, h)\|_2 \leq \frac{1}{\Theta} (2\omega(2\omega + 4\lambda\omega + \gamma)\|U\|_2^2 + \|h\|). $$

Next, if there is $\rho_h > 0$ such that
$$A(\rho_h) := \Theta\rho_h - 2\omega(2\omega + 4\lambda\omega + \gamma)\rho_h^2 = \|h\|, \quad (12)$$

then $\mathcal{R}(\cdot, h)$ maps $B(\rho_h)$ into itself. So it remains to study (12). In order to find $H > 0$ – the largest right-hand side of (12) when this equation has a solution $\rho_H > 0$, we need to solve $A(r) = H$ together with $DA(r) = 0$. This implies

$$\rho_H = \frac{\Theta}{4\omega(2\omega + 4\lambda\omega + \gamma)}$$

and (9). So assuming (7), (9), we know that (12) has a positive solution $\rho_h < \rho_H$. We take the smallest one which clearly has the form (10). So $\mathcal{R}(\cdot, h)$ maps $B(\rho_h)$ into itself and, moreover, by (4), (8)

$$\|\mathcal{R}(U_1, h) - \mathcal{R}(U_2, h)\|_2 \le \frac{4\omega(2\omega + 4\lambda\omega + \gamma)\|U_1 - U_2\|_2\rho_h}{\Theta}$$

for any $U_1, U_2 \in B(\rho_h)$. Hence, $\mathcal{R}(\cdot, h)$ is a contraction on $B(\rho_h)$ with a contraction constant

$$\frac{4\omega(2\omega + 4\lambda\omega + \gamma)\rho_h}{\Theta} < \frac{4\omega(2\omega + 4\lambda\omega + \gamma)\rho_H}{\Theta} = 1.$$

The proof of the existence and uniqueness is finished by the Banach fixed point theorem [2]. Next, let $h_1, h_2 \in X$ satisfy (9). Then $U(h_i) \in B(\rho_{h_i}) \subset B(\rho_{\tilde{H}})$ for $i = 1, 2$ and $\tilde{H} := \max\{\|h_1\|, \|h_2\|\}$. Note $\tilde{H}$ satisfies (9). By (4), (5) and (8), we derive

$$\begin{aligned} \|U(h_1) - U(h_2)\|_2 &= \|\mathcal{R}(U(h_1), h_1) - \mathcal{R}(U(h_2), h_2)\|_2 \\ &\le \|\mathcal{R}(U(h_1), h_1) - \mathcal{R}(U(h_2), h_1)\|_2 + \|\mathcal{R}(U(h_2), h_1) - \mathcal{R}(U(h_2), h_2)\|_2 \\ &\le \frac{4\omega(2\omega + 4\lambda\omega + \gamma)\|U(h_1) - U(h_2)\|_2\rho_{\tilde{H}} + \|h_1 - h_2\|}{\Theta} \end{aligned}$$

which implies (11). $\square$

**Remark 1.** In the following special cases, (7) holds and we can replace $\Theta$ with the corresponding $\Theta_i$ in the above considerations:

1. If $p \in 2\pi\mathbb{Z}$, $\lambda > \frac{1}{2}$, then

$$\Theta \ge \Theta_1 := \min\{1, \omega^2(2\lambda - 1)\} > 0.$$

2. Let $\gamma = 0$. If $p \in \pi\mathbb{Q}$, then we can write $p = \frac{p_1}{p_2}\pi$ for some $p_1 \in \mathbb{Z}$, $p_2 \in \mathbb{N}$, where $p_1, p_2$ are relatively prime (their only common divisor is 1). Moreover,

$$M_1 := \{\cos kp \mid k \in \mathbb{Z} \setminus \{0\}\} \subset \left\{\cos \frac{k\pi}{p_2} \mid k \in \{0, 1, \dots, 2p_2 - 1\}\right\} =: M_2. \quad (13)$$

Indeed, for each $k \in \mathbb{Z}$ there exist $i, j \in \mathbb{Z}$ such that $0 \le i \le 2p_2 - 1$ and $kp_1 = 2jp_2 + i$. Hence, for $\cos kp \in M_1$ we have
$$ \cos kp = \cos \frac{kp_1\pi}{p_2} = \cos \frac{i\pi}{p_2} \in M_2. $$

To find a better relationship between $M_1$ and $M_2$, we need to solve $\cos \frac{\tilde{k}_1 p_1}{p_2} \pi = \cos \frac{k_2 \pi}{p_2}$ for $\tilde{k}_1 \in \mathbb{Z} \setminus \{0\}$ and $k_2 \in \{0, 1, \dots, 2p_2 - 1\}$. This is equivalent to $\frac{k_1 p_1}{p_2} \pi = \frac{k_2 \pi}{p_2} + 2\pi l$ for $l \in \mathbb{Z}$ and $k_1 = \pm \tilde{k}_1$, i.e., $k_1 p_1 = k_2 + 2lp_2$. If $p_1$ is odd, then $p_1$ and $2p_2$ are relatively prime. Thus, we know [3] that there exist $\hat{k}, \hat{l} \in \mathbb{Z}$ such that $\hat{k}p_1 = 1 + 2\hat{l}p_2$. Consequently, $\hat{k}k_2p_1 = k_2 + 2\hat{l}k_2p_2$ and $M_1 = M_2$. If $p_1$ is even, then $p_1 = 2\bar{p}_1$ and $k_1 2\bar{p}_1 = k_2 + 2lp_2$, so $k_2 = 2\bar{k}_2$ and we have $k_1 \bar{p}_1 = \bar{k}_2 + lp_2$. Clearly $\bar{p}_1$ and $p_2$ are relatively prime. Thus, there exist $\bar{k}, \bar{l} \in \mathbb{Z}$ such that $\bar{k}\bar{p}_1 = 1 + \bar{l}p_2$. Then $\bar{k}\bar{k}_2\bar{p}_1 = \bar{k}_2 + \bar{l}\bar{k}_2p_2$. Hence $M_1 = \{\cos \frac{2k\pi}{p_2} | k \in \{0, 1, \dots, p_2 - 1\} \} \subsetneq M_2$.

$$ \omega^2 \notin \left\{ \frac{1}{k^2(1 - 2\lambda \cos kp)} \right\}_{k \in \mathbb{Z} \setminus \{0\}} \quad (14) $$

and

$$ 0 \notin 2\lambda M_2 - 1 \quad (15) $$

for $M_2$ defined in (13), we get the existence of $\delta > 0$ such that

$$ \sqrt{\left(\frac{1}{k^2} + \omega^2(2\lambda \cos kp - 1)\right)^2 + \frac{\gamma^2 \omega^2}{k^2}} \ge \omega^2 \left|\frac{1}{\omega^2 k^2} + 2\lambda \cos kp - 1\right| \ge \delta > 0 $$

for each $k \in \mathbb{Z} \setminus \{0\}$. Thus, if $p \in \pi \mathbb{Q}$ and (14), (15) are valid, then

$$ \Theta \ge \min \left\{ 1, \inf_{k \in \mathbb{Z} \setminus \{0\}} \omega^2 \left| \frac{1}{\omega^2 k^2} + 2\lambda \cos kp - 1 \right| \right\} =: \Theta_2 > 0. $$

Note that $\frac{1}{k^2(1-2\lambda \cos kp)} \to 0$ if $|k| \to \infty$. So if $p \in \pi \mathbb{Q}$ and (15) holds, there is at most a finite number of resonant modes $k_0 \in \mathbb{Z} \setminus \{0\}$ determined by equation

$$ \omega^2 = \frac{1}{k_0^2(1 - 2\lambda \cos k_0 p)}. \quad (16) $$

More precisely, for given $\lambda, \omega$ there is not more than $2(p_2 + 1)$ resonant modes $k_0$. Indeed, let $\{k_{0j}\}_{j=1}^J$ be an increasing sequence of positive resonant modes corresponding to the fixed $\lambda, \omega$. Then $\{1 - 2\lambda \cos k_{0j} p\}_{j=1}^J$ has to be decreasing, since

$$ k_{0j}^2 (1 - 2\lambda \cos k_{0j} p) = \frac{1}{\omega^2} = \text{const.} \quad \forall j = 1, \dots, J, $$

i.e., $\{\cos k_{0j} p\}_{j=1}^J$ is increasing. From properties of $\cos x$ and (13), it follows that the longest sequence of the form $\cos k_{0j} p$ has the values $\{\cos \frac{k\pi}{p_2} | k = p_2, \dots, 2p_2\}$ which has $p_2 + 1$ elements. Finally, to each $k_{0j}$ corresponds resonant mode $-k_{0j}$.
3. If $\gamma > 0$, condition (7) is satisfied even without the non-resonance condition (14), since in each resonant mode $k_0$ we have

$$ \sqrt{\left(\frac{1}{k_0^2} + \omega^2(2\lambda \cos k_0 p - 1)\right)^2 + \frac{\gamma^2 \omega^2}{k_0^2}} = \frac{\gamma \omega}{|k_0|} > 0 $$

and there is only the finite number of resonant modes. For each non-resonant mode $k$ it holds

$$ \sqrt{\left(\frac{1}{k^2} + \omega^2(2\lambda \cos kp - 1)\right)^2 + \frac{\gamma^2 \omega^2}{k^2}} \geq \left|\frac{1}{k^2} + \omega^2(2\lambda \cos kp - 1)\right| > 0. $$

Summarizing, if $\gamma > 0$, $p \in \pi\mathbb{Q}$ and (15) holds, condition (7) is fulfilled and

$$ \Theta \geq \min \left\{ 1, \min_{k \in M} \frac{\gamma \omega}{|k|}, \inf_{k \in \mathbb{Z}\setminus\{0\} \cup M} \omega^2 \left| \frac{1}{\omega^2 k^2} + 2\lambda \cos kp - 1 \right| \right\} =: \Theta_3 > 0 $$

where $M$ is the set of resonant modes $k_0$.

### 3. Simple resonances

In this section, we investigate the bifurcation of a small solution of (2) under the assumption of a simple resonance described in Remark 1(2) So, we suppose

(R) $p \in \pi\mathbb{Q}$, condition (15) is valid, for some $k_0 \in \mathbb{Z}\setminus\{0\}$ equation (16) holds and has the only solutions $\pm k_0$,

and we consider the equation

$$ A(U, \varepsilon) := KU - \mathcal{F}(U) - \varepsilon\mathcal{G}(U) = 0 \quad (17) $$

with

$$ \begin{align} KU &:= \omega^2(U''(z) - \lambda U''(z-p) - \lambda U''(z+p)) + U(z) \nonumber \\ \mathcal{F}(U) &:= -\omega^2((U^2(z))'' - \lambda(U^2(z-p))'' - \lambda(U^2(z+p))'') \nonumber \\ \mathcal{G}(U)(z) &:= -\gamma\omega(U'(z) + (U^2(z))') + h(z) \tag{18} \end{align} $$

for $U$ and $\varepsilon$ small.

We shall solve equation (17) via Lyapunov-Schmidt reduction method [4]. Clearly, $A(0, 0) = 0$ and $D_U A(0, 0) = K$. For simplicity we denote $\kappa := |k_0|$ and $R_1^\kappa := \mathcal{N}K, R_2^\kappa := \mathcal{R}K$ the null space and the range of the operator $K \in L(Z, X)$, respectively. Equation (6) with $\gamma = 0$ together with (R) yields that if $U \in Z$ with $U(z) = \sum_{k \in \mathbb{Z}} c_k e^{k_1 z}$, then

$$ KU(z) = \sum_{k \in \mathbb{Z} \setminus \{\pm k_0\}} (1 - \omega^2 k^2 + 2\lambda\omega^2 k^2 \cos kp) c_k e^{k_1 z}. $$
Consequently,

$$
R_1^\kappa = \{c e^{k_1 z} + \bar{c} e^{-k_1 z} \mid c \in \mathbb{C}\}, \quad R_2^\kappa = \left\{ U \in X \mid U(z) = \sum_{k \in \mathbb{Z} \setminus \{\pm k_0\}} c_k e^{k_1 z} \right\}.
$$

Let $\mathcal{Q} : X \to \mathbb{R}_2^\kappa$ be the projection onto $\mathbb{R}_2^\kappa \subset X$ given by

$$
\mathcal{Q}U(z) = \sum_{k \in \mathbb{Z} \setminus \{\pm k_0\}} c_k e^{k_1 z}.
$$

Now, we take the decomposition $U = U_1 + U_2$, $U_1 \in \mathbb{R}_1^\kappa$, $U_2 \in \mathbb{R}_2^\kappa \cap Z$ for $U \in Z$ and decou-
ple (17) to equations

$$
QA(U_1 + U_2, \varepsilon) = 0, \tag{19}
$$

$$
(\mathbb{I} - \mathcal{Q}) A(U_1 + U_2, \varepsilon) = 0. \tag{20}
$$

Applying the implicit function theorem to the first equation implies the existence of neighbour-
hoods $V_1 \times V_0 \subset \mathbb{R}_1^\kappa \times \mathbb{R}$ of the point $(0,0)$ in $\mathbb{R}_1^\kappa \times \mathbb{R}$, $V_2 \subset \mathbb{R}_2^\kappa \cap Z$ of $0$ in $\mathbb{R}_2^\kappa \cap Z$ and a
unique $C^\infty$-function $U_2: V_1 \times V_0 \to V_2$ such that $QA(U_1 + \tilde{U}_2, \varepsilon) = 0$ for $(U_1, \varepsilon) \in V_1 \times V_0$
and $\tilde{U}_2 \in V_2$ if and only if $\tilde{U}_2 = U_2(U_1, \varepsilon)$. Moreover, $U_2(0, 0) = 0$.

Thus, equation (20) has the form

$$
(\mathbb{I} - \mathcal{Q}) A(U_1 + U_2(U_1, \varepsilon), \varepsilon) = 0 \quad (21)
$$

for $U_1 \in V_1$ and $\varepsilon \in V_0$. Differentiating equation (19) we obtain

$$
D_{U_1} U_2(0, 0) = 0, \quad D_{\varepsilon} U_2(0, 0) = -\mathcal{K}^{-1} Q D_{\varepsilon} A(0, 0) = \mathcal{K}^{-1} Q \mathcal{G}(0). \quad (22)
$$

Hence

$$
U_2(U_1, \varepsilon) = \mathcal{K}^{-1} \mathcal{Q} \mathcal{G}(0) \varepsilon + \mathcal{O}(\|(U_1, \varepsilon)\|^2). \quad (23)
$$

Expanding $\mathcal{F}$ and $\mathcal{G}$ of (18) into Taylor series using the last expansion of $U_2(U_1, \varepsilon)$ yields

$$
\begin{align*}
\mathcal{F}(U_1 + U_2(U_1, \varepsilon)) &= \frac{1}{2} \mathrm{D}^2 \mathcal{F}(0) (U_1 + \mathrm{D}_{\varepsilon} U_2(0,0)\varepsilon)^2 + \mathcal{O}(\|(U_1, \varepsilon)\|^3), \\
\mathcal{G}(U_1 + U_2(U_1, \varepsilon)) &= \mathcal{G}(0) + \mathrm{D}\mathcal{G}(0) (U_1 + \mathrm{D}_{\varepsilon} U_2(0,0)\varepsilon) + \mathcal{O}(\|(U_1, \varepsilon)\|^2).
\end{align*}
$$

Note that $D^2\mathcal{F}(0)V^2 = 2\mathcal{F}(V)$ and $\varepsilon O(\|(U_1, \varepsilon)\|^2) = O(\|(U_1, \varepsilon)\|^3)$. Thus equation (21) is equivalent to

$$
(\mathbb{I} - \mathcal{Q})(\mathcal{F}(U_1 + \mathcal{D}_{\varepsilon} U_2(0,0)\varepsilon) + \varepsilon\mathcal{G}(0) + \varepsilon\mathcal{D}\mathcal{G}(0)(U_1 + \mathcal{D}_{\varepsilon}U_2(0,0)\varepsilon) + \mathcal{O}(\|(U_1, \varepsilon)\|^3)) = 0
$$

or, using $D_\epsilon U_2(0, 0)$ and $\mathcal{F}$ of (18),
$$ (\mathbb{I} - \mathcal{Q})(\omega^2(((U_1(z) + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z)\varepsilon)^2)'' - \lambda((U_1(z-p) + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z-p)\varepsilon)^2)'' \\ - \lambda((U_1(z+p) + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z+p)\varepsilon)^2)'' - \varepsilon\mathcal{Q}(0)(z) \\ - \varepsilon\mathcal{D}\mathcal{G}(0)(z)(U_1 + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)\varepsilon) + O(\|(U_1, \varepsilon)\|^3)) = 0. \quad (24) $$

Note that

$$ ((U_1(z) + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z)\varepsilon)^2)'' = 2(U_1'(z) + (\mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z))'\varepsilon)^2 \\ + 2(U_1(z) + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z)\varepsilon)(U_1''(z) + (\mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z))''\varepsilon) $$

for any $z \in \mathbb{R}$. Let $U_1(z) = c e^{k_0 z} + \tilde{c} e^{-k_0 z} \in R_1^\kappa$ and $h(z) = \sum_{k \in \mathbb{Z}} h_k e^{k_1 z} \in X$. Then

$$ U_1'(z) = \tilde{c} e^{k_0 z} + \tilde{c} e^{-k_0 z} \in R_1^\kappa, \quad \mathcal{K}^{-1} \mathcal{Q} \mathcal{G}(0)(z) = \sum_{k \in \mathbb{Z} \setminus \{\pm k_0\}} d_k e^{k_1 z} \in R_2^\kappa, \\ (\mathcal{K}^{-1} \mathcal{Q} \mathcal{G}(0)(z))' = \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0\}} \tilde{d}_k e^{k_1 z} \in R_2^\kappa \quad (25) $$

where

$$ \tilde{c} = \iota k_0 c, \quad d_k = \frac{h_k}{1 - \omega^2 k^2 + 2\lambda\omega^2 k^2 \cos k p}, \quad \tilde{d}_k = \iota k d_k. $$

Clearly, $d_{-k} = \overline{d_k}$, $\tilde{d}_{-k} = \overline{\tilde{d}_k}$ for each $k \in \mathbb{Z} \setminus \{0, \pm k_0\}$. It can be easily shown that

$$ (\mathbb{I} - \mathcal{Q})(U_1'(z))^2 = 0 = (\mathbb{I} - \mathcal{Q})(U_1(z)U_1''(z)), \\ (\mathbb{I} - \mathcal{Q})(U_1'(z))(\mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z))' = \tilde{c}\tilde{d}_{2k_0} e^{-k_0 z} + \tilde{c}\tilde{d}_{2k_0} e^{k_0 z}. \quad (26) $$

Next, setting $\tilde{d}_{\pm k_0} = \tilde{d}_0 = 0$,

$$ \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0\}} \tilde{d}_k e^{k_1 z} \quad \sum_{j \in \mathbb{Z} \setminus \{0, \pm k_0\}} \tilde{d}_j e^{j_1 z} = \sum_{l \in \mathbb{Z}} \sum_{k+j=l} \tilde{d}_k \tilde{d}_j e^{l_1 z}. $$

Hence the coefficient by $e^{k_0 z}$ in $((\mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z))')^2$ is

$$ b := \sum_{k+j=k_0} \tilde{d}_k \tilde{d}_j = \sum_{k \in \mathbb{Z}} \tilde{d}_k \tilde{d}_{k_0-k} = \sum_{k \in \mathbb{Z}} \overline{\tilde{d}_k \tilde{d}_{k-k_0}} = \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \overline{\tilde{d}_k \tilde{d}_{k-k_0}}. $$

Therefore

$$ (\mathbb{I} - \mathcal{Q})(((\mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z))')^2) = b e^{k_0 z} + \bar{b} e^{-k_0 z}. \quad (27) $$

Analogically, denoting

$$ \tilde{\overline{d}}_k := -k^2 d_k, \quad k \in \mathbb{Z} \setminus \{0, \pm k_0\}, \quad \tilde{c} := -k_0^2 c, \quad \tilde{b} := \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} d_k \overline{\tilde{d}_{k-k_0}} $$
we have

$$U_1''(z) = \tilde{c} e^{k_0 t z} + \bar{\tilde{c}} e^{-k_0 t z} \in R_1^\kappa, \quad (\mathcal{K}^{-1} Q G(0)(z))'' = \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0\}} \tilde{d}_k e^{k_0 t z} \in R_2^\kappa,$$

and the next identities can be proved

$$\begin{align}
(\mathbb{I} - Q)(U_1(z)(\mathcal{K}^{-1} Q G(0)(z)))'' &= \bar{c} \tilde{d}_{2k_0} e^{k_0 t z} + \overline{c \tilde{d}_{2k_0}} e^{-k_0 t z}, \nonumber \\
(\mathbb{I} - Q)(U_1''(z) \mathcal{K}^{-1} Q G(0)(z)) &= (\bar{\tilde{c}} \tilde{d}_{2k_0} + \tilde{c} \tilde{d}_0) e^{k_0 t z} + (\overline{\tilde{c} \tilde{d}_{2k_0}} + \overline{\tilde{c} \tilde{d}_0}) e^{-k_0 t z}, \nonumber \\
(\mathbb{I} - Q)(\mathcal{K}^{-1} Q G(0)(z)(\mathcal{K}^{-1} Q G(0)(z)))'' &= \tilde{b} e^{k_0 t z} + \overline{\tilde{b}} e^{-k_0 t z}. \tag{28}
\end{align}$$

Moreover, by definition of $\mathcal{G}$ in (18),

$$\begin{align}
(\mathbb{I} - Q)\mathcal{G}(0)(z) &= h_{k_0} e^{k_0 t z} + \overline{h_{k_0}} e^{-k_0 t z}, \nonumber \\
(\mathbb{I} - Q)\mathcal{D}\mathcal{G}(0)(U_1 + \mathcal{K}^{-1}\mathcal{Q}\mathcal{G}(0)(z)\varepsilon) &= -\gamma\omega(\tilde{c}e^{k_0 t z} + \overline{\tilde{c}}e^{-k_0 t z}). \tag{29}
\end{align}$$

Now, we put (26), (27), (28), (29) in (24) and collect coefficients by $e^{k_0 t z}$ and $e^{-k_0 t z}$ to obtain

$$H e^{k_0 t z} + \overline{H} e^{-k_0 t z} + O(\|(U_1, \varepsilon)\|^3) = 0 \quad (30)$$

where

$$H = 2\omega^2 \left( 2\bar{c}\tilde{d}_{2k_0} + \bar{c}\tilde{d}_{2k_0} + \tilde{c}\tilde{d}_{2k_0} + \tilde{c}d_0 \right) (1 - \lambda e^{-k_0 t p} - \lambda e^{k_0 t p})\varepsilon - (h_{k_0} - \gamma\omega\tilde{c})\varepsilon + 2\omega^2(b + \tilde{b})(1 - \lambda e^{-k_0 t p} - \lambda e^{k_0 t p})\varepsilon^2.$$

Equation (30) is equivalent to

$$H + O(\|(c, \bar{c}, \varepsilon)\|^3) = 0. \quad (31)$$

Using resonance condition (16), we get

$$(1 - \lambda e^{-k_0 t p} - \lambda e^{k_0 t p}) = 1 - 2\lambda \cos k_0 p = \frac{1}{\omega^2 k_0^2}.$$

Dropping tildes yields

$$2\bar{c}\tilde{d}_{2k_0} + \bar{c}\tilde{d}_{2k_0} + \tilde{c}\tilde{d}_{2k_0} + \tilde{c}d_0 = -k_0^2\bar{c}\tilde{d}_{2k_0} - k_0^2cd_0, \\ b + \tilde{b} = \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} k_0(k - k_0)d_k\overline{d_{k-k_0}}.$$

Thus we have to solve
$$ -(2(\bar{c}d_{2k_0} + cd_0) + h_{k_0} - \gamma\omega\bar{c})\varepsilon + \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} k_0(k-k_0)d_k \overline{d_{k-k_0}} \\ + O(\|(c, \bar{c}, \varepsilon)\|^3) = 0. $$

Since $c, h_k \in \mathbb{C}$, we have

$$ c = x + \iota y, \quad h_k = \mu_k + \iota v_k, \quad d_k = \frac{\mu_k + \iota v_k}{\delta_k} $$

where $x, y, \mu_k, v_k \in \mathbb{R}$ and $\delta_k = 1 - \omega^2 k^2 + 2\lambda\omega^2 k^2 \cos kx \in \mathbb{R}$. Obviously, $\delta_0 = 1, h_0 = d_0 = \mu_0 \in \mathbb{R}$ and $\mu_k = \mu_{-k}, v_k = -v_{-k} \forall k \in \mathbb{N}$. So, we get

$$ -\left( \frac{2(x - \iota y)(\mu_{2k_0} + \iota v_{2k_0})}{\delta_{2k_0}} + 2(x + \iota y)\mu_0 + \mu_{k_0} + \iota v_{k_0} - \gamma\omega\iota k_0(x + \iota y) \right)\varepsilon \\ + \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_k + \iota v_k)(\mu_{k-k_0} - \iota v_{k-k_0})}{\delta_k \delta_{k-k_0}} + O(\|(x, y, \varepsilon)\|^3) = 0 $$

or, separating the real and imaginary parts,

$$ -\left( \frac{2(x \mu_{2k_0} + y \nu_{2k_0})}{\delta_{2k_0}} + 2x\mu_0 + \mu_{k_0} + \gamma\omega k_0 y \right)\varepsilon \\ + \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_k \mu_{k-k_0} + \nu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + O(\|(x, y, \varepsilon)\|^3) = 0, \\ -\left( \frac{2(x \nu_{2k_0} - y \mu_{2k_0})}{\delta_{2k_0}} + 2y\mu_0 + \nu_{k_0} - \gamma\omega k_0 x \right)\varepsilon \\ + \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_{k-k_0} \nu_k - \mu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + O(\|(x, y, \varepsilon)\|^3) = 0. $$

Multiplying both equations by $(-1)$ and collecting terms by $x$ and $y$ we obtain

$$ \begin{aligned} & 2\left(\frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0\right)\varepsilon x + \left(\frac{2\nu_{2k_0}}{\delta_{2k_0}} + \gamma\omega k_0\right)\varepsilon y + \mu_{k_0}\varepsilon \\ & -\frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_k \mu_{k-k_0} + \nu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + O(\|(x, y, \varepsilon)\|^3) = 0, \\ & \left(\frac{2\nu_{2k_0}}{\delta_{2k_0}} - \gamma\omega k_0\right)\varepsilon x + 2\left(-\frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0\right)\varepsilon y + \nu_{k_0}\varepsilon \\ & -\frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_{k-k_0} \nu_k - \mu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + O(\|(x, y, \varepsilon)\|^3) = 0. \end{aligned} \quad (32) $$

Now we intend to derive the cubic $O(\|(U_1, 0)\|^3)$ of (30) as follows. We denote by $B := -\frac{1}{2}D^2\mathcal{F}(0)$ and $u_2(U_1) := U_2(U_1, 0)$. Then (19) and (20) are equivalent to
$$Q(Ku_2(U_1) + B(U_1 + u_2(U_1))^2) = 0, \quad (33)$$

$$(\mathbb{I} - Q)B(U_1 + u_2(U_1))^2 = 0. \quad (34)$$

First note that $u_2(U_1) = O(\|U_1\|^2)$ by (23), $(\mathbb{I} - Q)B U_1^2 = 0$ and then (34) has the form

$$\begin{aligned} 0 = (\mathbb{I} - Q)B(U_1 + u_2(U_1))^2 &= (\mathbb{I} - Q)B U_1^2 + 2(I - Q)B(U_1, u_2(U_1)) + (I - Q)B(u_2(U_1))^2 \\ &= (I - Q)B(U_1, D_{U_1}^2 u_2(0)U_1^2) + O(\|U_1\|^4), \end{aligned}$$

hence $O(\|(U_1, 0)\|^3) = (I - Q)B(U_1, D_{U_1}^2 u_2(0)U_1^2)$. To derive $D_{U_1}^2 u_2(0)$, we twice differentiate (33) to obtain

$$\mathcal{K}D_{U_1}u_2(U_1)V_1 + 2QB(U_1 + u_2(U_1), V_1 + D_{U_1}u_2(U_1)V_1) = 0$$

and then

$$\mathcal{K}D_{U_1}^2 u_2(0)(V_1, V_2) + 2QB(V_1, V_2) = 0,$$

hence

$$D_{U_1}^2 u_2(0) U_1^2 = -2 \mathcal{K}^{-1} Q B U_1^2 = 2 \mathcal{K}^{-1} Q \mathcal{F}(U_1),$$

consequently, we arrive at

$$O(\|(U_1, 0)\|^3) = (I - Q)B(U_1, D_{U_1}^2 u_2(0)U_1^2) = -(I - Q)D^2\mathcal{F}(0)(U_1, \mathcal{K}^{-1}Q\mathcal{F}(U_1)). \quad (35)$$

Let $U_1(z) = c e^{k_0 tz} + \bar{c} e^{-k_0 tz} \in \mathbb{R}_1^\kappa$. Then

$$\begin{aligned} \mathcal{K}^{-1} Q \mathcal{F}(U_1) &= \mathcal{K}^{-1} Q \left( 4\omega^2 k_0^2 (1 - 2\lambda \cos 2k_0 p) \left( c^2 e^{2k_0 tz} + \bar{c}^2 e^{-2k_0 tz} \right) \right) \\ &= \frac{1 - \delta_{2k_0}}{\delta_{2k_0}} c^2 e^{2k_0 tz} + \frac{1 - \delta_{2k_0}}{\delta_{2k_0}} \bar{c}^2 e^{-2k_0 tz}. \end{aligned}$$

Next we derive

$$\begin{aligned} D^2\mathcal{F}(0)(V, W) &= -2\omega^2((V(z)W(z))'' - \lambda(V(z-p)W(z-p))'' - \lambda(V(z+p)W(z+p))'') \\ &= -2\omega^2(V(z)W''(z) + V''(z)W(z) + 2V'(z)W'(z)) \\ &\quad - \lambda(V(z-p)W''(z-p) + V''(z-p)W(z-p) + 2V'(z)W'(z)) \\ &\quad - \lambda(V(z+p)W''(z+p) + V''(z+p)W(z+p) + 2V'(z)W'(z)). \end{aligned}$$

Denoting $\mu = \frac{1 - \delta_{2k_0}}{\delta_{2k_0}}$ we get
$$
\begin{align*}
U_1(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)]'' &= -4k_0^2 U_1(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)] \\
&= -4k_0^2 (\mu c^3 e^{3k_0 tz} + \mu \bar{c}^3 e^{-3k_0 tz} + \mu c^2 \bar{c} e^{k_0 tz} + \mu c \bar{c}^2 e^{-k_0 tz}),
\end{align*}
$$

$$
\begin{align*}
U_1''(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)] &= -k_0^2 U_1(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)] \\
&= -k_0^2 (\mu c^3 e^{3k_0 tz} + \mu \bar{c}^3 e^{-3k_0 tz} + \mu c^2 \bar{c} e^{k_0 tz} + \mu c \bar{c}^2 e^{-k_0 tz}),
\end{align*}
$$

$$
\begin{align*}
U_1'(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)]' &= (k_0 t c e^{k_0 t z} - k_0 t \bar{c} e^{-k_0 t z}) (2k_0 t \mu c^2 e^{2k_0 t z} - 2k_0 t \mu \bar{c}^2 e^{-2k_0 t z}) \\
&= -2k_0^2 \mu c^3 e^{3k_0 t z} - 2k_0^2 \mu \bar{c}^3 e^{-3k_0 t z} + 2k_0^2 \mu c^2 \bar{c} e^{k_0 t z} + 2k_0^2 \mu c \bar{c}^2 e^{-k_0 t z}.
\end{align*}
$$

So

$$
(I - \mathcal{Q}) (U_1(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)]'' + 2U_1'(z)[\mathcal{K}^{-1} \mathcal{Q} \mathcal{F}(U_1)(z)]') = 0.
$$

Hence

$$
\begin{align*}
-(I - \mathcal{Q})\mathrm{D}^2\mathcal{F}(0)(U_1, \mathcal{K}^{-1}\mathcal{Q}\mathcal{F}(U_1)) &= 2\omega^2(\mu c^2 \bar{c} e^{k_0 t z} (-k_0^2) + \mu c \bar{c}^2 e^{-k_0 t z} (-k_0^2) \\
&\quad - \lambda \mu c^2 \bar{c} e^{k_0 t (z-p)} (-k_0^2) - \lambda \mu c \bar{c}^2 e^{-k_0 t (z-p)} (-k_0^2) - \lambda \mu c^2 \bar{c} e^{k_0 t (z+p)} (-k_0^2) \\
&\quad - \lambda \mu c \bar{c}^2 e^{-k_0 t (z+p)} (-k_0^2)) \\
&= -2k_0^2 \omega^2 \mu |c|^2 (U_1(z) - \lambda U_1(z-p) - \lambda U_1(z+p)) \\
&= -2k_0^2 \omega^2 \mu |c|^2 (1 - \lambda e^{k_0 t p} - \lambda e^{-k_0 t p}) U_1(z) \\
&= -2\mu |c|^2 U_1(z) = -\frac{2(1 - \delta_{2k_0})}{\delta_{2k_0}} |c|^2 (c e^{k_0 t z} + \bar{c} e^{-k_0 t z}).
\end{align*}
$$

Thus, by (35) and the last identity, in (31) we have

$$
O(\|(c, \bar{c}, 0)\|^3) = - \frac{2c|c|^2(1 - \delta_{2k_0})}{\delta_{2k_0}}.
$$

Summarizing, the bifurcation equation (32) has a form

$$
\begin{align*}
& 2 \left( \frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) \varepsilon x + \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} + \gamma \omega k_0 \right) \varepsilon y + \mu_{k_0} \varepsilon + \frac{2(1 - \delta_{2k_0})x(x^2 + y^2)}{\delta_{2k_0}} \\
& - \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_k \mu_{k-k_0} + \nu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + \varepsilon O(\|(x, y, \varepsilon)\|^2) = 0, \\
& \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} - \gamma \omega k_0 \right) \varepsilon x + 2 \left( -\frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) \varepsilon y + \nu_{k_0} \varepsilon + \frac{2(1 - \delta_{2k_0})y(x^2 + y^2)}{\delta_{2k_0}} \\
& - \frac{2\varepsilon^2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_{k-k_0} \nu_k - \mu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}} + \varepsilon O(\|(x, y, \varepsilon)\|^2) = 0. \tag{36}
\end{align*}
$$
Now, we scale (36)

$$x \leftrightarrow \varepsilon x, \quad y \leftrightarrow \varepsilon y, \quad \varepsilon \leftrightarrow \varepsilon^3, \quad \gamma \leftrightarrow \gamma \quad (37)$$

and divide both equations by $\varepsilon^3$. This leads to

$$\begin{align}
\mu_{k_0} + \frac{2(1 - \delta_{2k_0})x(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) &= 0 \nonumber \\
\nu_{k_0} + \frac{2(1 - \delta_{2k_0})y(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) &= 0 \tag{38}
\end{align}$$

and the following result holds.

**Theorem 2.** Assume (R) and $(\mu_{k_0}, \nu_{k_0}) \neq (0, 0)$. Then equation (17) has a solution $U \in Z$ close to 0 for any $\varepsilon \neq 0$ sufficiently small.

**Proof.** Let us denote

$$\tilde{H}_1(x, y, \varepsilon) := \left( \mu_{k_0} + \frac{2(1 - \delta_{2k_0})x(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon), \nu_{k_0} + \frac{2(1 - \delta_{2k_0})y(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) \right).$$

Implicit function theorem applied to equation (38), i.e. $\tilde{H}_1(x, y, \varepsilon) = 0$ gives the existence of neighbourhoods $W_0 \subset V_0 \subset \mathbb{R}$ of the point 0 in $\mathbb{R}$, $W \subset \mathbb{R}^2$ of $-\sqrt[3]{\frac{\delta_{2k_0}}{2(1-\delta_{2k_0})(\mu_{k_0}^2+\nu_{k_0}^2)}}(\mu_{k_0}, \nu_{k_0})$ in $\mathbb{R}^2$ and a unique pair of continuous functions $x(\varepsilon)$, $y(\varepsilon)$ defined in $W_0$ such that $\tilde{H}_1(\tilde{x}, \tilde{y}, \varepsilon) = 0$ for $\varepsilon \in W_0$, $(\tilde{x}, \tilde{y}) \in W$ if and only if $\tilde{x} = x(\varepsilon)$ and $\tilde{y} = y(\varepsilon)$. Moreover, $(x(0), y(0)) = -\sqrt[3]{\frac{\delta_{2k_0}}{2(1-\delta_{2k_0})(\mu_{k_0}^2+\nu_{k_0}^2)}}(\mu_{k_0}, \nu_{k_0})$. After scaling backwards in (37), we have the desired solution of equation (17). $\square$

Now, assume $\varepsilon h$ instead of $h$ in (17). Then we scale (36)

$$x \leftrightarrow \varepsilon x, \quad y \leftrightarrow \varepsilon y, \quad h_k \leftrightarrow \varepsilon h_k, \quad \gamma \leftrightarrow \gamma \quad (39)$$

and divide both equations by $\varepsilon^2$. This leads to

$$\begin{align}
\gamma \omega k_0 y + \mu_{k_0} + O(\varepsilon) &= 0 \nonumber \\
-\gamma \omega k_0 x + \nu_{k_0} + O(\varepsilon) &= 0 \tag{40}
\end{align}$$

and the following result.

**Theorem 3.** Assume (R), $\gamma \neq 0$ and $\varepsilon h$ instead of $h$ in (17) with $(\mu_{k_0}, \nu_{k_0}) \neq (0, 0)$. Then equation (17) has a solution $U \in Z$ close to 0 for any $\varepsilon \neq 0$ sufficiently small.
**Proof.** Let us denote

$$ \tilde{H}_2(x, y, \varepsilon) := (\gamma \omega k_0 y + \mu_{k_0} + O(\varepsilon), -\gamma \omega k_0 x + v_{k_0} + O(\varepsilon)). $$

Implicit function theorem applied to equation (40), i.e. $\tilde{H}_2(x, y, \varepsilon) = 0$ gives the existence of neighbourhoods $W_0 \subset V_0 \subset \mathbb{R}$ of the point 0 in $\mathbb{R}$, $W \subset \mathbb{R}^2$ of $(\frac{v_{k_0}}{\gamma \omega k_0}, -\frac{\mu_{k_0}}{\gamma \omega k_0})$ in $\mathbb{R}^2$ and a unique pair of continuous functions $x(\varepsilon)$, $y(\varepsilon)$ defined in $W_0$ such that $\tilde{H}_2(\tilde{x}, \tilde{y}, \varepsilon) = 0$ for $\varepsilon \in W_0$, $(\tilde{x}, \tilde{y}) \in W$ if and only if $\tilde{x} = x(\varepsilon)$ and $\tilde{y} = y(\varepsilon)$. Moreover, $x(0) = \frac{v_{k_0}}{\gamma \omega k_0}$, $y(0) = -\frac{\mu_{k_0}}{\gamma \omega k_0}$. After scaling backwards in (39), we have the desired solution of equation (17). $\square$

Finally, we deal with a gap between Theorems 2 and 3, namely that $h$ is not scaled and $\mu_{k_0} = v_{k_0} = 0$. Now, we scale (36)

$$ x \leftrightarrow \varepsilon x, \quad y \leftrightarrow \varepsilon y, \quad \varepsilon \leftrightarrow \pm\varepsilon^2, \quad \gamma \leftrightarrow \gamma \quad (41) $$

and divide both equations by $\varepsilon^3$. This leads to

$$ \begin{aligned} & \pm 2 \left( \frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) x \pm \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} + \gamma\omega k_0 \right) y + \frac{2(1 - \delta_{2k_0})x(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) = 0, \\ & \pm \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} - \gamma\omega k_0 \right) x \pm 2 \left( -\frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) y + \frac{2(1 - \delta_{2k_0})y(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) = 0 \end{aligned} \quad (42) $$

and the following result.

**Theorem 4.** Assume $(\mathbb{R})$ and $(\mu_{k_0}, v_{k_0}) = (0, 0)$. Then equation (17) has a solution $U \in \mathbb{Z}$ close to 0 for any $\varepsilon \neq 0$ sufficiently small provided

$$ \delta_{2k_0}^2 (4\mu_0^2 + \gamma^2 \omega^2 k_0^2) - 4(\mu_{2k_0}^2 + v_{2k_0}^2) \neq 0. \quad (43) $$

The order of this solution is $O(\varepsilon)$. If, in addition, it holds

$$ \delta_{2k_0}^2 (4\mu_0^2 + \gamma^2 \omega^2 k_0^2) - 4(\mu_{2k_0}^2 + v_{2k_0}^2) < 0, \quad (44) $$

then (17) has another solution. This one is of order $O(\sqrt{|\varepsilon|})$.

**Proof.** Let us denote

$$ \tilde{H}_3(x, y, \varepsilon) := \left( \pm 2 \left( \frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) x \pm \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} + \gamma\omega k_0 \right) y + \frac{2(1 - \delta_{2k_0})x(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon), \right. \\ \left. \pm \left( \frac{2\nu_{2k_0}}{\delta_{2k_0}} - \gamma\omega k_0 \right) x \pm 2 \left( -\frac{\mu_{2k_0}}{\delta_{2k_0}} + \mu_0 \right) y + \frac{2(1 - \delta_{2k_0})y(x^2 + y^2)}{\delta_{2k_0}} + O(\varepsilon) \right). $$

First we note that $\tilde{H}_3(x, y, \varepsilon) = A(x, y) + B(x, y) + O(\varepsilon)$ where $A$ is linear and $B$ is cubic. Since (43) is equivalent to $\det A \neq 0$, implicit function theorem applied to equation (42),
i.e. $\tilde{H}_3(x, y, \varepsilon) = 0$ gives the existence of neighbourhoods $W_0 \subset V_0 \subset \mathbb{R}$ of the point 0 in $\mathbb{R}$, $W \subset \mathbb{R}^2$ of $(0, 0)$ in $\mathbb{R}^2$ and a unique pair of continuous functions $x(\varepsilon)$, $y(\varepsilon)$ defined in $W_0$ such that $\tilde{H}_3(\tilde{x}, \tilde{y}, \varepsilon) = 0$ for $\varepsilon \in W_0$, $(\tilde{x}, \tilde{y}) \in W$ if and only if $\tilde{x} = x(\varepsilon)$ and $\tilde{y} = y(\varepsilon)$. Moreover, $(x(0), y(0)) = (0, 0)$, so $x(\varepsilon) = O(\varepsilon)$, $y(\varepsilon) = O(\varepsilon)$, which by scaling (41) gives $O(\varepsilon)$ solution for (17). On the other hand, if (44) holds then the Brouwer degree $\deg(\tilde{H}_3, B_1, (0, 0)) = \operatorname{sgn} \det A = -1$ while $\deg(\tilde{H}_3, B_2, (0, 0)) = 1$ for small and large balls $B_1, B_2$ in $\mathbb{R}^2$ centred at $(0, 0)$. This gives another solution of (42) of order $O(1)$ placed in the set $B_2 \setminus B_1$. After scaling backwards in (41), we have the desired $O(\sqrt{|\varepsilon|})$ solution of (17). $\square$

**Remark 2.** The above $O(\varepsilon)$ solution for (17) can be derived by putting Taylor expansions $x(\varepsilon) = x'(0)\varepsilon + O(\varepsilon^2)$, $y(\varepsilon) = y'(0)\varepsilon + O(\varepsilon^2)$ into (36) and comparing $\varepsilon^2$-terms. So we derive

$$ 
\begin{aligned}
& 2 \left( \frac{\mu_2 k_0}{\delta_{2k_0}} + \mu_0 \right) x'(0) + \left( \frac{2\nu_2 k_0}{\delta_{2k_0}} + \gamma \omega k_0 \right) y'(0) \\
&= \frac{2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_k \mu_{k-k_0} + \nu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}}, \\
& \left( \frac{2\nu_2 k_0}{\delta_{2k_0}} - \gamma \omega k_0 \right) x'(0) + 2 \left( -\frac{\mu_2 k_0}{\delta_{2k_0}} + \mu_0 \right) y'(0) \\
&= \frac{2}{k_0^2} \sum_{k \in \mathbb{Z} \setminus \{0, \pm k_0, 2k_0\}} \frac{k_0(k - k_0)(\mu_{k-k_0} \nu_k - \mu_k \nu_{k-k_0})}{\delta_k \delta_{k-k_0}},
\end{aligned} 
\quad (45) $$

which determines the first order approximation of this solution. Hence by (23) and (25), we get

$$ U(z) = U_1(z) + U_2(U_1, \varepsilon) = \varepsilon \left( c'(0) e^{k_0 t z} + \overline{c'(0)} e^{-k_0 t z} + \sum_{k \in \mathbb{Z} \setminus \{\pm k_0\}} d_k e^{k_0 t z} \right) + O(\varepsilon^2) \quad (46) $$

for $c'(0) = x'(0) + y'(0)t$.

**Remark 3.** 1. Now we specify the $O(\sqrt{|\varepsilon|})$ solution (solutions) of Theorem 4. Nontrivial roots of equation

$$ \tilde{H}_3(x, y, 0) = 0 \quad (47) $$

can be explicitly calculated as follows. Denoting

$$ 
\begin{aligned}
a &:= \pm 2 \left( \frac{\mu_2 k_0}{\delta_{2k_0}} + \mu_0 \right), & b &:= \pm \left( \frac{2\nu_2 k_0}{\delta_{2k_0}} + \gamma \omega k_0 \right), \\
c &:= \pm \left( \frac{2\nu_2 k_0}{\delta_{2k_0}} - \gamma \omega k_0 \right), & d &:= \pm 2 \left( -\frac{\mu_2 k_0}{\delta_{2k_0}} + \mu_0 \right), & \nu &:= \frac{2(1 - \delta_{2k_0})}{\delta_{2k_0}}
\end{aligned} 
\quad (48) $$

with the same fixed sign in all $a, b, c, d$ depending on the sign of unscaled $\varepsilon$, this equation is equivalent to the system
$$
\begin{align}
\tilde{a}x + \tilde{b}y + x(x^2 + y^2) &= 0 \nonumber \\
\tilde{c}x + \tilde{d}y + y(x^2 + y^2) &= 0 \tag{49}
\end{align}
$$

for $\tilde{a} = a/\upsilon$, $\tilde{b} = b/\upsilon$, $\tilde{c} = c/\upsilon$, $\tilde{d} = d/\upsilon$. Note $\upsilon \neq 0$, i.e., $\delta_{2k_0} \neq 1$. It is easy to see that $y = 0$
whenever $x = 0$ and vice versa. Hence we assume $x \neq 0 \neq y$. Dividing the first equation by $x$, the second equation by $y$, and comparing the resulting equations, one derives

$$
\tilde{a} + \tilde{b} \frac{y}{x} = \tilde{c} \frac{x}{y} + \tilde{d}.
$$

Therefrom $x = q_{1,2}y$ for

$$
q_{1,2} = \frac{(\tilde{a} - \tilde{d}) \pm \sqrt{(\tilde{a} - \tilde{d})^2 + 4\tilde{b}\tilde{c}}}{2\tilde{c}}
$$

(here $q_1$ corresponds to “+”-sign and $q_2$ to “-”-sign). As a consequence, putting $x$ in the second equation of (49) one obtains

$$
y_{1,2} = \pm \sqrt{-\frac{\tilde{c}q_{1,2} + \tilde{d}}{1 + q_{1,2}^2}}
$$

(again $y_1$ corresponds to “+”, $y_2$ to “-”), i.e., one has four points $(q_i y_j(q_i), y_j(q_i))$, $i, j = 1, 2$.
Condition (44) means $ad - bc < 0$, i.e., $\tilde{a}\tilde{d} - \tilde{b}\tilde{c} < 0$. So if (44) holds, then

$$
\sqrt{(\tilde{a} - \tilde{d})^2 + 4\tilde{b}\tilde{c}} = \sqrt{(\tilde{a} + \tilde{d})^2 - 4(\tilde{a}\tilde{d} - \tilde{b}\tilde{c})} > |\tilde{a} + \tilde{d}|
$$

Accordingly,

$$
\tilde{c}q_{1,2} + \tilde{d} = \frac{1}{2} \left( \tilde{a} + \tilde{d} \pm \sqrt{(\tilde{a} - \tilde{d})^2 + 4\tilde{b}\tilde{c}} \right) \quad (50)
$$

which is negative only for the minus sign (without any respect to the sign of unscaled $\epsilon$). Thus
the nontrivial solutions of (47) are

$$
(q_2 y_{1,2}(q_2), y_{1,2}(q_2)). \tag{51}
$$

On the other hand, assuming

$$
0 < \delta_{2k_0}^2 (4\mu_0^2 + \gamma^2 \omega^2 k_0^2) - 4(\mu_{2k_0}^2 + \nu_{2k_0}^2) < 4\mu_0^2, \quad (52)
$$

equation (47) still has nontrivial solutions. Indeed, in this case

$$
0 < \sqrt{(\tilde{a} - \tilde{d})^2 + 4\tilde{b}\tilde{c}} = \sqrt{(\tilde{a} + \tilde{d})^2 - 4(\tilde{a}\tilde{d} - \tilde{b}\tilde{c})} < |\tilde{a} + \tilde{d}|
$$
so $q_{1,2}$ are both real, again. However, (50) is negative for both signs on suppose that $\tilde{a} + \tilde{d} < 0$, i.e., $\pm \frac{\mu_0}{v} < 0$ where the sign depends on the sign of unscaled $\epsilon$, e.g. if originally $\epsilon > 0$ and $\frac{\mu_0}{v} < 0$, equation (47) has four distinct nontrivial solutions

$$ (q_i y_j(q_i), y_j(q_i)), \quad i, j = 1, 2, \qquad (53) $$

while for $\frac{\mu_0}{v} > 0$, (47) has no nontrivial solutions.

In conclusion, if condition (44) holds, there might exist two solutions of (17) of order $O(\sqrt{|\epsilon|})$ for $\epsilon \neq 0$ small, and if (52) is valid, there might exist four solutions of (17) of order $O(\sqrt{|\epsilon|})$ for $\epsilon \operatorname{sgn} \frac{\mu_0 \delta_{2k_0}}{1-\delta_{2k_0}} < 0$ and no solutions of order $O(\sqrt{|\epsilon|})$ for $\epsilon \operatorname{sgn} \frac{\mu_0 \delta_{2k_0}}{1-\delta_{2k_0}} > 0$. We recall by Theorem 4, that there exists a small solution in all cases of lower order than $O(\sqrt{|\epsilon|})$, namely of order $O(\epsilon)$. So we might have 3, 5 and 1 small solutions, respectively. To apply the implicit function theorem and so to confirm the existence of these solutions, one has to show that the derivative of $\tilde{H}_3$ at points (51) or (53) is surjective, which is rather awkward for general $a, b, c, d$, and so we do not go into details in this paper. But surely these solutions generically exist. Indeed, the determinant of the Jacobian of the left hand side of (49) at any $(qy, y)$ is

$$ -\tilde{b}\tilde{c} + \tilde{a}\tilde{d} + (3\tilde{a} + \tilde{d} + q(-2(\tilde{b} + \tilde{c}) + (\tilde{a} + 3\tilde{d})q)) y^2 + 3(1+q^2)^2 y^4. \qquad (54) $$

Now we fix all parameters only $\mu_0$ is variable. Then $q = q_{1,2}$, $\tilde{b}$, $\tilde{c}$ are constant (independent of $\mu_0$) and $y_{1,2}^2$, $\tilde{a}$, $\tilde{d}$ linearly depend on $\mu_0$: $\tilde{a} = \pm \frac{2}{v}(r_1 + \mu_0)$, $\tilde{d} = \pm \frac{2}{v}(-r_1 + \mu_0)$ and $y^2 = y_{1,2}^2 = \frac{-\tilde{c}vq \pm 2(r_1 - \mu_0)}{v(1+q^2)}$ for constants $r_1, r_2$ independent of $\mu_0$. Inserting these expressions into (54), we have

$$ \frac{1}{v^2(1+q^2)} (a_{\pm}(r_1, q, \tilde{b}, \tilde{c}, v) + b_{\pm}(r_1, q, \tilde{b}, \tilde{c}, v)\mu_0) \qquad (55) $$

for polynomials

$$ \begin{aligned} a_+(r_1, q, \tilde{b}, \tilde{c}, v) &= \tilde{c}(-\tilde{b} + (\tilde{b} + 5\tilde{c})q^2 + 3\tilde{c}q^4) v^2 - 4q(\tilde{b} + \tilde{c}(5+2q^2)) v r_1 + 16r_1^2, \\ b_+(r_1, q, \tilde{b}, \tilde{c}, v) &= 4(q(\tilde{b} + \tilde{c}(2+q^2)) v - 4r_1), \\ a_-(r_1, q, \tilde{b}, \tilde{c}, v) &= c(-\tilde{b} + (\tilde{b} + 5\tilde{c})q^2 + 3\tilde{c}q^4) v^2 + 4r_1(q(\tilde{b} + c(5+2q^2)) v + 4r_1), \\ b_-(r_1, q, \tilde{b}, \tilde{c}, v) &= -4(q(\tilde{b} + \tilde{c}(2+q^2)) v + 4r_1). \end{aligned} $$

Now it is clear that $a_{\pm}$, $b_{\pm}$ are generically nonzero and then (55) is nonzero for generic $\mu_0$. This shows that the above solutions exist in generic cases.

2. Let $p = \frac{\pi}{4}$, $\lambda = \frac{12}{37\sqrt{2}}$. Then by (13) we have

$$ M_2 = \left\{ \cos \frac{k\pi}{4} \mid k = \{0, 1, \dots, 7\} \right\} = \left\{ 0, \pm \frac{\sqrt{2}}{2}, \pm 1 \right\} \qquad (56) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. The graph of functions $\frac{37}{1225}$ and $\frac{1}{x^2(1-2\lambda \cos xp)}$ with values at nonzero integers depicted with asterisks.</figcaption></figure>

and condition (15) is fulfilled. If $\omega^2 = \frac{37}{25}$, then $k_0 = \pm 1$ satisfies the resonance condition (16). Note that

$$ \frac{1}{k^2(1 - 2\lambda \cos kp)} \le \frac{1}{k^2(1 - 2\lambda)} = \frac{37}{k^2(37 - 12\sqrt{2})} < \frac{37}{13k^2} < \frac{37}{25} $$

whenever $|k| \ge 2$. Thus $k_0 = \pm 1$ is the only resonant couple corresponding to $\omega^2 = \frac{37}{25}$. This is the case of a simple resonance discussed in this section.

On the other side, if $\omega^2 = \frac{37}{1225}$, then there are at least two couples $k_0 = \pm 5, \pm 7$ such that (16) is fulfilled. Since

$$ \frac{37}{k^2(37 - 12\sqrt{2})} < \frac{37}{20k^2} < \frac{37}{1225} $$

whenever $|k| \ge 8$, it is sufficient to evaluate $\frac{1}{k^2(1-2\lambda \cos kp)}$ at $k = \{1, \dots, 7\}$ to determine the exact number of resonant modes. From Fig. 1, one can see that for $\omega^2 = \frac{37}{1225}$ there are two resonant couples $k_0 = \pm 5, \pm 7$. This is a double resonance which can be investigated analogically to the present section.

## 4. Examples

In this section we present several examples to illustrate the above results.

**Example 1.** First we study the non-resonant case from Section 2. We consider the equation

$$ \omega^2 \left( U(z) - \lambda U \left( z - \frac{\pi}{4} \right) - \lambda U \left( z + \frac{\pi}{4} \right) \right)'' + \gamma \omega U'(z) + U(z) + \omega^2 \left( U^2(z) - \lambda U^2 \left( z - \frac{\pi}{4} \right) - \lambda U^2 \left( z + \frac{\pi}{4} \right) \right)'' + \gamma \omega (U^2(z))' - \varepsilon \cos z = 0 \quad (57) $$
for $z \in \mathbb{R}$ and parameters $\gamma, \lambda, \omega > 0$, $\frac{1}{2} \neq \lambda \neq \frac{\sqrt{2}}{2}$, and apply Theorem 1. Note that equation (57) is in the form of (2) with $p = \frac{\pi}{4}$, $h(z) = \varepsilon \cos z$, and $M_2$ of (13) is equivalent to (56), i.e., condition (15) is satisfied. Hence using Remark 1(3) and Theorem 1, we obtain the existence of a unique solution $U(\varepsilon, z) = O(\varepsilon)$ of equation (57) for any $\varepsilon$ sufficiently close to 0. Since $U(\varepsilon, z)$ is $C^\infty$-smooth, we can expand $U(\cdot, z)$ into Taylor series $U(\varepsilon, z) = \varepsilon \tilde{U}_1(z) + O(\varepsilon^2)$, and look for the function $\tilde{U}_1$. It is easy to see that $\tilde{U}_1$ is uniquely determined by the equation

$$ \omega^2 \left( \tilde{U}_1(z) - \lambda \tilde{U}_1 \left( z - \frac{\pi}{4} \right) - \lambda \tilde{U}_1 \left( z + \frac{\pi}{4} \right) \right)'' + \gamma \omega \tilde{U}_1'(z) + \tilde{U}_1(z) - \cos z = 0. $$

We look for its solution in the form of $\tilde{U}_1(z) = A \cos z + B \sin z$, and, by comparing coefficients at $\cos z$ and $\sin z$, we obtain

$$ -\gamma \omega A + (1 + \omega^2(\sqrt{2}\lambda - 1))B = 0, \quad (1 + \omega^2(\sqrt{2}\lambda - 1))A + \gamma \omega B = 1. $$

Therefrom,

$$ A = \frac{1 + \omega^2(\sqrt{2}\lambda - 1)}{\gamma^2\omega^2 + (1 + \omega^2(\sqrt{2}\lambda - 1))^2}, \quad B = \frac{\gamma\omega}{\gamma^2\omega^2 + (1 + \omega^2(\sqrt{2}\lambda - 1))^2}. $$

In particular, if

$$ \lambda = \frac{12}{37\sqrt{2}}, \quad \omega = \frac{\sqrt{37}}{5}, \qquad (58) $$

then $A = 0$, $B = \frac{5\sqrt{37}}{37\gamma}$, i.e.,

$$ U(\varepsilon, z) = \frac{5\sqrt{37}}{37\gamma}\varepsilon \sin z + O(\varepsilon^2). \qquad (59) $$

**Example 2.** Now, moving onto resonant cases represented by equation (17), we consider the equation

$$ \begin{aligned} & \omega^2 \left( U(z) - \lambda U \left( z - \frac{\pi}{4} \right) - \lambda U \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega U'(z) + U(z) \\ & + \omega^2 \left( U^2(z) - \lambda U^2 \left( z - \frac{\pi}{4} \right) - \lambda U^2 \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega (U^2(z))' - \varepsilon \cos z = 0, \end{aligned} \qquad (60) $$

i.e. when $\gamma$ in (57) is replaced by $\varepsilon\gamma$. Proceeding with point 2 of Remark 3 for parameters (58), then we obtain that $k_0 = \pm 1$, $\delta_{2k_0} = -\frac{123}{25}$ and $\nu = \frac{2(1-\delta_{2k_0})}{\delta_{2k_0}} = -\frac{296}{123}$.

Since $h(z) = \cos z = \frac{1}{2}e^{iz} + \frac{1}{2}e^{-iz}$ in (18), then $\mu_{k_0} = \frac{1}{2}$ and $\nu_{k_0} = 0$. So Theorem 2 can be applied to yield

$$ (x(0), y(0)) = -\sqrt[3]{\frac{\delta_{2k_0}}{2(1-\delta_{2k_0})(\mu_{k_0}^2 + \nu_{k_0}^2)}} (\mu_{k_0}, \nu_{k_0}) = \left(\sqrt[3]{\frac{123}{74}}, 0\right) \doteq (0.592281, 0). $$
Using (22), (25) and (37), we derive that the solution from Theorem 2 of (60) with parameters (58) is given by

$$ 
\begin{aligned}
& \sqrt[3]{\varepsilon c(\sqrt[3]{\varepsilon})} e^{i z} + \sqrt[3]{\varepsilon \bar{c}(\sqrt[3]{\varepsilon})} e^{-i z} + U_2 \left( \sqrt[3]{\varepsilon c(\sqrt[3]{\varepsilon})} e^{i z} + \sqrt[3]{\varepsilon \bar{c}(\sqrt[3]{\varepsilon})} e^{-i z} \right) \varepsilon \\
& = \sqrt[3]{\varepsilon} \left( c(0) e^{i z} + \bar{c}(0) e^{-i z} + O \left( \sqrt[3]{\varepsilon} \right) \right) = \sqrt[3]{\varepsilon} \left( 2x(0) \cos z + O \left( \sqrt[3]{\varepsilon} \right) \right) \\
& = \sqrt[3]{\frac{123 \varepsilon}{74}} \cos z + O \left( \sqrt[3]{\varepsilon^2} \right) \doteq \sqrt[3]{\varepsilon} 1.18456 \cos z + O \left( \sqrt[3]{\varepsilon^2} \right)
\end{aligned} 
\quad (61) $$

for a smooth function $c(\cdot)$ with $c(0) = x(0) + y(0)_t$ and $\varepsilon \neq 0$ small.

**Example 3.** Similarly, Theorem 3 can be applied in the equation

$$ 
\begin{aligned}
& \omega^2 \left( U(z) - \lambda U \left( z - \frac{\pi}{4} \right) - \lambda U \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega U'(z) + U(z) \\
& + \omega^2 \left( U^2(z) - \lambda U^2 \left( z - \frac{\pi}{4} \right) - \lambda U^2 \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega (U^2(z))' - \varepsilon^2 \cos z = 0,
\end{aligned} 
\quad (62) $$

with

$$ \lambda = \frac{12}{37\sqrt{2}}, \quad \omega = \frac{\sqrt{37}}{5}, \quad \gamma = 0.1. \quad (63) $$

Note that $\cos z$-term in (60) has been replaced by $\varepsilon \cos z$ in our case here, i.e., $h(z)$ in (18) is replaced by $\varepsilon h(z)$. Then using the proof of Theorem 3, we obtain that $(x(0), y(0)) = (0, -\frac{25}{\sqrt{37}}k_0)$, since $k_0 = \pm 1$. Using (22), (25) and (39), we derive that the solution from Theorem 3 of (62) with parameters (63) is given by

$$ 
\begin{aligned}
& \varepsilon c(\varepsilon) e^{i k_0 z} + \varepsilon \bar{c}(\varepsilon) e^{-i k_0 z} + U_2(\varepsilon c(\varepsilon) e^{i k_0 z} + \varepsilon \bar{c}(\varepsilon) e^{-i k_0 z}, \varepsilon) \\
& = \varepsilon(c(0) e^{i k_0 z} + \bar{c}(0) e^{-i k_0 z} + O(\varepsilon)) = \varepsilon(-2y(0) \sin k_0 z + O(\varepsilon)) \\
& = \frac{50\varepsilon}{\sqrt{37}} k_0 \sin k_0 z + O(\varepsilon^2) \doteq 8.21995\varepsilon \sin z + O(\varepsilon^2).
\end{aligned} 
\quad (64) $$

**Example 4.** Next, we consider

$$ 
\begin{aligned}
& \omega^2 \left( U(z) - \lambda U \left( z - \frac{\pi}{4} \right) - \lambda U \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega U'(z) + U(z) \\
& + \omega^2 \left( U^2(z) - \lambda U^2 \left( z - \frac{\pi}{4} \right) - \lambda U^2 \left( z + \frac{\pi}{4} \right) \right)'' + \varepsilon \gamma \omega (U^2(z))' \\
& - \varepsilon(\mu_0 + 2(\cos 2z - \sin 2z)) = 0
\end{aligned} 
\quad (65) $$

with parameters (63), i.e., equation (17) with the particular driving function $h(z) = \mu_0 + 2(\cos 2z - \sin 2z)$. So, $\mu_{k_0} = \nu_{k_0} = 0$ and $\mu_{2k_0} = \nu_2 = -\nu_{-2} = 1$. Hence, Theorems 2 and 3 cannot be used. However, using Theorem 4 one can obtain that (43) is of the form
$$-8 + \frac{15129}{625} \left( \frac{37}{2500} + 4\mu_0^2 \right) \neq 0,$$

i.e. $|\mu_0| \neq \frac{\sqrt{11940227}}{12300} \doteq 0.280932$. Moreover, (44) is equivalent to

$$|\mu_0| < \frac{\sqrt{11940227}}{12300}. \qquad (66)$$

First, let $k_0 = 1$ and consider $\varepsilon > 0$, i.e. use “+”-sign in (48). Then we can derive that $q_2 \doteq -0.294153$ and $y_{1,2} = \pm\sqrt{0.214884 + 0.764897\mu_0}$. Note that $0 < 0.214884 + 0.764897\mu_0 \le 0.429768$. Using formula (54), we can estimate the Jacobian of $\tilde{H}_3(q_2y_{1,2}(q_2), y_{1,2}(q_2), 0)$ as

$$0 < (0.218047 + 0.776154\mu_0)v^2 \leq 2.52553, \qquad (67)$$

since the formula (54) is equal to $\frac{1}{v^2} \det D\tilde{H}_3(qy, y, 0)$. Therefore, we have 3 different small solutions for $k_0 = 1$, $\varepsilon > 0$ small. Using $k_0 = -1$, we obtain $q_2 \doteq 0.294153$, $y_{1,2} = \pm\sqrt{0.214884 + 0.764897\mu_0}$ and by (54), estimate (67) follows. Hence for $k_0 = \pm 1$, $\varepsilon > 0$ small and $\mu_0$ satisfying (66), by (51) we obtain nontrivial solutions of (47). Scaling backwards in (41), we obtain the following small solutions of (65) with parameters (63)

$$\begin{aligned} & \sqrt{\varepsilon c}(\sqrt{\varepsilon})e^{tk_0z} + \sqrt{\varepsilon \tilde{c}}(\sqrt{\varepsilon})e^{-tk_0z} + U_2(\sqrt{\varepsilon c}(\sqrt{\varepsilon})e^{tk_0z} + \sqrt{\varepsilon \tilde{c}}(\sqrt{\varepsilon})e^{-tk_0z}, \varepsilon) \\ &= \sqrt{\varepsilon}(2x(0)\cos k_0z - 2y(0)\sin k_0z) + O(\varepsilon) = 2\sqrt{\varepsilon}y(0)(q_2\cos k_0z - \sin k_0z) + O(\varepsilon) \\ &= \pm 2\sqrt{\varepsilon}(0.214884 + 0.764897\mu_0)(0.294153\cos z + \sin z) + O(\varepsilon). \end{aligned} \qquad (68)$$

Moreover, there is another (smallest) solution of (65) with parameters (63) of order $U(z) = O(\varepsilon)$. Using (45), we obtain $x'(0) = y'(0) = 0$ and by (46), the solution is

$$U(\varepsilon) = \varepsilon \left( \mu_0 - \frac{50}{123} \cos 2z + \frac{50}{123} \sin 2z \right) + O(\varepsilon^2). \qquad (69)$$

Considering $k_0 = 1$ and $\varepsilon < 0$, i.e. using “-”-sign in (48), we derive $q_2 \doteq 1.83348$ and $y_{1,2} = \pm\sqrt{0.0535297 - 0.190543\mu_0}$. Note that $0 < 0.0535297 - 0.190543\mu_0 \le 0.107059$. Formula (54) yields

$$0 < 0.218047 - 0.776154\mu_0 \leq 0.436093, \qquad (70)$$

which implies that we also have 3 different small solutions for $k_0 = 1$, $\varepsilon < 0$ small and $\mu_0$ satisfying (66). Note that if $k_0 = -1$, $\varepsilon < 0$, we obtain $q_2 \doteq -1.83348$, $y_{1,2} = \pm\sqrt{0.0535297 - 0.190543\mu_0}$ and, by (54), estimate (70). The 3 solutions of (65) with parameters (63) have the forms (69) and

$$\begin{aligned} & \sqrt{-\varepsilon c}(\sqrt{-\varepsilon})e^{tk_0z} + \sqrt{-\varepsilon \tilde{c}}(\sqrt{-\varepsilon})e^{-tk_0z} + U_2(\sqrt{-\varepsilon c}(\sqrt{-\varepsilon})e^{tk_0z} + \sqrt{-\varepsilon \tilde{c}}(\sqrt{-\varepsilon})e^{-tk_0z}, \varepsilon) \\ &= \sqrt{-\varepsilon}(2x(0)\cos k_0z - 2y(0)\sin k_0z) + O(\varepsilon) = 2\sqrt{-\varepsilon}y(0)(q_2\cos k_0z - \sin k_0z) + O(\varepsilon) \\ &= \pm 2\sqrt{-\varepsilon}(0.0535297 - 0.190543\mu_0)(1.83348\cos z - \sin z) + O(\varepsilon). \end{aligned}$$
Next, (52) implies

$$0.280932 \doteq \frac{\sqrt{11\,940\,227}}{12\,300} < |\mu_0| < \frac{\sqrt{\frac{11\,940\,227}{74}}}{1400} \doteq 0.286921. \quad (71)$$

From the above computations, we can verify that formula (54) has one of the forms $0.218047 \pm 0.776154\mu_0 \neq 0$. Since sgn $v = -1$, we have that: if $\epsilon$ is small, $\epsilon \operatorname{sgn} \mu_0 > 0$ and $\mu_0$ fulfils (71), then there are 5 different small solutions $(4 \times O(\sqrt{|\epsilon|})$ and $1 \times O(\epsilon)$), if $\epsilon \operatorname{sgn} \mu_0 < 0$ there is only 1 small solution of order $O(\sqrt{|\epsilon|})$, actually (69). For instance, if $\epsilon > 0$ and $\mu_0$ satisfies (71), we obtain

$$ \begin{array}{c|l|l} k_0 & q_i & y_i \\ \hline 1 & q_1 = 1.83348 & y_{1,2}(q_1) = \pm\sqrt{-0.0535297 + 0.190543\mu_0} \\ & q_2 = -0.294153 & y_{1,2}(q_2) = \pm\sqrt{0.214884 + 0.764897\mu_0} \\ -1 & q_1 = -1.83348 & y_{1,2}(q_1) = \pm\sqrt{-0.0535297 + 0.190543\mu_0} \\ & q_2 = 0.294153 & y_{1,2}(q_2) = \pm\sqrt{0.214884 + 0.764897\mu_0} \end{array} $$

Hence, equation (65) with parameters (63) has the solutions (69) and

$$ \pm 2\sqrt{\epsilon(-0.0535297 + 0.190543\mu_0)}(1.83348 \cos z - \sin z) + O(\epsilon), \\ \pm 2\sqrt{\epsilon(0.214884 + 0.764897\mu_0)}(0.294153 \cos z + \sin z) + O(\epsilon) $$

for $\mu_0 > 0$, and the only solution (69) of order $O(\sqrt{\epsilon})$ for $\mu_0 < 0$.

## 5. Numerical results

To illustrate the theoretical results obtained in the previous sections, we have solved the governing equations (1) and (2) numerically. The advance-delay equation (2) is solved using a pseudo-spectral method, i.e. we express the solution $U$ in the Fourier series

$$U(z) = \sum_{j=1}^{J} \left[ B_j \cos\left((j-1)\tilde{k}z\right) + C_j \sin\left(j\tilde{k}z\right) \right], \quad (72)$$

where $\tilde{k} = 2\pi/L$ and $-L/2 < z < L/2$. Substituting the series in (2) and projecting it onto the Fourier modes in a reminiscence of the Galerkin truncation method, we obtain a system of algebraic equations for the Fourier coefficients $B_j$ and $C_j$, which are then solved using, e.g., a Newton-Raphson method. Typically, we use $L = 2\pi$ and large $J$ that will be indicated in each calculation below.

The stability of a solution obtained from (2) is determined numerically through calculating its Floquet multipliers, which are eigenvalues of the monodromy matrix. As a first-order system, the linearized equation of (1) is
<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) The oscillation amplitude of periodic solutions of (2) as a function of the driving amplitude $h_0$ for $\omega = 1.23$ with $\gamma = 0$ (black thick lines) and $\gamma = 0.001$ (blue thin lines). (b) Solution profiles at the driving amplitudes indicated in the legend for vanishing $\gamma$. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \dot{u}_n = v_n, \qquad (73) $$

$$ \Delta[\dot{v}_n(1+2U_n)] = -\gamma v_n(1+2U_n) - u_n(1+2\gamma\dot{U}_n) - 2\Delta[\ddot{U}_n u_n + 2\dot{U}_n v_n], \quad (74) $$

where $\Delta[\bullet_n] = \bullet_n - \lambda \bullet_{n-1} - \lambda \bullet_{n+1}$ and $U_n(t) = U(\omega t + np)$ (cf. (72)) is a periodic solution of (1). The linear system is integrated using a Runge-Kutta method of order four with periodic boundary conditions over the period $T = 2\pi/\omega$ and with the site index $n = 1, 2, \dots, N$. The $i$th column of the monodromy matrix $M$ is vector $[u_1(T), \dots, u_N(T), v_1(T), \dots, v_N(T)]^t$, that corresponds to the initial condition $[u_1(0), \dots, u_N(0), v_1(0), \dots, v_N(0)]^t$ equal to the $i$th column vector of the identity matrix $I_{2N}$. A periodic solution is asymptotically stable if all Floquet multipliers except the trivial Floquet multiplier are strictly smaller than one in modulus. The (in)stability result obtained from calculating the monodromy matrix is also confirmed by integrating (1). Note that the physically relevant interval for the coupling parameter is $|\lambda| < 1/2$ [6]. It is because there will be a band of unstable multipliers when $|\lambda| > 1/2$, implying that all solutions of the system will be modulationally unstable. Following Section 4, we set $\lambda = 12/(37\sqrt{2})$ and $p = \pi/4$. Computations for different parameter values have been performed, where we did not see any significant difference.

First, we consider the case of

$$ h(z) = h_0 \cos z. \qquad (75) $$

From (16) for $\gamma = 0$ the first resonance of the system is $k = \pm 1$, which corresponds to $\omega = \sqrt{37}/5$. In Fig. 2, we depict the bifurcation diagrams of periodic solutions for a value of $\omega = 1.23 > \sqrt{37}/5$. On the vertical axis is the oscillation amplitude $A$ of the periodic solutions defined as

$$ A = \frac{1}{2} (\max(U(z)) - \min(U(z))) $$

as a function of the driving amplitude $h_0$.
<figure><img src="image_4.png" /><figcaption>Fig. 3. The same as Fig. 2(a) for $\omega = \sqrt{37}/5$ with (a) $\gamma = 0$ and (b) $\gamma = 0.1$. In both panels, the solid black and the dashed blue lines are the numerics and the analytical results (61) (a) and (59), (64) (b), respectively.</figcaption></figure>

To check the convergence and dependence of the method on the number of modes, we computed the bifurcation diagrams for several values of $J$. Interestingly we observed that there is always a critical value of $h_0$ beyond which the bifurcation diagrams depend significantly on the parameter. The diagrams for three values of $J$ are shown in Fig. 2(a). Using the figure, we conclude that to determine the validity region of our numerical method, we need to compute bifurcation diagrams using at least two different values of $J$. The ‘breaking point’ of the method is when the curves start branching out. Studying the solution profiles near the breaking point as shown in Fig. 2(b), the breaking point seemingly corresponds to an extreme point of $U(z)$ having discontinuous first derivative. We therefore infer that our equation (2) may have peaked-periodic solutions.

As shown as black thick lines in Fig. 2(a), when $\gamma = 0$ there is one saddle-node bifurcation that occurs as we vary $h_0$. Note that the bifurcation diagram has a mirror symmetric (with respect to the vertical axis $h_0 = 0$) counterpart, i.e. the system is symmetric with the transformation $h_0 \to -h_0$, $U \to -U$. When $\gamma$ is turned on to a non-zero value, the bifurcation curve interacts with its symmetry and at the intersection point $h_0 = 0$ for non-vanishing $A$ an additional turning-point is formed, i.e. there is a curve-splitting-and-merging. In that case, as shown as blue thin lines in Fig. 2(b), there are now two saddle-node bifurcations.

When $\omega$ is decreased toward the resonance frequency, the critical drive of $h_0$ for the occurrence of the first saddle-node bifurcation, i.e. the bifurcation that also exists when $\gamma = 0$, decreases towards zero. In the limit $\omega = \sqrt{37}/5$, we present in Fig. 3(a) the bifurcation curve of periodic solutions at the resonance. The numerics (shown as solid black curves) are in good agreement with the analytical result (dashed line) in Example 2, i.e. (61), that when $\gamma = 0$ the slope of the existence curve is singular at $h_0 = 0$.

In panel (b) of the same figure, we depict the bifurcation curve of the periodic solutions for the same parameter values, but with $\gamma \neq 0$. In agreement with the analytical result in Examples 1 and 3, i.e. (59) and (64), adding dissipation will regularize the slope at the origin.

We have also considered the multiple resonance case discussed in Section 4. For the numerics, we take the drive

$$h(z) = h_0 (\mu_0 + 2(\cos(2z) - \sin(2z))), \quad (76)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. The same as Fig. 2(a), but for the periodic drive (76). The dashed lines are the analytical amplitudes from (68) and (69). Here, $\gamma = 0$.</figcaption></figure>

with $\mu_0 = 1/4$, which satisfies the condition (66). From Example 4, we know that in this resonance case, there are three small solutions given by (68) and (69). Shown in Fig. 4 are numerically computed bifurcation curves of two of the three periodic solutions.

In panel (a), we show only one bifurcation diagram of the solutions (68). It is because the amplitudes of the two solutions are almost the same. In panel (b), we depict the smaller solution with amplitudes that scale as $O(h_0)$. One can note that in both plots our analytical results from Lyapunov-Schmidt reduction agree with the numerics.

We have also determined the stability of the solutions obtained numerically and found that in the non-resonant condition, the periodic solutions are generally stable for small enough driving amplitude $h_0$.

Shown in Fig. 5 are the Floquet multipliers of some periodic solutions, calculated using $N = 200$. Panels (a)–(b) depict the multipliers of the periodic solutions corresponding to the bifurcation curve in Fig. 3(b) for $h_0 = 0.01$ and 0.04, respectively. As all the multipliers in the first panel are inside the unit circle, we can deduce that the solution is stable. However, as $h_0$ increases further, some multipliers leave the circle as shown in panel (b) and induce modulational instability. Panel (c) of the same figure corresponds to a solution along the existence curve in Fig. 3(a), i.e. with $h_0 = 0.01$ and $\gamma = 0$. It is clear that the solution is unstable. For this resonant case, we found that all solutions for any driving amplitude are modulationally unstable.

We have also simulated the dynamics of the unstable solutions by solving the governing equation (1). We observed that the typical dynamics of the instability are in the form of solution blow-up, similarly to that in [1].

As a conclusion, we found that the qualitative behaviour of the system (1) is similar to that of forced FPU lattices discussed in [8], despite the significant difference in the coupling terms as explained in Section 1.

## Acknowledgments

The work of MA has been co-financed from resources of the operational program “Educa-
tion and Lifelong Learning” of the European Social Fund and the National Strategic Refer-
ence Framework (NSRF) 2007–2013 within the framework of the Action State Scholarships
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Floquet multipliers of periodic solutions for some parameter values (see the text). The dashed curve is the unit circle, showing as guide to the eye.</figcaption></figure>

Foundation's (IKY) mobility grants programme for the short term training in recognized scientific/research centres abroad for candidate doctoral or postdoctoral researchers in Greek universities or research centres. MF is partially supported by Grants VEGA-MS 1/0071/14, VEGA-SAV 2/0029/13 and by the Slovak Research and Development Agency under the contract No. APVV-14-0378. MP is supported by Grant VEGA-SAV 2/0029/13. The work of VR has been co-financed by the European Union (European Social Fund ESF) and Greek national
funds through the Operational Program Education and Lifelong Learning of the National Strategic Reference Framework (NSRF) Research Funding Program D.534 MIS: 379337: THALES Investing in knowledge society through the European Social Fund. VR and HS acknowledge partial support from the London Mathematical Society through a visitor grant.

## References

[1] M. Agaoglou, V.M. Rothos, H. Susanto, G.P. Veldes, D.J. Frantzeskakis, Bifurcation results for travelling waves in nonlinear magnetic metamaterials, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 24 (2014) 1450147.

[2] M.S. Berger, Nonlinearity and Functional Analysis, Academic Press, New York, 1977.

[3] G. Birkhoff, S. Mac Lane, A Survey of Modern Algebra, A K Peters, Ltd., New York, 2008.

[4] C. Chicone, Ordinary Differential Equations with Applications, Texts in Applied Mathematics, vol. 34, Springer, New York, 2006.

[5] T.J. Cui, D. Smith, R. Liu (Eds.), Metamaterials: Theory, Design, and Applications, Springer, New York, 2009.

[6] J. Diblík, M. Fečkan, M. Pospíšil, V.M. Rothos, H. Susanto, Travelling waves in nonlinear magnetic metamaterials, in: R. Carretero, J. Cuevas, D. Frantzeskakis, N. Karachalios, P.G. Kevrekidis, F. Palmero (Eds.), Localized Excitations in Nonlinear Complex Systems, Nonlinear Systems and Complexity, Springer International Publishing, Switzerland, 2014, pp. 335–358.

[7] C. Enkrich, M. Wegener, S. Linden, S. Burger, L. Zschiedrich, F. Schmidt, J.F. Zhou, Th. Koschny, C.M. Soukoulis, Magnetic metamaterials at telecommunication and visible frequencies, Phys. Rev. Lett. 95 (2005) 203901.

[8] M. Fečkan, M. Pospíšil, V.M. Rothos, H. Susanto, Periodic travelling waves of forced FPU lattices, J. Dynam. Differential Equations 25 (2013) 795–820.

[9] M.J. Freire, R. Marqués, F. Medina, M.A.G. Laso, F. Martín, Planar magnetoinductive wave transducers: theory and applications, Appl. Phys. Lett. 85 (2004) 4439.

[10] G. Gallavotti, The Fermi–Pasta–Ulam Problem: A Status Report, Lecture Notes in Phys., vol. 728, Springer, Berlin, Heidelberg, 2008.

[11] A.N. Grigorenko, A.K. Geim, H.F. Gleeson, Y. Zhang, A.A. Firsov, I.Y. Khrushchev, J. Petrovic, Nanofabricated media with negative permeability at visible frequencies, Nature 438 (2005) 335–338.

[12] A. Ishikawa, T. Tanaka, S. Kawata, Negative magnetic permeability in the visible light region, Phys. Rev. Lett. 95 (2005) 237401.

[13] J.D. Jackson, Classical Electrodynamics, third ed., John Wiley and Sons, New Jersey, 1999.

[14] N. Katsarakis, G. Konstantinidis, A. Kostopoulos, R.S. Penciu, T.F. Gundogdu, M. Kafesaki, E.N. Economou, Th. Koschny, C.M. Soukoulis, Magnetic response of split-ring resonators in the far-infrared frequency regime, Optim. Lett. 30 (2005) 1348–1350.

[15] I. Kourakis, N. Lazarides, G.P. Tsironis, Self-focusing and envelope pulse generation in nonlinear magnetic metamaterials, Phys. Rev. E 75 (2007) 067601.

[16] N. Lazarides, M. Eleftheriou, G.P. Tsironis, Discrete breathers in nonlinear magnetic metamaterials, Phys. Rev. Lett. 97 (2006) 157406.

[17] N.M. Litchinitser, V.M. Shalaev, Optical metamaterials: invisibility in visible and nonlinearities in reverse, in: C. Denz, S. Flach, Yu.S. Kivshar (Eds.), Nonlinearities in Periodic Structures and Metamaterials, Springer, New York, 2010, pp. 217–240.

[18] S. Longhi, Gap solitons in metamaterials, Waves Random Complex Media 15 (2005) 119–126.

[19] R. Marques, F. Martin, M. Sorolla, Metamaterials with Negative Parameters. Theory, Design, and Microwave Applications, John Wiley and Sons, New Jersey, 2008.

[20] D. Schurig, J.J. Mock, B.J. Justice, S.A. Cummer, J.B. Pendry, A.F. Starr, D.R. Smith, Metamaterial electromagnetic cloak at microwave frequencies, Science 314 (2006) 977–980.

[21] E. Shamonina, V.A. Kalinin, K.H. Ringhofer, L. Solymar, Magnetoinductive waves in one, two, and three dimensions, J. Appl. Phys. 92 (2002) 6252.

[22] O. Sydoruk, O. Zhuromskyy, E. Shamonina, L. Solymar, Phonon-like dispersion curves of magnetoinductive waves, Appl. Phys. Lett. 87 (2005) 072501.

[23] R.R.A. Syms, E. Shamonina, L. Solymar, Positive and negative refraction of magnetoinductive waves in two dimensions, Eur. Phys. J. B 46 (2005) 301–308.

[24] N.L. Tsitsas, T.P. Horikis, Y. Shen, P.G. Kevrekidis, N. Whitaker, D.J. Frantzeskakis, Short pulse equations and localized structures in frequency band gaps of nonlinear metamaterials, Phys. Lett. A 374 (2010) 1384–1388.
[25] N.L. Tsitsas, N. Rompotis, I. Kourakis, P.G. Kevrekidis, D.J. Frantzeskakis, Higher-order effects and ultrashort solitons in left-handed metamaterials, Phys. Rev. E 79 (2009) 037601.

[26] G.P. Veldes, The nonlinear magneto-inductive lattice model, unpublished.

[27] T.J. Yen, W.J. Padilla, N. Fang, D.C. Vier, D.R. Smith, J.B. Pendry, D.N. Basov, X. Zhang, Terahertz magnetic response from artificial materials, Science 303 (2004) 1494–1496.

[28] J. Zhou, Th. Koschny, M. Kafesaki, E.N. Economou, J.B. Pendry, C.M. Soukoulis, Saturation of the magnetic response of split-ring resonators at optical frequencies, Phys. Rev. Lett. 95 (2005) 223902.