Travelling waves in nonlinear magneto-inductive lattices

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


I. INTRODUCTION
In this work, we consider [27] where γ ≥ 0, λ, ω > 0, p = 0 are parameters and h ∈ C(R, R) is such that 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 [21] (see also [18] 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,29] and demonstrated experimentally in, e.g., [7,11,14,28].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,[22][23][24].Magneto-inductive metamaterials consisting of periodic arrays of SRRs have been made in one-, two-and three dimensions [20].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 [25] and higher-order nonlinear Schrödinger equations [26].In the second class, metamaterials are modelled by arrays of coupled oscillators, i.e. lattice equations, such as a nonlinear Klein-Gordon equation [19] and coupled Klein-Gordon equations [15,17].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 slit-ring resonators that compose the magneto-inductive materials [27].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 γ, λ and h.We also study the modulational stability of the periodic solutions by computing Floquet multipliers of the linearised 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 II 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 III, 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 IV.In Section V we solve the governing equations numerically.Comparisons with the analytical results are presented where we obtain good agreement.

II. EXISTENCE OF SMALL PERIODIC SOLUTIONS
Putting u n (t) = U (ωt + pn), z = ωt + pn in (1), we obtain We take Banach spaces with the norms respectively.It is easy to verify that Z Y X are compact embeddings and The following lemma is clear.

By setting
We have the next result.
Proof.Property ( 5) is obvious.Next, we use for c = −p, 0, p, and and so K ∈ L(Z, X) with is a constant depending on γ, λ, ω and p, then we also have Now we can prove the following existence result on (2) when all parameters except h are fixed.
Theorem 1. Assume (7) along with Then equation (2) has a unique solution where Moreover, U (h) can be approximated by an iteration process.Finally, it holds for any h 1 , h 2 ∈ X satisfying (9).
Proof.We rewrite (2) as a parameterized fixed point problem in Z.We already know that R : Z × X → Z is continuous and by ( 3), (8) such that Next, if there is ρ h > 0 such that then R(•, h) maps B(ρ 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 ρ H > 0, we need to solve A(r) = H together with DA(r) = 0.This implies and (9).So assuming (7), (9), we know that (12) has a positive solution ρ h < ρ H .We take the smallest one which clearly has the form (10).So R(•, h) maps B(ρ h ) into itself and, moreover, by ( 4), ( 8) The proof of the existence and uniqueness is finished by the Banach fixed point theorem [2].Next, let (9).By (4), ( 5) and ( 8), we derive which implies (11).
Remark 1.In the following special cases, (7) holds and we can replace Θ with the corresponding Θ i in the above considerations: 2. Let γ = 0.If p ∈ πQ, then we can write p = p1 p2 π for some p 1 ∈ Z, p 2 ∈ N, where p 1 , p 2 are relatively prime (their only common divisor is 1).Moreover, 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 To find a better relationship between M 1 and M 2 , we need to solve cos k1p1 p2 π = cos k2π p2 for k 1 ∈ Z\{0} and k 2 ∈ {0, 1, . . ., 2p 2 − 1}.This is equivalent to k1p1 p2 π = k2π p2 + 2πl for l ∈ Z and and for M 2 defined in (13), we get the existence of δ > 0 such that for each k ∈ Z\{0}.Thus, if p ∈ πQ and ( 14), ( 15) are valid, then Note that So if p ∈ πQ and (15) holds, there is at most a finite number of resonant modes k 0 ∈ Z\{0} determined by equation 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 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π p2 | 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 (7) is satisfied even without the non-resonance condition (14), since in each resonant mode k 0 we have and there is only the finite number of resonant modes.For each non-resonant mode k it holds Summarizing, if γ > 0, p ∈ πQ and (15) holds, condition (7) is fulfilled and where M is the set of resonant modes k 0 .

III. 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 with for U and ε 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 κ := |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 Now, we take the decomposition Applying the implicit function theorem to the first equation implies the existence of neighbourhoods Thus, equation ( 20) has the form Hence Expanding F and G of ( 18) into Taylor series using the last expansion of U 2 (U 1 , ε) yields ).Thus equation ( 21) is equivalent to or, using D ε U 2 (0, 0) and F of ( 18), Note that Next, setting d±k0 = d0 = 0, Therefore Analogically, denoting we have and the next identities can be proved Moreover, by definition of G in ( 18), Now, we put ( 26), ( 27), ( 28), ( 29) in (24) and collect coefficients by e k0ız and e −k0ız to obtain where Equation ( 30) is equivalent to Using resonance condition ( 16), we get Thus we have to solve or, separating the real and imaginary parts, Multiplying both equations (−1) and collecting terms by x and y we obtain 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 (19) and ( 20) are equivalent to ) by ( 23), (I − Q)BU 2 1 = 0 and then (34) has the form 0 ).To derive D 2 U1 u 2 (0), we twice differentiate (33) to obtain Next we derive Hence Thus, by (35) and the last identity, in (31) we have Summarizing, the bifurcation equation (32) has a form Now, we scale (36) and divide both equations by ε 3 .This leads to and the following result holds.
1. Now we specify the O( |ε|) solution (solutions) of Theorem 4. Nontrivial roots of equation can be explicitly calculated as follows.Denoting Therefrom x = q 1,2 y for (here q 1 corresponds to "+"-sign and q 2 to "−"-sign).As a consequence, putting x in the second equation of (49) one obtains (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. Accordingly, which is negative only for the minus sign (without any respect to the sign of unscaled ε).Thus the nontrivial solutions of (47) are On the other hand, assuming equation (47) still has nontrivial solutions.Indeed, in this case 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 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. 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 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 : for constants r 1 , r 2 independent of µ 0 .Inserting these expressions into (54), we have for polynomials 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.

IV. EXAMPLES
In this section we present several examples to illustrate the above results.
Example 1.First we study the non-resonant case from Section II.We consider the equation for z ∈ R and parameters γ, λ, ω > 0, 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 (15) 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.
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 where ) 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 (1).Note that the physically relevant interval for the coupling parameter is |λ| < 1/2 [6].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 IV, 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 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 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 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 ω is decreased toward the resonance frequency, the critical drive of h 0 for the occurrence of the first saddlenode 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) is 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 IV.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 (1).We observed that the typical dynamics of the instability is 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 I.
48)with the same fixed sign in all a, b, c, d depending on the sign of unscaled ε, this equation is equivalent to the systemax + 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., δ 2k0 = 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. 2 .
FIG. 2. (a)The oscillation amplitude of periodic solutions of (2) as a function of the driving amplitude h0 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 γ.

FIG. 5 .
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.