M Agaoglou

M Fečkan

M Pospíšil

V M Rothos
email: rothos@auth.gr

H Susanto
email: hsusanto@essex.ac.uk.

Travelling waves in nonlinear magneto-inductive lattices

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

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.

Introduction

In this work, we consider [START_REF] Veldes | The nonlinear magneto-inductive lattice model[END_REF]

d 2 dt 2 (u n -λu n-1 -λu n+1 ) + γ d dt u n + u n + d 2 dt 2 (u 2 n -λu 2 n-1 -λu 2 n+1 ) + γ d dt u 2 n -h(ωt + pn) = 0 ( 1 )
where γ ≥ 0, λ, ω > 0, p = 0 are parameters and h ∈ C(R, R) is such that

h(z) = k∈Z h k e kız , k∈Z |h k | < ∞.
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 [START_REF]Metamaterials: Theory, Design, and Applications[END_REF]. 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 [START_REF] Schurig | Metamaterial electromagnetic cloak at microwave frequencies[END_REF] (see also [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] 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 [START_REF] Ishikawa | Negative magnetic permeability in the visible light region[END_REF][START_REF] Zhou | Saturation of the magnetic response of split-ring resonators at optical frequencies[END_REF] and demonstrated experimentally in, e.g., [START_REF] Enkrich | Magnetic metamaterials at telecommunication and visible frequencies[END_REF][START_REF] Grigorenko | Nanofabricated media with negative permeability at visible frequencies[END_REF][START_REF] Katsarakis | Magnetic response of split-ring resonators in the far-infrared frequency regime[END_REF][START_REF] Yen | Terahertz magnetic response from artificial materials[END_REF]. 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 [START_REF] Freire | Planar magnetoinductive wave transducers: theory and applications[END_REF][START_REF] Shamonina | Magnetoinductive waves in one, two, and three dimensions[END_REF][START_REF] Sydoruk | Phonon-like dispersion curves of magnetoinductive waves[END_REF][START_REF] Syms | Positive and negative refraction of magnetoinductive waves in two dimensions[END_REF]. Magneto-inductive metamaterials consisting of periodic arrays of SRRs have been made in one-, two-and three dimensions [START_REF] Marques | Metamaterials with Negative Parameters. Theory, Design, and Microwave Applications[END_REF]. 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 [START_REF] Jackson | Classical Electrodynamics[END_REF].

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 [START_REF] Tsitsas | Short pulse equations and localized structures in frequency band gaps of nonlinear metamaterials[END_REF] and higher-order nonlinear Schrödinger equations [START_REF] Tsitsas | Higher-order effects and ultrashort solitons in left-handed metamaterials[END_REF]. In the second class, metamaterials are modelled by arrays of coupled oscillators, i.e. lattice equations, such as a nonlinear Klein-Gordon equation [START_REF] Longhi | Gap solitons in metamaterials[END_REF] and coupled Klein-Gordon equations [START_REF] Kourakis | Self-focusing and envelope pulse generation in nonlinear magnetic metamaterials[END_REF][START_REF] Lazarides | Discrete breathers in nonlinear magnetic metamaterials[END_REF]. The governing equation (1) falls into the second approach with γ representing the loss coefficient, λ 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 [START_REF] Veldes | The nonlinear magneto-inductive lattice model[END_REF]. Note that the nonlinearity is different from our previous work [START_REF] Agaoglou | Bifurcation results for travelling waves in nonlinear magnetic metamaterials[END_REF][START_REF] Diblík | Travelling waves in nonlinear magnetic metamaterials[END_REF], 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 γ , λ 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 [START_REF] Agaoglou | Bifurcation results for travelling waves in nonlinear magnetic metamaterials[END_REF] the nonlinearity in the coupling terms between the sites is akin to that in the Fermi-Pasta-Ulam lattices [START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem: A Status Report[END_REF]. 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 [START_REF] Fečkan | Periodic travelling waves of forced FPU lattices[END_REF]. 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.

Existence of small periodic solutions

Putting u n (t) = U(ωt + pn), z = ωt + pn in (1), we obtain

ω 2 (U (z) -λU (z -p) -λU (z + p)) + γ ωU (z) + U(z) + ω 2 (U 2 (z) -λU 2 (z -p) -λU 2 (z + p)) + γ ω(U 2 (z)) -h(z) = 0. (2)
We take Banach spaces

X := U ∈ C(R, R) | U(z) = k∈Z c k e kız , k∈Z |c k | < ∞ , Y := ⎧ ⎨ ⎩ U ∈ C 1 (R, R) | U(z) = k∈Z c k e kız , k∈Z\{0} |k||c k | < ∞ ⎫ ⎬ ⎭ , Z := ⎧ ⎨ ⎩ U ∈ C 2 (R, R) | U(z) = k∈Z c k e kız , k∈Z\{0} k 2 |c k | < ∞ ⎫ ⎬ ⎭
with the norms

U := k∈Z |c k |, U 1 := |c 0 | + k∈Z\{0} |k||c k |, U 2 := |c 0 | + k∈Z\{0} k 2 |c k |, respectively. It is easy to verify that Z Y X are compact embeddings and U ≤ U 1 ∀U ∈ Y , U 1 ≤ U 2 ∀U ∈ Z.
The following lemma is clear.

Lemma 1. If U 1 , U 2 ∈ X then U 1 U 2 ∈ X and U 1 U 2 ≤ U 1 U 2 .
By setting

KU := ω 2 U (z) -λU (z -p) -λU (z + p) + γ ωU (z) + U(z), F(U, h) := -ω 2 (U 2 (z)) -λ(U 2 (z -p)) -λ(U 2 (z + p)) -γ ω(U 2 (z)) + h(z), equation (2)
has the form

KU = F(U, h) .
We have the next result.

Lemma 2. Function

F : Z × X → X fulfils F(U, h) ≤ 2ω(2ω + 4λω + γ ) U 2 2 + h , (3)
F(U 1 , h) -F(U 2 , h) ≤ 2ω(2ω + 4λω + γ ) U 1 -U 2 2 ( U 1 2 + U 2 2 ), (4)
F(U, h 1 ) -F(U, h 2 ) ≤ h 1 -h 2 ( 5
)
for any U, U 1 , U 2 ∈ Z and h, h 1 , h 2 ∈ X.
Proof. Property ( 5) is obvious. Next, we use

U = U 1 ≤ 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) ≤ 4 U(z + c) 2 2 for c = -p, 0, p. Moreover, U(z ± p) = k∈Z |c k || e kıp | = k∈Z |c k | = U for U(z) = k∈Z c k e kız .
Similarly we derive

(U 2 (z)) = 2 U(z)U (z) ≤ 2 U 2 2 .
Hence, property (3) is obtained. Analogically, we get

(U 2 1 (z + c) -U 2 2 (z + c)) ≤ (U 1 (z + c) -U 2 (z + c)) (U 1 (z + c) + U 2 (z + c)) + 2 (U 1 (z + c) -U 2 (z + c)) (U 1 (z + c) + U 2 (z + c)) + (U 1 (z + c) -U 2 (z + c))(U 1 (z + c) + U 2 (z + c)) ≤ 4 U 1 -U 2 2 U 1 + U 2 2
for c = -p, 0, p, and

(U 2 1 (z) -U 2 2 (z)) ≤ (U 1 (z) -U 2 (z)) (U 1 (z) + U 2 (z)) + (U 1 (z) -U 2 (z))(U 1 (z) + U 2 (z)) ≤ 2 U 1 -U 2 2 U 1 + U 2 2 since a 2 -b 2 = (a -b)(a + b) for any a, b ∈ C. Now, (4) follows easily. 2 Next, if U ∈ Z with U(z) = k∈Z c k e kız then KU(z) = k∈Z 1 -ω 2 k 2 + 2λω 2 k 2 cos kp + ıγ ωk c k e kız ( 6
)
and so K ∈ L(Z, X) with

K L(Z,X) ≤ 1 + ω 2 (1 + 2λ) + γ ω. If := inf k∈Z\{0} ⎧ ⎨ ⎩ 1, 1 k 2 + ω 2 (2λ cos kp -1) 2 + γ 2 ω 2 k 2 ⎫ ⎬ ⎭ > 0 ( 7 )
is a constant depending on γ , λ, ω and p, then we also have

K -1 ∈ L(X, Z) ⊂ L(X, Y ) ⊂ L(X). So K -1 : X → Z is continuous such that K -1 L(X,Z) ≤ 1 . ( 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 < 2 8ω(2ω + 4λω + γ ) . ( 9
)
Then equation (2) has a unique solution

U(h) ∈ B(ρ h ) in a closed ball B(ρ h ) := {U ∈ Z | U 2 ≤ ρ h },
where

ρ h = - 2 -8ω(2ω + 4λω + γ ) h 4ω(2ω + 4λω + γ ) . ( 10
)
Moreover, U(h) can be approximated by an iteration process. Finally, it holds

U(h 1 ) -U(h 2 ) 2 ≤ h 1 -h 2 2 -8ω(2ω + 4λω + γ ) max{ h 1 , h 2 } ( 11
)
for any h 1 , h 2 ∈ X satisfying [START_REF] Freire | Planar magnetoinductive wave transducers: theory and applications[END_REF].

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

U = R(U, h) := K -1 F(U, h)
in Z. We already know that R : Z × X → Z is continuous and by (3), [START_REF] Fečkan | Periodic travelling waves of forced FPU lattices[END_REF] such that

R(U, h) 2 ≤ 1 2ω(2ω + 4λω + γ ) U 2 2 + h .
Next, if there is ρ h > 0 such that

A(ρ h ) := ρ h -2ω(2ω + 4λω + γ )ρ 2 h = h , (12)
then R(•, h) maps B(ρ h ) into itself. So it remains to study [START_REF] Ishikawa | Negative magnetic permeability in the visible light region[END_REF]. In order to find H > 0 -the largest right-hand side of [START_REF] Ishikawa | Negative magnetic permeability in the visible light region[END_REF] when this equation has a solution ρ H > 0, we need to solve A(r) = H together with DA(r) = 0. This implies

ρ H = 4ω(2ω + 4λω + γ )
and [START_REF] Freire | Planar magnetoinductive wave transducers: theory and applications[END_REF]. So assuming [START_REF] Enkrich | Magnetic metamaterials at telecommunication and visible frequencies[END_REF], [START_REF] Freire | Planar magnetoinductive wave transducers: theory and applications[END_REF], we know that (12) has a positive solution ρ h < ρ H . We take the smallest one which clearly has the form [START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem: A Status Report[END_REF]. So R(•, h) maps B(ρ h ) into itself and, moreover, by ( 4), ( 8)

R(U 1 , h) -R(U 2 , h) 2 ≤ 4ω(2ω + 4λω + γ ) U 1 -U 2 2 ρ h for any U 1 , U 2 ∈ B(ρ h ). Hence, R(•, h) is a contraction on B(ρ h ) with a contraction constant 4ω(2ω + 4λω + γ )ρ h < 4ω(2ω + 4λω + γ )ρ H = 1.
The proof of the existence and uniqueness is finished by the Banach fixed point theorem [START_REF] Berger | Nonlinearity and Functional Analysis[END_REF].

Next, let h 1 , h 2 ∈ X satisfy (9). Then U(h i ) ∈ B(ρ h i ) ⊂ B(ρ H ) for i = 1, 2 and H := max { h 1 , h 2 }.
Note H satisfies [START_REF] Freire | Planar magnetoinductive wave transducers: theory and applications[END_REF]. By (4), ( 5) and (8), we derive

U(h 1 ) -U(h 2 ) 2 = R(U (h 1 ), h 1 ) -R(U (h 2 ), h 2 ) 2 ≤ R(U (h 1 ), h 1 ) -R(U (h 2 ), h 1 ) 2 + R(U (h 2 ), h 1 ) -R(U (h 2 ), h 2 ) 2 ≤ 4ω(2ω + 4λω + γ ) U(h 1 ) -U(h 2 ) 2 ρ H + h 1 -h 2
which implies [START_REF] Grigorenko | Nanofabricated media with negative permeability at visible frequencies[END_REF]. 2

Remark 1. In the following special cases, (7) holds and we can replace with the corresponding i in the above considerations:

1. If p ∈ 2πZ, λ > 1 2 , then ≥ 1 := min{1, ω 2 (2λ -1)} > 0.
2. Let γ = 0. If p ∈ πQ, then we can write p = p 1 p 2 π for some p 1 ∈ Z, p 2 ∈ N, where p 1 , p 2 are relatively prime (their only common divisor is 1). Moreover,

M 1 := {cos kp | k ∈ Z\{0}} ⊂ cos kπ p 2 | k ∈ {0, 1, . . . , 2p 2 -1} =: M 2 . ( 13
)
Indeed, for each k ∈ Z there exist i, j ∈ Z such that 0 ≤ i ≤ 2p 2 -1 and kp 1 = 2jp 2 + i. Hence, for cos kp ∈ M 1 we have

cos kp = cos kp 1 π p 2 = cos iπ p 2 ∈ M 2 .
To find a better relationship between M 1 and M 2 , we need to solve cos

k 1 p 1 p 2 π = cos k 2 π p 2 for k 1 ∈ Z\{0} and k 2 ∈ {0, 1, . . . , 2p 2 -1}. This is equivalent to k 1 p 1 p 2 π = k 2 π p 2 + 2πl for l ∈ Z and k 1 = ± 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 [START_REF] Birkhoff | A Survey of Modern Algebra[END_REF] that there exist k, l ∈ Z such that kp 1 = 1 + 2 lp 2 . Consequently, kk 2 p 1 = k 2 + 2 lk 2 p 2 and M 1 = M 2 . If p 1 is even, then p 1 = 2 p1 and k 1 2 p1 = k 2 + 2lp 2 , so k 2 = 2 k2 and we have k 1 p1 = k2 + lp 2 . Clearly p1 and p 2 are relatively prime. Thus, there exist k, l ∈ Z such that

k p1 = 1 + lp 2 . Then k k2 p1 = k2 + l k2 p 2 . Hence M 1 = cos 2kπ p 2 | k ∈ {0, 1, . . . , p 2 -1} M 2 . Nevertheless, assuming ω 2 / ∈ 1 k 2 (1 -2λ cos kp) k∈Z\{0} ( 14
)
and

0 / ∈ 2λM 2 -1 (15)
for M 2 defined in [START_REF] Jackson | Classical Electrodynamics[END_REF], we get the existence of δ > 0 such that

1 k 2 + ω 2 (2λ cos kp -1) 2 + γ 2 ω 2 k 2 ≥ ω 2 1 ω 2 k 2 + 2λ cos kp -1 ≥ δ > 0
for each k ∈ Z\{0}. Thus, if p ∈ πQ and ( 14), [START_REF] Kourakis | Self-focusing and envelope pulse generation in nonlinear magnetic metamaterials[END_REF] are valid, then

≥ min 1, inf k∈Z\{0} ω 2 1 ω 2 k 2 + 2λ cos kp -1 =: 2 > 0.
Note that

1 k 2 (1-2λ cos kp) → 0 if |k| → ∞.
So if p ∈ πQ and (15) holds, there is at most a finite number of resonant modes k 0 ∈ Z\{0} determined by equation

ω 2 = 1 k 2 0 (1 -2λ cos k 0 p) . ( 16
)
More precisely, for given λ, ω there is not more than 2(p 2 + 1) resonant modes k 0 . Indeed, let {k 0j } J j =1 be an increasing sequence of positive resonant modes corresponding to the fixed λ, ω. Then {1 -2λ cos k 0j p} J j =1 has to be decreasing, since

k 2 0j (1 -2λ cos k 0j p) = 1 ω 2 = const. ∀j = 1, . . . , J,
i.e., {cos k 0j p} J j =1 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 kπ p 2 | k = p 2 , . . . , 2p 2 } which has p 2 + 1 elements. Finally, to each k 0j corresponds resonant mode -k 0j .

3. If γ > 0, condition [START_REF] Enkrich | Magnetic metamaterials at telecommunication and visible frequencies[END_REF] is satisfied even without the non-resonance condition [START_REF] Katsarakis | Magnetic response of split-ring resonators in the far-infrared frequency regime[END_REF], since in each resonant mode k 0 we have

1 k 2 0 + ω 2 (2λ cos k 0 p -1) 2 + γ 2 ω 2 k 2 0 = γ ω |k 0 | > 0
and there is only the finite number of resonant modes. For each non-resonant mode k it holds

1 k 2 + ω 2 (2λ cos kp -1) 2 + γ 2 ω 2 k 2 ≥ 1 k 2 + ω 2 (2λ cos kp -1) > 0.
Summarizing, if γ > 0, p ∈ πQ and ( 15) holds, condition ( 7) is fulfilled and

≥ min 1, min k∈M γ ω |k| , inf k∈Z\({0}∪M) ω 2 1 ω 2 k 2 + 2λ cos kp -1 =: 3 > 0
where M is the set of resonant modes k 0 .

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 ∈ πQ, condition ( 15) is valid, for some k 0 ∈ Z\{0} equation ( 16) holds and has the only solutions ±k 0 , and we consider the equation

A(U, ε) := KU -F(U ) -εG(U ) = 0 (17)
with

KU := ω 2 (U (z) -λU (z -p) -λU (z + p)) + U(z) F(U ) := -ω 2 ((U 2 (z)) -λ(U 2 (z -p)) -λ(U 2 (z + p)) ) G(U )(z) := -γ ω(U (z) + (U 2 (z)) ) + h(z) ( 18
)
for U and ε small. We shall solve equation ( 17) via Lyapunov-Schmidt reduction method [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF]. Clearly, A(0, 0) = 0 and D U A(0, 0) = K. For simplicity we denote κ := |k 0 | and R κ 1 := N K, R κ 2 := RK the null space and the range of the operator K ∈ L(Z, X), respectively. Equation ( 6) with γ = 0 together with (R) yields that if U ∈ Z with U(z) = k∈Z c k e kız , then

KU(z) = k∈Z\{±k 0 } (1 -ω 2 k 2 + 2λω 2 k 2 cos kp)c k e kız . Consequently, R κ 1 = c e κız + c e -κız | c ∈ C , R κ 2 = ⎧ ⎨ ⎩ U ∈ X | U(z) = k∈Z\{±k 0 } c k e kız ⎫ ⎬ ⎭ . Let Q : X → R κ 2 be the projection onto R κ 2 ⊂ X given by QU (z) = k∈Z\{±k 0 } c k e kız .
Now, we take the decomposition

U = U 1 + U 2 , U 1 ∈ R κ 1 , U 2 ∈ R κ 2 ∩ Z for U ∈ Z and decou- ple (17) to equations QA(U 1 + U 2 , ε) = 0, ( 19
)
(I -Q)A(U 1 + U 2 , ε) = 0. ( 20
)
Applying the implicit function theorem to the first equation implies the existence of neighbourhoods

V 1 × V 0 ⊂ R κ 1 × R of the point (0, 0) in R κ 1 × R, V 2 ⊂ R κ 2 ∩ Z of 0 in R κ 2 ∩ Z and a unique C ∞ -function U 2 : V 1 × V 0 → V 2 such that QA(U 1 + U 2 , ε) = 0 for (U 1 , ε) ∈ V 1 × V 0 and U 2 ∈ V 2 if and only if U 2 = U 2 (U 1 , ε). Moreover, U 2 (0, 0) = 0.
Thus, equation ( 20) has the form

(I -Q)A(U 1 + U 2 (U 1 , ε), ε) = 0 ( 21
)
for

U 1 ∈ V 1 and ε ∈ V 0 .
Differentiating equation ( 19) we obtain

D U 1 U 2 (0, 0) = 0, D ε U 2 (0, 0) = -K -1 QD ε A(0, 0) = K -1 QG(0). ( 22
)
Hence

U 2 (U 1 , ε) = K -1 QG(0)ε + O( (U 1 , ε) 2 ). ( 23
)
Expanding F and G of ( 18) into Taylor series using the last expansion of U 2 (U 1 , ε) yields

F(U 1 + U 2 (U 1 , ε)) = 1 2 D 2 F(0)(U 1 + D ε U 2 (0, 0)ε) 2 + O( (U 1 , ε) 3 ), G(U 1 + U 2 (U 1 , ε)) = G(0) + DG(0)(U 1 + D ε U 2 (0, 0)ε) + O( (U 1 , ε) 2 ). Note that D 2 F(0)V 2 = 2F (V ) and εO( (U 1 , ε) 2 ) = O( (U 1 , ε) 3
). Thus equation ( 21) is equivalent to

(I -Q)(F(U 1 + D ε U 2 (0, 0)ε) + εG(0) + εDG(0)(U 1 + D ε U 2 (0, 0)ε) + O( (U 1 , ε) 3 )) = 0
or, using D ε U 2 (0, 0) and F of ( 18),

(I -Q)(ω 2 (((U 1 (z) + K -1 QG(0)(z)ε) 2 ) -λ((U 1 (z -p) + K -1 QG(0)(z -p)ε) 2 ) -λ((U 1 (z + p) + K -1 QG(0)(z + p)ε) 2 ) ) -εG(0)(z) -εDG(0)(z)(U 1 + K -1 QG(0)ε) + O( (U 1 , ε) 3 )) = 0. ( 24
)
Note that

((U 1 (z) + K -1 QG(0)(z)ε) 2 ) = 2(U 1 (z) + (K -1 QG(0)(z)) ε) 2 + 2(U 1 (z) + K -1 QG(0)(z)ε)(U 1 (z) + (K -1 QG(0)(z)) ε) for any z ∈ R. Let U 1 (z) = c e k 0 ız + c e -k 0 ız ∈ R κ 1 and h(z) = k∈Z h k e kız ∈ X. Then U 1 (z) = c e k 0 ız + c e -k 0 ız ∈ R κ 1 , K -1 QG(0)(z) = k∈Z\{±k 0 } d k e kız ∈ R κ 2 , (K -1 QG(0)(z)) = k∈Z\{0,±k 0 } dk e kız ∈ R κ 2 ( 25
)
where

c = ık 0 c, d k = h k 1 -ω 2 k 2 + 2λω 2 k 2 cos kp , dk = ıkd k . Clearly, d -k = d k , d-k = dk for each k ∈ Z\{0, ±k 0 }.
It can be easily shown that

(I -Q)(U 1 (z)) 2 = 0 = (I -Q)(U 1 (z)U 1 (z)), (I -Q)(U 1 (z)(K -1 QG(0)(z)) ) = c d2k 0 e -k 0 ız + c d2k 0 e k 0 ız . ( 26
)
Next, setting d±k 0 = d0 = 0, Hence the coefficient by e k 0 ız in ((

K -1 QG(0)(z)) ) 2 is b := k+j =k 0 dk dj = k∈Z dk dk 0 -k = k∈Z dk dk-k 0 = k∈Z\{0,±k 0 ,2k 0 } dk dk-k 0 .
Therefore

(I -Q)(((K -1 QG(0)(z)) ) 2 ) = b e k 0 ız +b e -k 0 ız . ( 27
)
Analogically, denoting

dk := -k 2 d k , k ∈ Z\{0, ±k 0 }, c := -k 2 0 c, b := k∈Z\{0,±k 0 ,2k 0 } d k dk-k 0 we have U 1 (z) = c e k 0 ız + c e -k 0 ız ∈ R κ 1 , (K -1 QG(0)(z)) = k∈Z\{0,±k 0 } dk e kız ∈ R κ 2 ,
and the next identities can be proved

(I -Q)(U 1 (z)(K -1 QG(0)(z)) ) = c d2k 0 e k 0 ız +c d2k 0 e -k 0 ız , (I -Q)(U 1 (z)K -1 QG(0)(z)) = cd 2k 0 + cd 0 e k 0 ız + cd 2k 0 + cd 0 e -k 0 ız , (I -Q)(K -1 QG(0)(z)(K -1 QG(0)(z)) ) = b e k 0 ız + b e -k 0 ız . ( 28
)
Moreover, by definition of G in ( 18),

(I -Q)G(0)(z) = h k 0 e k 0 ız +h k 0 e -k 0 ız , (I -Q)DG(0)(U 1 + K -1 QG(0)(z)ε) = -γ ω(c e k 0 ız + c e -k 0 ız ). ( 29
)
Now, we put ( 26), ( 27), ( 28), ( 29) in [START_REF] Tsitsas | Short pulse equations and localized structures in frequency band gaps of nonlinear metamaterials[END_REF] and collect coefficients by e k 0 ız and e -k 0 ız to obtain

H e k 0 ız +H e -k 0 ız +O( (U 1 , ε) 3 ) = 0 (30)
where

H = 2ω 2 2 c d2k 0 + c d2k 0 + cd 2k 0 + cd 0 (1 -λ e -k 0 ıp -λ e k 0 ıp )ε -(h k 0 -γ ω c)ε + 2ω 2 (b + b)(1 -λ e -k 0 ıp -λ e k 0 ıp )ε 2 .
Equation ( 30) is equivalent to

H + O( (c, c, ε) 3 ) = 0. ( 31
)
Using resonance condition ( 16), we get

(1 -λ e -k 0 ıp -λ e k 0 ıp ) = 1 -2λ cos k 0 p = 1 ω 2 k 2 0 .
Dropping tildes yields

2 c d2k 0 + c d2k 0 + cd 2k 0 + cd 0 = -k 2 0 cd 2k 0 -k 2 0 cd 0 , b + b = k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )d k d k-k 0 .
Thus we have to solve

-2(cd 2k 0 + cd 0 ) + h k 0 -γ ω c ε + 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )d k d k-k 0 + O( (c, c, ε) 3 ) = 0. Since c, h k ∈ C, we have c = x + ıy, h k = μ k + ıν k , d k = μ k + ıν k δ k where x, y, μ k , ν k ∈ R and δ k = 1 -ω 2 k 2 + 2λω 2 k 2 cos kp ∈ R. Obviously, δ 0 = 1, h 0 = d 0 = μ 0 ∈ R and μ k = μ -k , ν k = -ν -k ∀k ∈ N. So, we get - 2(x -ıy)(μ 2k 0 + ıν 2k 0 ) δ 2k 0 + 2(x + ıy)μ 0 + μ k 0 + ıν k 0 -γ ωık 0 (x + ıy) ε + 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k + ıν k )(μ k-k 0 -ıν k-k 0 ) δ k δ k-k 0 + O( (x, y, ε) 3 ) = 0
or, separating the real and imaginary parts,

- 2(xμ 2k 0 + yν 2k 0 ) δ 2k 0 + 2xμ 0 + μ k 0 + γ ωk 0 y ε + 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k μ k-k 0 + ν k ν k-k 0 ) δ k δ k-k 0 + O( (x, y, ε) 3 ) = 0, - 2(xν 2k 0 -yμ 2k 0 ) δ 2k 0 + 2yμ 0 + ν k 0 -γ ωk 0 x ε + 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k-k 0 ν k -μ k ν k-k 0 ) δ k δ k-k 0 + O( (x, y, ε) 3 ) = 0.
Multiplying both equations by (-1) and collecting terms by x and y we obtain

2 μ 2k 0 δ 2k 0 + μ 0 εx + 2ν 2k 0 δ 2k 0 + γ ωk 0 εy + μ k 0 ε - 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k μ k-k 0 + ν k ν k-k 0 ) δ k δ k-k 0 + O( (x, y, ε) 3 ) = 0, 2ν 2k 0 δ 2k 0 -γ ωk 0 εx + 2 - μ 2k 0 δ 2k 0 + μ 0 εy + ν k 0 ε - 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k-k 0 ν k -μ k ν k-k 0 ) δ k δ k-k 0 + O( (x, y, ε) 3 ) = 0. ( 32
)
Now we intend to derive the cubic O( (U 1 , 0) 3 ) of (30) as follows. We denote by B := -1 2 D 2 F(0) and u 2 (U 1 ) := U 2 (U 1 , 0). Then [START_REF] Marques | Metamaterials with Negative Parameters. Theory, Design, and Microwave Applications[END_REF] and ( 20) are equivalent to

Q(Ku 2 (U 1 ) + B(U 1 + u 2 (U 1 )) 2 ) = 0, ( 33
) (I -Q)B(U 1 + u 2 (U 1 )) 2 = 0. (34) First note that u 2 (U 1 ) = O( U 1 2
) by ( 23), (I -Q)BU 2 1 = 0 and then (34) has the form

0 = (I -Q)B(U 1 + u 2 (U 1 )) 2 = (I -Q)BU 2 1 + 2(I -Q)B(U 1 , u 2 (U 1 )) + (I -Q)B(u 2 (U 1 )) 2 = (I -Q)B(U 1 , D 2 U 1 u 2 (0)U 2 1 ) + O( U 1 4 ), hence O( (U 1 , 0) 3 ) = (I -Q)B(U 1 , D 2 U 1 u 2 (0)U 2 1 ). To derive D 2 U 1 u 2 (0), we twice differenti- ate (33) to obtain KD 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 KD 2 U 1 u 2 (0)(V 1 , V 2 ) + 2QB(V 1 , V 2 ) = 0, hence D 2 U 1 u 2 (0)U 2 1 = -2K -1 QBU 2 1 = 2K -1 QF(U 1 ),
consequently, we arrive at

O( (U 1 , 0) 3 ) = (I -Q)B(U 1 , D 2 U 1 u 2 (0)U 2 1 ) = -(I -Q)D 2 F(0)(U 1 , K -1 QF(U 1 )). ( 35
) Let U 1 (z) = c e k 0 ız + c e -k 0 ız ∈ R κ 1 . Then K -1 QF(U 1 ) = K -1 Q 4ω 2 k 2 0 (1 -2λ cos 2k 0 p) c 2 e 2k 0 ız + c2 e -2k 0 ız = 1 -δ 2k 0 δ 2k 0 c 2 e 2k 0 ız + 1 -δ 2k 0 δ 2k 0 c2 e -2k 0 ız .
Next we derive

D 2 F(0)(V , W ) = -2ω 2 ((V (z)W (z)) -λ(V (z -p)W (z -p)) -λ(V (z + p)W (z + p)) ) = -2ω 2 V (z)W (z) + V (z)W (z) + 2V (z)W (z) -λ V (z -p)W (z -p) + V (z -p)W (z -p) + 2V (z)W (z) -λ V (z + p)W (z + p) + V (z + p)W (z + p) + 2V (z)W (z) . Denoting μ = 1-δ 2k 0 δ 2k 0 we get U 1 (z)[K -1 QF(U 1 )(z)] = -4k 2 0 U 1 (z)[K -1 QF(U 1 )(z)]
= -4k 2 0 μc 3 e 3k 0 ız +μ c3 e -3k 0 ız +μc 2 c e k 0 ız +μc c2 e -k 0 ız ,

U 1 (z)[K -1 QF(U 1 )(z)] = -k 2 0 U 1 (z)[K -1 QF(U 1 )(z)]
= -k 2 0 μc 3 e 3k 0 ız +μ c3 e -3k 0 ız +μc 2 c e k 0 ız +μc c2 e -k 0 ız , U 1 (z)[K -1 QF(U 1 )(z)] = k 0 ıc e k 0 ız -k 0 ı c e -k 0 ız 2k 0 ıμc 2 e 2k 0 ız -2k 0 ıμ c2 e -2k 0 ız = -2k 2 0 μc 3 e 3k 0 ız -2k 2 0 μ c3 e -3k 0 ız +2k 2 0 μc 2 c e k 0 ız +2k 2 0 μc c2 e -k 0 ız .

So (I -Q) U 1 (z)[K -1 QF(U 1 )(z)] + 2U 1 (z)[K -1 QF(U 1 )(z)] = 0.
Hence

-(I -Q)D 2 F(0)(U 1 , K -1 QF(U 1 )) = 2ω 2 μc 2 c e k 0 ız (-k 2 0 ) + μc c2 e -k 0 ız (-k 2 0 ) -λμc 2 c e k 0 ı(z-p) (-k 2 0 ) -λμc c2 e -k 0 ı(z-p) (-k 2 0 ) -λμc 2 c e k 0 ı(z+p) (-k 2 0 ) -λμc c2 e -k 0 ı(z+p) (-k 2 0 ) = -2k 2 0 ω 2 μ|c| 2 U 1 (z) -λU 1 (z -p) -λU 1 (z + p) = -2k 2 0 ω 2 μ|c| 2 (1 -λ e k 0 ıp -λ e -k 0 ıp )U 1 (z) = -2μ|c| 2 U 1 (z) = - 2(1 -δ 2k 0 ) δ 2k 0 |c| 2 c e k 0 ız + c e -k 0 ız .
Thus, by (35) and the last identity, in (31) we have

O( (c, c, 0) 3 ) = - 2c|c| 2 (1 -δ 2k 0 ) δ 2k 0 .
Summarizing, the bifurcation equation (32) has a form

2 μ 2k 0 δ 2k 0 + μ 0 εx + 2ν 2k 0 δ 2k 0 + γ ωk 0 εy + μ k 0 ε + 2(1 -δ 2k 0 )x(x 2 + y 2 ) δ 2k 0 - 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k μ k-k 0 + ν k ν k-k 0 ) δ k δ k-k 0 + εO( (x, y, ε) 2 ) = 0, 2ν 2k 0 δ 2k 0 -γ ωk 0 εx + 2 - μ 2k 0 δ 2k 0 + μ 0 εy + ν k 0 ε + 2(1 -δ 2k 0 )y(x 2 + y 2 ) δ 2k 0 - 2ε 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k-k 0 ν k -μ k ν k-k 0 ) δ k δ k-k 0 + εO( (x, y, ε) 2 ) = 0. ( 36
)
Now, we scale (36)

x ←→ εx, y ←→ εy, ε ←→ ε 3 , γ ←→ γ (37)
and divide both equations by ε 3 . This leads to

μ k 0 + 2(1 -δ 2k 0 )x(x 2 + y 2 ) δ 2k 0 + O(ε) = 0 ν k 0 + 2(1 -δ 2k 0 )y(x 2 + y 2 ) δ 2k 0 + O(ε) = 0 ( 38
)
and the following result holds.

Theorem 2. Assume (R) and (μ k 0 , ν k 0 ) = (0, 0). Then equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] has a solution U ∈ Z close to 0 for any ε = 0 sufficiently small.

Proof. Let us denote

H 1 (x, y, ε) := μ k 0 + 2(1 -δ 2k 0 )x(x 2 + y 2 ) δ 2k 0 + O(ε), ν k 0 + 2(1 -δ 2k 0 )y(x 2 + y 2 ) δ 2k 0 + O(ε) .
Implicit function theorem applied to equation (38), i.e. H 1 (x, y, ε) = 0 gives the existence of

neighbourhoods W 0 ⊂ V 0 ⊂ R of the point 0 in R, W ⊂ R 2 of -3 δ 2k 0 2(1-δ 2k 0 )(μ 2 k 0 +ν 2 k 0 ) (μ k 0 , ν k 0 ) in R 2
and a unique pair of continuous functions x(ε), y(ε) defined in W 0 such that H 1 ( x, ỹ, ε) = 0 for ε ∈ W 0 , ( x, ỹ) ∈ W if and only if x = x(ε) and ỹ = y(ε). Moreover, (x(0), y(0

)) = -3 δ 2k 0 2(1-δ 2k 0 )(μ 2 k 0 +ν 2 k 0 ) (μ k 0 , ν k 0 ).
After scaling backwards in (37), we have the desired solution of equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF]. 2 Now, assume εh instead of h in [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF]. Then we scale (36)

x ←→ εx, y ←→ εy, h k ←→ εh k , γ ←→ γ ( 39
)
and divide both equations by ε 2 . This leads to

γ ωk 0 y + μ k 0 + O(ε) = 0 -γ ωk 0 x + ν k 0 + O(ε) = 0 ( 40
)
and the following result.

Theorem 3. Assume (R), γ = 0 and εh instead of h in [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] with (μ k 0 , ν k 0 ) = (0, 0). Then equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] has a solution U ∈ Z close to 0 for any ε = 0 sufficiently small.

Proof. Let us denote

H 2 (x, y, ε) := γ ωk 0 y + μ k 0 + O(ε), -γ ωk 0 x + ν k 0 + O(ε) .
Implicit function theorem applied to equation (40), i.e. H 2 (x, y, ε) = 0 gives the existence of neighbourhoods

W 0 ⊂ V 0 ⊂ R of the point 0 in R, W ⊂ R 2 of ( ν k 0 γ ωk 0 , - μ k 0 γ ωk 0 ) in R 2
and a unique pair of continuous functions x(ε), y(ε) defined in W 0 such that H 2 ( x, ỹ, ε) = 0 for ε ∈ W 0 , ( x, ỹ) ∈ W if and only if x = x(ε) and ỹ = y(ε). Moreover, x(0) = ν k 0 γ ωk 0 , y(0) = -μ k 0 γ ωk 0 . After scaling backwards in (39), we have the desired solution of equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF]. 2

Finally, we deal with a gap between Theorems 2 and 3, namely that h is not scaled and μ k 0 = ν k 0 = 0. Now, we scale (36)

x ←→ εx, y ←→ εy, ε ←→ ±ε 2 , γ ←→ γ (41)
and divide both equations by ε 3 . This leads to

±2 μ 2k 0 δ 2k 0 + μ 0 x ± 2ν 2k 0 δ 2k 0 + γ ωk 0 y + 2(1 -δ 2k 0 )x(x 2 + y 2 ) δ 2k 0 + O(ε) = 0, ± 2ν 2k 0 δ 2k 0 -γ ωk 0 x ± 2 - μ 2k 0 δ 2k 0 + μ 0 y + 2(1 -δ 2k 0 )y(x 2 + y 2 ) δ 2k 0 + O(ε) = 0 (42)
and the following result.

Theorem 4. Assume (R) and (μ k 0 , ν k 0 ) = (0, 0). Then equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] has a solution U ∈ Z close to 0 for any ε = 0 sufficiently small provided

δ 2 2k 0 (4μ 2 0 + γ 2 ω 2 k 2 0 ) -4(μ 2 2k 0 + ν 2 2k 0 ) = 0. ( 43
)
The order of this solution is O(ε). If, in addition, it holds

δ 2 2k 0 (4μ 2 0 + γ 2 ω 2 k 2 0 ) -4(μ 2 2k 0 + ν 2 2k 0 ) < 0, ( 44
)
then [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] has another solution. This one is of order O( √ |ε|).

Proof. Let us denote

H 3 (x, y, ε) := ± 2 μ 2k 0 δ 2k 0 + μ 0 x ± 2ν 2k 0 δ 2k 0 + γ ωk 0 y + 2(1 -δ 2k 0 )x(x 2 + y 2 ) δ 2k 0 + O(ε), ± 2ν 2k 0 δ 2k 0 -γ ωk 0 x ± 2 - μ 2k 0 δ 2k 0 + μ 0 y + 2(1 -δ 2k 0 )y(x 2 + y 2 ) δ 2k 0 + O(ε) .
First we note that H 3 (x, y, ε) = A(x, y) + B(x, y) + O(ε) where A is linear and B is cubic. Since (43) is equivalent to det A = 0, implicit function theorem applied to equation (42),

i.e. H 3 (x, y, ε) = 0 gives the existence of neighbourhoods W 0 ⊂ V 0 ⊂ R of the point 0 in R, W ⊂ R 2 of (0, 0) in R 2 and a unique pair of continuous functions x(ε), y(ε) defined in W 0 such that H 3 ( x, ỹ, ε) = 0 for ε ∈ W 0 , ( x, ỹ) ∈ W if and only if x = x(ε) and ỹ = y(ε). Moreover, (x(0), y(0)) = (0, 0), so 36) and comparing ε 2 -terms. So we derive

x(ε) = O(ε), y(ε) = O(ε),
(ε) = x (0)ε + O(ε 2 ), y(ε) = y (0)ε + O(ε 2 ) into (
2 μ 2k 0 δ 2k 0 + μ 0 x (0) + 2ν 2k 0 δ 2k 0 + γ ωk 0 y (0) = 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k μ k-k 0 + ν k ν k-k 0 ) δ k δ k-k 0 , 2ν 2k 0 δ 2k 0 -γ ωk 0 x (0) + 2 - μ 2k 0 δ 2k 0 + μ 0 y (0) = 2 k 2 0 k∈Z\{0,±k 0 ,2k 0 } k 0 (k -k 0 )(μ k-k 0 ν k -μ k ν k-k 0 ) δ k δ k-k 0 , (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 , ε) = ε ⎛ ⎝ c (0) e k 0 ız +c (0) e -k 0 ız + k∈Z\{±k 0 } d k e kız ⎞ ⎠ + O(ε 2 ) (46)
for c (0) = x (0) + y (0)ı. Therefrom x = q 1,2 y for

a := ±2 μ 2k 0 δ 2k 0 + μ 0 , b:= ± 2ν 2k 0 δ 2k 0 + γ ωk 0 , c := ± 2ν 2k 0 δ 2k 0 -γ ωk 0 , d:= ±2 - μ 2k 0 δ 2k 0 + μ 0 , υ:= 2(1 -δ 2k 0 ) δ 2k 0 ( 48
q 1,2 = ( a -d) ± ( a -d) 2 + 4 b c 2 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 = ± - cq 1,2 + d 1 + q 2 1,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 adbc < 0, i.e., a db c < 0. So if (44) holds, then

( a -d) 2 + 4 b c = ( a + d) 2 -4( a d -b c) > | a + d|.
Accordingly,

cq 1,2 + d = 1 2 a + d ± ( a -d) 2 + 4 b c (50)
which is negative only for the minus sign (without any respect to the sign of unscaled ε). Thus the nontrivial solutions of (47) are

(q 2 y 1,2 (q 2 ), y 1,2 (q 2 )). ( 51
)
On the other hand, assuming

0 < δ 2 2k 0 (4μ 2 0 + γ 2 ω 2 k 2 0 ) -4(μ 2 2k 0 + ν 2 2k 0 ) < 4μ 2 0 , ( 52
)
equation ( 47) still has nontrivial solutions. Indeed, in this case

0 < ( a -d) 2 + 4 b c = ( a + d) 2 -4( a d -b c) < | a + d|,
so q 1,2 are both real, again. However, (50) is negative for both signs on suppose that a + d < 0, i.e., ± μ 0 υ < 0 where the sign depends on the sign of unscaled ε, e.g. if originally ε > 0 and μ 0 υ < 0, equation (47) has four distinct nontrivial solutions (q i y j (q i ), y j (q i )), i, j = 1, 2, (53) while for μ 0 υ > 0, (47) has no nontrivial solutions. In conclusion, if condition (44) holds, there might exist two solutions of ( 17) of order O( √ |ε|)

for ε = 0 small, and if (52) is valid, there might exist four solutions of ( 17) of order O( √ |ε|) for ε sgn

μ 0 δ 2k 0 1-δ 2k 0
< 0 and no solutions of order O( √ |ε|) for ε sgn

μ 0 δ 2k 0 1-δ 2k 0
> 0. We recall by Theorem 4, that there exists a small solution in all cases of lower order than O( √ |ε|), namely of order O(ε). 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 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

-b c + a d + 3 a + d + q(-2( b + c) + ( a + 3 d)q) y 2 + 3 1 + q 2 2 y 4 . ( 54
)
Now we fix all parameters only μ 0 is variable. Then q = q 1,2 , b, c are constant (independent of μ 0 ) and y 2 1,2 , a, d linearly depend on μ 0 :

a = ± 2 υ (r 1 + μ 0 ), d = ± 2 υ (-r 1 + μ 0 ) and y 2 = y 2 1,2 = -cυq±2(r 1 -μ 0 ) υ(1+q 2 )
for constants r 1 , r 2 independent of μ 0 . Inserting these expressions into (54), we have

1 υ 2 (1 + q 2 ) a ± (r 1 , q, b, c, υ) + b ± (r 1 , q, b, c, υ)μ 0 ( 55
)
for polynomials

a + (r 1 , q, b, c, υ) = c -b + ( b + 5 c)q 2 + 3 cq 4 υ 2 -4q b + c 5 + 2q 2 υr 1 + 16r 2 1 , b + (r 1 , q, b, c, υ) = 4 q b + c 2 + q 2 υ -4r 1 , a -(r 1 , q, b, c, υ) = c -b + ( b + 5 c)q 2 + 3 cq 4 υ 2 + 4r 1 q b + c 5 + 2q 2 υ + 4r 1 , b -(r 1 , q, b, c, υ) = -4 q b + c 2 + q 2 υ + 4r 1 .
Now it is clear that a ± , b ± are generically nonzero and then (55) is nonzero for generic μ 0 . This shows that the above solutions exist in generic cases.

Let

p = π 4 , λ = 12 37 √ 2
. Then by [START_REF] Jackson | Classical Electrodynamics[END_REF] we have with values at nonzero integers depicted with asterisks. and condition ( 15) is fulfilled. If ω 2 = 37 25 , then k 0 = ±1 satisfies the resonance condition [START_REF] Lazarides | Discrete breathers in nonlinear magnetic metamaterials[END_REF]. Note that

M 2 = cos kπ 4 | k = {0, 1, . . . , 7} = 0, ± √ 2 2 , ±1 (56)
1 k 2 (1 -2λ cos kp) ≤ 1 k 2 (1 -2λ) = 37 k 2 (37 -12 √ 2) < 37 13k 2 < 37 25
whenever |k| ≥ 2. Thus k 0 = ±1 is the only resonant couple corresponding to ω 2 = 37 25 . This is the case of a simple resonance discussed in this section.

On the other side, if ω 2 = 37 1225 , then there are at least two couples k 0 = ±5, ±7 such that ( 16) is fulfilled. Since 37

k 2 (37 -12 √ 2) < 37 20k 2 < 37 1225 whenever |k| ≥ 8, it is sufficient to evaluate 1 k 2 (1-2λ cos kp) at k = {1, .
. . , 7} to determine the exact number of resonant modes. From Fig. 1, one can see that for ω 2 = 37 1225 there are two resonant couples k 0 = ±5, ±7. This is a double resonance which can be investigated analogically to the present section.

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

ω 2 U(z) -λU z - π 4 -λU z + π 4 + γ ωU (z) + U(z) + ω 2 U 2 (z) -λU 2 z - π 4 -λU 2 z + π 4 + γ ω(U 2 (z)) -ε cos z = 0 (57) for z ∈ R and parameters γ, λ, ω > 0, 1 2 = λ = √ 2
2 , and apply Theorem 1. Note that equation (57) is in the form of (2) with p = π 4 , h(z) = ε cos z, and M 2 of ( 13) is equivalent to (56), i.e., condition [START_REF] Kourakis | Self-focusing and envelope pulse generation in nonlinear magnetic metamaterials[END_REF] is satisfied. Hence using Remark 1(3) and Theorem 1, we obtain the existence of a unique solution U(ε, z) = O(ε) of equation (57) for any ε sufficiently close to 0. Since U(ε, z) is C ∞ -smooth, we can expand U(•, z) into Taylor series U(ε, z) = ε U 1 (z) + O(ε 2 ), and look for the function U 1 . It is easy to see that U 1 is uniquely determined by the equation

ω 2 U 1 (z) -λ U 1 z - π 4 -λ U 1 z + π 4 + γ ω U 1 (z) + U 1 (z) -cos z = 0.
We look for its solution in the form of U 1 (z) = A cos z + B sin z, and, by comparing coefficients at cos z and sin z, we obtain

-γ ωA + (1 + ω 2 ( √ 2λ -1))B = 0, (1 + ω 2 ( √ 2λ -1))A + γ ωB = 1.
Therefrom,

A = 1 + ω 2 ( √ 2λ -1) γ 2 ω 2 + (1 + ω 2 ( √ 2λ -1)) 2 , B = γ ω γ 2 ω 2 + (1 + ω 2 ( √ 2λ -1)) 2 .
In particular, if

λ = 12 37 √ 2 , ω = √ 37 5 , ( 58
) then A = 0, B = 5 √ 37 37γ , i.e., U(ε, z) = 5 √ 37 37γ ε sin z + O(ε 2 ). ( 59
)
Example 2. Now, moving onto resonant cases represented by equation ( 17), we consider the equation

ω 2 U(z) -λU z - π 4 -λU z + π 4 + εγ ωU (z) + U(z) + ω 2 U 2 (z) -λU 2 z - π 4 -λU 2 z + π 4 + εγ ω(U 2 (z)) -ε cos z = 0, (60)
i.e. when γ in (57) is replaced by εγ . Proceeding with point 2 of Remark 3 for parameters (58), then we obtain that k 0 = ±1, δ 2k 0 = -123 25 and υ = [START_REF] Longhi | Gap solitons in metamaterials[END_REF], then μ k 0 = 1 2 and ν k 0 = 0. So Theorem 2 can be applied to yield

2(1-δ 2k 0 ) δ 2k 0 = -296 123 . Since h(z) = cos z = 1 2 e ız + 1 2 e -ız in
(x(0), y(0)) = -3 δ 2k 0 2(1 -δ 2k 0 )(μ 2 k 0 + ν 2 k 0 ) (μ k 0 , ν k 0 ) = 3 123 74 1 2 , 0 . = (0.592281, 0). √ εc( 3 √ ε) e ız + 3 √ ε c( 3 √ ε) e -ız +U 2 3 √ εc( 3 √ ε) e ız + 3 √ ε c( 3 √ ε) e -ız , ε = 3 √ ε c(0) e ız + c(0) e -ız +O 3 √ ε = 3 √ ε 2x(0) cos z + O 3 √ ε = 3 123ε 74 cos z + O 3 ε 2 . = 3 √ ε1.18456 cos z + O 3 ε 2 (61)
for a smooth function c(•) with c(0) = x(0) + y(0)ı and ε = 0 small. Example 3. Similarly, Theorem 3 can be applied in the equation

ω 2 U(z) -λU z - π 4 -λU z + π 4 + εγ ωU (z) + U(z) + ω 2 U 2 (z) -λU 2 z - π 4 -λU 2 z + π 4 + εγ ω(U 2 (z)) -ε 2 cos z = 0, ( 62
) with λ = 12 37 √ 2 , ω = √ 37 5 , γ = 0.1. (63)
Note that cos z-term in (60) has been replaced by ε cos z in our case here, i.e., h(z) in ( 18) is replaced by εh(z). Then using the proof of Theorem 3, we obtain that (x(0), y(0)) = 0, -25 √ 37 k 0 , since k 0 = ±1. Using ( 22), ( 25) and (39), we derive that the solution from Theorem 3 of (62) with parameters (63) is given by εc(ε) e ık 0 z +ε c(ε) e -ık 0 z +U 2 (εc(ε) e ık 0 z +ε c(ε) e -ık 0 z , ε)

= ε(c(0) e ık 0 z + c(0) e -ık 0 z +O(ε)) = ε(-2y(0) sin k 0 z + O(ε)) = 50ε √ 37 k 0 sin k 0 z + O(ε 2 ) . = 8.21995ε sin z + O(ε 2 ). ( 64
)
Example 4. Next, we consider

ω 2 U(z) -λU z - π 4 -λU z + π 4 + εγ ωU (z) + U(z) + ω 2 U 2 (z) -λU 2 z - π 4 -λU 2 z + π 4 + εγ ω(U 2 (z)) -ε(μ 0 + 2(cos 2z -sin 2z)) = 0 (65)
with parameters (63), i.e., equation [START_REF] Litchinitser | Optical metamaterials: invisibility in visible and nonlinearities in reverse[END_REF] with the particular driving function h(z) = μ 0 + 2(cos 2zsin 2z). So, μ k 0 = ν k 0 = 0 and μ 2k 0 = ν 2 = -ν -2 = 1. Hence, Theorems 2 and 3 cannot be used. First, let k 0 = 1 and consider ε > 0, i.e. use "+"-sign in (48). Then we can derive that q 2 . = -0.294153 and y 1,2 = ± √ 0.214884 + 0.764897μ 0 . Note that 0 < 0.214884 + 0.764897μ 0 ≤ 0.429768. Using formula (54), we can estimate the Jacobian of H 3 (q 2 y 1,2 (q 2 ), y 1,2 (q 2 ), 0) as 0 < (0.218047

+ 0.776154μ 0 )v 2 ≤ 2.52553, (67)
since the formula (54) is equal to 1 v 2 det D H 3 (qy, y, 0). Therefore, we have 3 different small solutions for k 0 = 1, ε > 0 small. Using k 0 = -1, we obtain q 2 . = 0.294153, y 1,2 = ± √ 0.214884 + 0.764897μ 0 and by (54), estimate (67) follows. Hence for k 0 = ±1, ε > 0 small and μ 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)

√ εc( √ ε) e ık 0 z + √ ε c( √ ε) e -ık 0 z +U 2 ( √ εc( √ ε) e ık 0 z + √ ε c( √ ε) e -ık 0 z , ε) = √ ε(2x(0) cos k 0 z -2y(0) sin k 0 z) + O(ε) = 2 √ εy(0)(q 2 cos k 0 z -sin k 0 z) + O(ε) = ±2 ε(0.214884 + 0.764897μ 0 )(0.294153 cos z + sin z) + O(ε). ( 68
)
Moreover, there is another (smallest) solution of (65) with parameters (63) of order U(z) = O(ε). Using (45), we obtain x (0) = y (0) = 0 and by (46), the solution is

U(ε) = ε μ 0 - 50 123 cos 2z + 50 123 sin 2z + O(ε 2 ). ( 69
)
Considering k 0 = 1 and ε < 0, i.e. using "-"-sign in (48), we derive q 2 . = 1.83348 and y 1,2 = ± √ 0.0535297 -0.190543μ 0 . Note that 0 < 0.0535297 -0.190543μ 0 ≤ 0.107059. Formula (54) yields

0 < 0.218047 -0.776154μ 0 ≤ 0.436093, (70)
which implies that we also have 3 different small solutions for k 0 = 1, ε < 0 small and μ 0 satisfying (66). Note that if k 0 = -1, ε < 0, we obtain q 2 . = -1.83348, y 1,2 = ± √ 0.0535297 -0.190543μ 0 and, by (54), estimate (70). The 3 solutions of (65) with parameters (63) have the forms (69) and From the above computations, we can verify that formula (54) has one of the forms 0.218047 ± 0.776154μ 0 = 0. Since sgn υ = -1, we have that: if ε is small, ε sgn μ 0 > 0 and μ 0 fulfils (71), then there are 5 different small solutions (4 × O( √ |ε|) and 1 × O(ε)), if ε sgn μ 0 < 0 there is only 1 small solution of order O( √ |ε|), actually (69). For instance, if ε > 0 and μ 0 satisfies (71), we obtain for μ 0 > 0, and the only solution (69) of order O( √ ε) for μ 0 < 0.

√ -εc( √ -ε) e ık 0 z + √ -ε c( √ -ε) e -ık 0 z +U 2 √ -εc( √ -ε) e ık 0 z + √ -ε c( √ -ε) e -ık 0 z , ε = √ -ε(2x(0) cos k 0 z -2y(0) sin k 0 z) + O(ε) = 2 √ -εy(0)(q 2 cos k 0 z -sin k 0 z) + O(ε) = ±2 -ε(0.0535297 -0.190543μ 0 )(
k 0 q i y i 1 q 1 = 1.83348 y 1,2 (q 1 ) = ± √ -0.0535297 + 0.190543μ 0 q 2 = -0.294153 y 1,2 (q 2 ) = ± √ 0.214884 + 0.764897μ 0 -1 q 1 = -1.83348 y 1,2 (q 1 ) = ± √ -0.0535297 + 0.190543μ 0 q 2 = 0.294153 y 1,2 (q 2 ) = ± √ 0.214884 + 0.

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) = J j =1 B j cos (j -1) kz + C j sin j kz , ( 72
)
where k = 2π/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π 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

un = v n , ( 73
) [ vn (1 + 2U n )] = -γ v n (1 + 2U n ) -u n (1 + 2γ Un ) -2 [ Ün u n + 2 Un v n ], (74)
where

[• n ] = • n -λ • n-1 -λ• n+1 and U n (t) = U(ω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π/ω and with the site index n = 1, 2, . . . , N . The ith column of the monodromy matrix M is vector [u 1 (T ), . . . , u N (T ), v 1 (T ), . . . , v N (T )] t , that corresponds to the initial condition [u 1 (0), . . . , u N (0), v 1 (0), . . . , v N (0)] t equal to the ith 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 [START_REF] Agaoglou | Bifurcation results for travelling waves in nonlinear magnetic metamaterials[END_REF]. Note that the physically relevant interval for the coupling parameter is |λ| < 1/2 [START_REF] Diblík | Travelling waves in nonlinear magnetic metamaterials[END_REF].

It is because there will be a band of unstable multipliers when |λ| > 1/2, implying that all solutions of the system will be modulationally unstable. Following Section 4, we set λ = 12/ 37 √ 2 and p = π/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. ( 75
)
From ( 16) for γ = 0 the first resonance of the system is k = ±1, which corresponds to ω = √ 37/5. In Fig. 2, we depict the bifurcation diagrams of periodic solutions for a value of ω = 1.23 > √ 37/5. On the vertical axis is the oscillation amplitude A of the periodic solutions defined as

A = 1 2 (max (U (z)) -min (U (z)))
as a function of the driving amplitude h 0 . 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 γ = 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 → -h 0 , U → -U . When γ 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 turningpoint 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 ω 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 γ = 0, decreases towards zero. In the limit ω = √ 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 γ = 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 γ = 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 (μ 0 + 2(cos(2z)sin(2z))) , ( with μ 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 γ = 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 [START_REF] Agaoglou | Bifurcation results for travelling waves in nonlinear magnetic metamaterials[END_REF]. We observed that the typical dynamics of the instability are in the form of solution blow-up, similarly to that in [START_REF] Agaoglou | Bifurcation results for travelling waves in nonlinear magnetic metamaterials[END_REF].

As a conclusion, we found that the qualitative behaviour of the system (1) is similar to that of forced FPU lattices discussed in [START_REF] Fečkan | Periodic travelling waves of forced FPU lattices[END_REF], despite the significant difference in the coupling terms as explained in Section 1.

  k∈Z\{0,±k 0 } dk e kız j ∈Z\{0,±k 0 } dj e jız = l∈Z k+j =l dk dj e lız .

Remark 3 . 1 .

 31 Now we specify the O( √ |ε|) solution (solutions) of Theorem 4. Nontrivial roots of equation H 3 (x, y, 0) = 0 (47) can be explicitly calculated as follows. Denoting

  ) with the same fixed sign in all a, b, c, d depending on the sign of unscaled ε, this equation is equivalent to the system ax + by + x(x2 + y 2 ) = 0 cx + dy + y(x2 + y 2 ) = 0 (49) for a = a/υ, b = b/υ, c = c/υ, d = d/υ. Note υ = 0, i.e., δ 2k 0 = 1. It is easy to see that y = 0 whenever x = 0 and vice versa. Hence we assume x = 0 = y. Dividing the first equation by x, the second equation by y, and comparing the resulting equations, one derives

Fig. 1 .

 1 Fig. 1. The graph of functions 37 1225 and 1 x 2 (1-2λ cos xp)

  764897μ 0 Hence, equation (65) with parameters (63) has the solutions (69) and ±2 ε(-0.0535297 + 0.190543μ 0 )(1.83348 cos zsin z) + O(ε), ±2 ε(0.214884 + 0.764897μ 0 )(0.294153 cos z + sin z) + O(ε)

Fig. 2 .

 2 Fig. 2. (a) The oscillation amplitude of periodic solutions of (2) as a function of the driving amplitude h 0 for ω = 1.23 with γ = 0 (black thick lines) and γ = 0.001 (blue thin lines). (b) Solution profiles at the driving amplitudes indicated in the legend for vanishing γ . (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Fig. 3 .

 3 Fig. 3. The same as Fig. 2(a) for ω = √ 37/5 with (a) γ = 0 and (b) γ = 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.

Fig. 4 .

 4 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, γ = 0.

  However, using Theorem 4 one can obtain that (43) is of the form

				-8 +	15 129 625	37 2500	+ 4μ 2 0 = 0,
	i.e. |μ 0 | =	√	11 940 227 12 300	. = 0.280932. Moreover, (44) is equivalent to
					|μ 0 | <	√	11 940 227 12 300	.	(66)

Acknowledgments

The work of MA has been co-financed from resources of the operational program "Education and Lifelong Learning" of the European Social Fund and the National Strategic Reference Framework (NSRF) 2007-2013 within the framework of the Action State Scholarships 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. 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.