Discrete Breathers in Nonlinear Magnetic Metamaterials

Discrete breathers in nonlinear magnetic metamaterials. N. Lazarides 1,2 , M. Eleftheriou 1,3 and G. P. Tsironis 1 arXiv:cond-mat/0605674v2 [cond-mat.mtrl-sci] 11 Sep 2006 1 Department of Physics, University of Crete, and Institute of Electronic Structure and Laser, Foundation of Research and Technology, P. O. Box 2208, 71003 Heraklion, Greece 2 Department of Electrical Engineering, Technological Educational Institute of Crete, P. O. Box 140, Stavromenos, 71500, Heraklion, Crete, Greece 3 Department of Music Technology and Acoustics, Technological Educational Institute of Crete, E. Daskalaki, Perivolia, 74100 Rethymno, Crete, Greece Magnetic metamaterials composed of split-ring resonators or U −type elements may exhibit discreteness effects in THz and optical frequencies due to weak coupling. We consider a model onedimensional metamaterial formed by a discrete array of nonlinear split-ring resonators with each ring interacting with its nearest neighbours. On-site nonlinearity and weak coupling among the individual array elements result in the appearence of discrete breather excitations or intrinsic localized modes, both in the energy-conserved and the dissipative system. We analyze discrete single and multibreather excitations, as well as a special breather configuration forming a magnetization domain wall and investigate their mobility and the magnetic properties their presence induces in the system. PACS numbers: 63.20.Pw, 75.30.Kz, 78.20.Ci Artificial non-magnetic materials exhibiting magnetic properties in the Terahertz and optical frequencies have been recently predicted theoretically[1, 2] and demonstrated experimentally[3, 4, 5]. The key element for most of these magnetic metamaterials (MMs) is either the split-ring resonator (SRR) or its U −shaped modification[6]. The realization of MMs at such (and possibly higher) frequencies will affect substantially THz optics and their applications in devices of compact cavities, adaptive lenses, tunable mirrors, isolators and converters. Moreover, MMs with negative magnetic response can be combined with plasmonic wires that exhibit negative permittivity, producing left-handed materials (LHM), i.e. metamaterials with negative magnetic permeability µ and dielectric permitivity ǫ leading to a negative index of refraction [7, 8, 9, 10, 11]. In the present work we focus entirely on MMs that have the additional features of being nonlinear as well as discrete. While nonlinearity results in self-focusing, discreteness induces localization and, as a result, the combination of both leads in the generation of nonlinearly localized modes of the type of Discrete Breather (DB)[12, 13, 14, 15, 16]. These modes act like stable impurity modes that are dynamically generated and may alter propagation and emission properties of the system. We consider a planar one-dimensional (1D) array of N identical SRRs with their axes perpendicular to the plane; each unit is equivalent to an RLC oscillator with self-inductance L, Ohmic resistance R and capacitance C. The units become nonlinear due to the Kerr dielectric that fills their gap and
has permittivity equal to ǫ(|E|2 ) = ǫ0 ǫℓ + α|E|2 /Ec2 , where E is the electric component of the applied electromagnetic field, Ec is a characteristic electric field, ǫℓ the linear permittivity, ǫ0 the permittivity of the vacuum, and α = +1 (α = −1) corresponding to self-focusing (-defocusing) nonlinearity [17, 18, 19]. As a result, the SRRs aquire a field-dependent capacitance C(|E|2 ) = ǫ(|Eg |2 )A/dg , where A is the area of the cross-section of the SRR wire, Eg is the electric field induced along the SRR slit, and dg is the size of the slit. Since C(Un ) = dQn /dUn , the charge Qn stored in the capacitor of the n-th SRR is   Un2 Un , n = 1, 2, ..., N, (1) Qn = Cℓ 1 + α 3ǫℓ Uc2 where Cℓ = ǫ0 ǫℓ A/dg is the linear capacitance, Un = dg Egn is the voltage across the slit of the nth SRR, and Uc = dg Ec . Neighbouring SRRs are coupled due to magnetic dipole-dipole interaction through their mutual inductance M , which decays as the cube of the distance. For weak coupling between SRRs in a planar configuration, it is a good approximation to consider only nearest neighboring SRR interactions. Then, the dynamics of Qn and the current In circulating in the nth SRR is described by L dIn + RIn + f (Qn ) = M dt  dQn = In (2) dt  dIn−1 dIn+1 + dt dt + E, (3) where E is the electromotive force (emf) induced in each SRR due to the applied field, and f (Qn ) = Un is given implicitly from Eq. (1). The value of E at a given instant is proportional to the magnetic field component of the applied field perpendicular to the SRR plane, and/or the electric field component parallel to the side of the SRRs which contains the slit[20]. Using the relations ωℓ−2 = LCℓ , τ = tωℓ , Ic = Uc ωℓ Cℓ , Qc = Cℓ Uc , E = Uc ε, In = Ic in , Qn = Qc qn , Eqs. (2) and (3) can be normalized to dqn = in , dτ (4) 2 d (λ in−1 − in + λ in+1 ) = γ in − ε(τ ) + f (qn ), (5) dτ 1 1.5 where γ = RCℓ ωℓ is the loss coefficient, and λ = M/L is the coupling parameter. In the following, we use periodic boundary conditions (i.e., iN +1 = i1 , i0 = iN ) except otherwise stated. Analytical solution of Eq. (1) for un = f (qn ) with the conditions of un being real and un (qn = 0) = 0, gives the approximate expression α 3 q +3 f (qn ) ≃ qn − 3 ǫℓ n  α 3 ǫℓ 2 0.5 1 0.5 i n 0 0 −0.5 −0.5 −1 −1.5 0 qn5 + O(qn7 ) 0 (6) R qn That is, the on-site potential V (qn ) = 0 f (qn′ ) dqn′ is soft for focusing nonlinearity and hard for defocusing nonlinearity. Substituting qn = A cos (kDn − Ωτ ) into the linearized Eqs. (4) and (5) with ε = 0 and γ = 0, we obtain the frequency spectrum of linear excitations 7 25 50 14 n −1 t FIG. 1: (color online). Time evolution of a Hamiltonian breather for approximately two periods for λ = 0.02, Tb = 6.69, α = +1, ǫℓ = 2 and N = 50. 29 0.4 , (7) where Ω = ω/ωℓ is the normalized frequency, D is the separation of neighbouring SRR centers (unit cell size), and k the wavenumber (−π ≤ k D ≤ π). The parameter λ can be calculated numerically for any SRR geometry, since the magnetic field of the current circulating the SRR is well known. Here we estimate λ with a simple model[1], neglecting the effects of nonlinearity and coupling on the resonance frequency [18]. For not very small array dimensions, the inductance of a circular SRR of radius a with circular cross-section of diameter h, is L = µ0 a[ln(16a/h) − 1.5], where µ0 is the permeability of the vacuum. For a squared SRR with square cross-section with side length ℓ = 5 µm, t = w = dg = 1 µm the SRR depth, width, and slit size, respectively, length of unit cell D = 7 µm[4], and using that ℓ′ = 4(ℓ − w) − dg ispthe length of the axis of the wire, a = ℓ′ /2π, h = 4 w t/π, we arrive at L ≃ 6 × 10−12 H. For this L, we evaluate the capacitance necessary for providing √ the resonance frequency for a single SRR, fr ≃ 1/2π LCℓ = 6.2 T Hz, consistent with the available experimental information[4], to be Cℓ ≃ 11 × 10−17 F . Consider two neighbouring SRRs (1 and 2) in an array of circular SRRs of radius a with circular cross-section of diameter h. The flux Φ2 threading SRR 2 due to the induced magnetic field in SRR 1 B1 (r) ≃ µ0 SI1 /4πr3 +O((a/r)3 ), where I1 is the induced current in SRR 1, S = πa2 is the SRR area and r is the distance from its center (r ∼ D), is approximately Φ2 ≃ B1 (r = D)S. Then, M = Φ2 /I1 ≃ µ0 S 2 /4πD3 and λ ≃ (π/4)(a/D)3 /[ln(16a/h) − 1.5]. For an array of squared SRRs with square cross-section with√dimensions as in [4] we obtain λ ≃ (ℓ′ /D)3 /32π 2 [ln(4ℓ′ / πwt) − 1.5] ≃ 0.02. For silver made SRRs, whose conductivity and skin depth are σ ≃ 6.15 × 107 S/m and √ δ ∼ 20 nm, respectively, we obtain R = 2a/σhδ = ℓ′ /2σδ πwt ≃ 3.44, and γ ≃ 0.01. 28 0.2 XE Ωk = [1 − 2 λ cos(k D)] −1/2 27 0 26 25 0 -0.2 100 200 t 300 400 -0.4 500 0 10 20 n 30 40 50 FIG. 2: Moving Hamiltonian breathers. Right: breather amplitude for Tb = 6.69, α = +1, ǫℓ = 2, N = 50, and λ = 0.062 Left: space-time evolution of the center of energy XE for the breather shown in the right figure, for λp = 0.1 (solid line); 0.2 (dotted line); 0.3 (dashed line). We consider first the lossless case without applied field (γ = 0, ε = 0). Then, Eqs. (4) and (5) can be derived from the Hamiltonian  X 1 H= (8) q̇n2 + V (qn ) − λ q̇n q̇n+1 . 2 n For Hamiltonian systems DBs may be constructed from the anti-continuous (AC) limit[14], where all oscillators are uncoupled (λ → 0), obeying identical dynamical equations. Fixing the amplitude of one of them (say the one located at n = nb ) to a specific value qb , with the corresponding current ib = 0, we can determine the oscillation period Tb . An initial condition with qn = 0 for any n 6= nb , qnb = qb , and q̇n = in = 0 for any n, represents a trivial DB. Continuation of this solution for λ 6= 0 using the Newton’s method [14], results in numerically exact DBs up to λmax , where the linear excitation frequency band (which expands with increasing λ) reaches the DB frequency ωb = 2π/Tb . The linear stability of Hamiltonian DBs is addressed through the eigenvalues of the monodromy matrix (Floquet coefficients). Fig. 1 shows the time evolution of a typical, linearly stable, highly localized, DB excitation (λmax ∼ 0.067 for the chosen parameters). In this figure, plotted vs. time t and array site 3 n, is the normalized current in circulating the nth SRR. Another trivial DB can be obtained for qn = qb for any n 6= nb , qnb = 0 and q̇n = in = 0 for any n, corresponding to what we could call a ”dark” DB in analogy with the dark soliton in nonlinear continuous systems. Such a DB can be continued up to λ ∼ 0.025 but it is linearly unstable except for very small λ. In order to investigate the mobility of these DBs we followed the procedure described be Chen et al[16]. According to this work, in order to generate a (steady state) moving DB, having obtained a static DB (q0 , i0 = 0) by Newton’s method, we integrate Eqs. (4) and (5) using as initial condition (q(τ = 0), i(τ = 0)) = (q0 , i0 = 0) + λp (0, δi), where λp is the perturbation strength, and the perturbation vector δi corresponds to the current part of the (normalized to unity) pinning mode eigenvector. The resulting DB motion is followed by plotting the instantaneous center of localization of energy XE of the DB for several values of λp and a λ value close to λmax (Fig. 2); the parameter XE is defined as XE = N X n=1 1 (9) PN where En is the energy at site n and Etot = n=1 En . We note that Hamiltonian DBs move slowly through the lattice as a result of the perturbation. Their velocity decreases with increasing λ, although not uniformly with λp ; in particular, the DB velocity as a function of λ decreases faster as λp increases. In order to generate DBs for the forced and damped system we start by solving Eqs. (4) and (5) in the AC limit[14] with emf ε(τ ) = ε0 sin(Ωτ ). We identify two different amplitude attractors of the single SRR oscillator, with amplitudes qh ≃ 1.6086 and qℓ ≃ 0.2866 for the high and low amplitude attractor, respectively. Subsequently, we fix the amplitude of one of the oscillators (say the one at n = nb ) to qh and all the others to qℓ (in are all set to zero). Using this configuration as initial condition, we turn on adiabatically the coupling λ. The initial condition can be continued for λ 6= 0 leading to dissipative DB formation [14]. The time evolution of a typical dissipative DB is shown in Fig. 3. Both the DB and the background are oscillating with different amplitudes (high and low, respectively). This should be compared to the Hamiltonian DB in Fig. 1, where the background is always zero. By interchanging qh and qℓ in the initial conditions, we obtain another DB oscillating with low amplitude, while the background oscillates with high amplitude (Fig. 4). With appropriate initial conditions we can also obtain multi-breathers where two or more sites oscillate with high (low) amplitude, while the other ones with low (high) amplitude. Next, we fix the amplitude of half of the SRRs in the array (say those for n > N/2) to qh and the others to qℓ , and integrate Eqs. (4) and (5) from the AC limit with open-ended boundary 0.5 0.5 i n 0 0 −0.5 −0.5 −1 −1.5 0 0 −1 7 25 50 14 n t FIG. 3: (color online). Time evolution of a one-site dissipative breather during approximatelly two periods for Tb = 6.82, λ = 0.02, γ = 0.01, ε0 = 0.04, α = +1, ǫℓ = 2 and N = 50. 1.5 1 i n · En /Etot , 1 1.5 2 0.5 n0 0 −2 0 −0.5 0 25 −1 7 −1.5 n 50 14 t FIG. 4: (color online). Time evolution of a one-site dissipative breather of the second type (see text) during approximatelly two periods, for λ = 0.01 and the other parameters as in Fig. 3. conditions (i.e., iN +1 = i0 = 0). In this way we obtain an oscillating domain-wall, as shown in Fig. 5. The Hamiltonian DBs investigated are linearly stable, ensuring that they are not affected by small amplitude perturbations. On the other hand, the dissipative DBs are attractors for initial conditions in the corresponding basin of attraction and are robust against different kinds of small perturbations[15]. We analysed numerically and confirmed the stability of the DBs presented above under various kinds of perturbations; we followed the purturbed DB evolution for long time intervals (over 2 × 104 Tb ) without observing any significant change in the DB shape[21]. The dissipative system, which includes forcing due to the applied field, offers the possibility to study its magnetic response. Assume that the emf is induced by the magnetic component H = H0 cos(ωt), of the applied field. Then, at least for uniform solutions (In = I), the magnetization M = SI/D3 can be defined. In the direction perpendicular to the SRRs plane, the general 4 and much larger in magnitude than that. Thus, for large enough SRR arrays, one may obtain domain-wall DBs connecting domains of the array with positive and negative µ. 1.5 2 1 1 0.5 i n0 0 −0.5 −1 50 −2 0 25 7 t 14 0 −1 −1.5 n FIG. 5: (color online). Time evolution of a domain-wall breather during approximatelly two periods. Parameters as in Fig. 3. 2 (a) 1 0 −1 −2 0 7 t 14 6 4 2 0 −2 −4 −6 (b) 0 7 t 14 FIG. 6: (color online). Time evolution of κin (τ ) (red curve), compared with cos(Ωτ ) (black curve), and their sum (green curve), for two SRRs of the domain-wall DB excitation. (a) low amplitude current oscillation SRR (n = 15); (b) high amplitude current oscillation SRR (n = 35). Parameters as in Fig. 3. relation B = µ0 (H + M) gives B = B0 (cos(Ωτ ) + κi(τ )), (10) where B0 = ε0 Uc /SΩωℓ , and κ = µ0 S 2 Ω/ε0 D3 L. For the material parameters used above κ ≃ 3. From Eq. (10) negative magnetic response appears whenever the second term in the parentheses is larger in magnitude than the first one, and has opposite sign. Then, one may assign a negative µ to the medium. Without nonlinearity, the unique, uniform solution for the SRR array gives positive response below the resonance frequency (∼ ωℓ ). Nonlinearity allows the existence of multiple stable states, which makes it possible to obtain either positive or negative µ below ωℓ , depending on the state of the system. Moreover, exploiting DB excitations, MMs with domains of opposite sign magnetic responses can be created. In Fig. 6 we show the time evolution of cos(Ωτ ) and κin (τ ) as well as their sum, for two SRRs of the domain-wall DB (Fig. 5), relatively far from the domain-wall and the ends. The SRR with low amplitude current oscillation (n = 15) shows positive (paramagnetic) response. In contrast, the SRR with high amplitude current oscillation (n = 35) shows extreme diamagnetic (negative) response, since κin (τ ) is almost out of phase with cos(Ωτ ) In conclusion, a 1D planar array of nonlinear SRRs coupled through nearest-neighbor mutual inductancies was investigated numerically. The existence of DBs of various types, for both the energy-conserved and the dissipative system, was demonstrated. We found that longer range interaction does not affect the DB properties substantially[21]. We found that Hamiltonian DBs may be set into uniform motion under a small perturbation. We also obtained a special DB solution (magnetization domain-wall), which separates domains of the array with different magnetization. Multiple magnetization states are possible in this system due to nonlinearity, which allows either for negative or positive µ below resonance. Moreover, one can exploit multibreathers and domainwall DBs to create MMs with domains of opposite sign magnetic response. Discreteness effects may appear in SRR arrays with dimensions close to those reported in [4], even though the field wavelength is much larger than the array dimensions. We acknowledge support from the grant ”PYTHAGORAS II” (KA. 2102/TDY 25) of the Greek Ministry of Education and the European Union. [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] J. Zhou et al., Phys. Rev. Lett. 95, 223902 (2005). A. Ishikawa et al., Phys. Rev. Lett. 95, 237401 (2005). T. J. Yen et al., Science 303, 1494 (2004). N. Katsarakis et al., Opt. Lett. 30, 1348 (2005). A. N. Grigorenko et al., Nature 438, 335 (2005); C. Enkrich et al., Phys. Rev. Lett. 95, 203901 (2005). A. K. Sarychev and V. M. Shalaev, SPIE Proc. 5508, 128 (2004). V. P. Drachev et al., Laser Phys. Lett. 3, 49 (2006); V. M. Shalaev et al., Opt. Lett. 30, 3356 (2005). D. Smith et al., Phys. Rev. Lett. 84, 4184 (2000). R. Shelby et al., Science 292, 77 (2001). C. G. Parazzoli et al., Phys. Rev. Lett. 90, 107401 (2003). J. B. Pendry et al., Phys. Rev. Lett. 76, 4773 (1996); J. B. Pendry et al., IEEE Trans. Microwave Theory Tech. 47, 2075 (1999). A. J. Sievers and S. Takeno, Phys. Rev. Lett. 61, 970 (1988). R. S. MacKay and S. Aubry, Nonlinearity 7, 1623 (1994). J. L. Marin and S. Aubry, Nonlinearity 9, 1501 (1996). J. L. Marin et al., Phys. Rev. E 63, 066603 (2001). D. Chen et al., Phys. Rev. Lett. 77, 4776 (1996). A. A. Zharov et al., Phys. Rev. Lett. 91, 037401 (2003). S. O’Brien et al., Phys. Rev. B 69, 241101(R) (2004). N. Lazarides, and G. P. Tsironis, Phys. Rev. E 71, 036614 (2005). N. Katsarakis et al., Appl. Phys. Lett. 84, 2943 (2004). M. Eleftheriou, N. Lazarides and G. P. Tsironis, work in progress.