Chimera states in networks of locally and non-locally coupled SQUIDs
CChimera states in networks of locally and non-locally coupled SQUIDs
J. Hizanidis, N. Lazarides, G. P. Tsironis
Department of Physics, University of Crete, P. O. Box 2208, 71003 Heraklion, Greece (Dated: February 7, 2019)Planar and linear arrays of SQUIDs (superconducting quantum interference devices), operate asnonlinear magnetic metamaterials in microwaves. Such
SQUID metamaterials are paradigmaticsystems that serve as a test-bed for simulating several nonlinear dynamics phenomena. SQUIDs arehighly nonlinear oscillators which are coupled together through magnetic dipole-dipole forces due totheir mutual inductance; that coupling falls-off approximately as the inverse cube of their distance,i. e., it is non-local. However, it can be approximated by a local (nearest-neighbor) coupling which inmany cases suffices for capturing the essentials of the dynamics of SQUID metamaterials. For eithertype of coupling, it is numerically demonstrated that chimera states as well as other spatially non-uniform states can be generated in SQUID metamaterials under time-dependent applied magneticflux for appropriately chosen initial conditions. The mechanism for the emergence of these states isdiscussed in terms of the multistability property of the individual SQUIDs around their resonancefrequency and the attractor crowding effect in systems of coupled nonlinear oscillators. Interestingly,generation and control of chimera states in SQUID metamaterials can be achieved in the presenceof a constant (dc) flux gradient with the SQUID metamaterial initially at rest.
Keywords: SQUID metamaterials, magnetic metamaterials, chimera states, attractor crowding,synchronization-desynchronization transition
I. INTRODUCTION
The notion of metamaterials refers to artificially struc-tured media designed to achieve properties not availablein natural materials. Originally they were comprisingsubwavelength resonant elements, such as the celebratedsplit-ring resonator (SRR). The latter, in its simplest ver-sion, is just a highly conducting metallic ring with a slit,that can be regarded as an effectively resistive - inductive- capacitive (
RLC ) electrical circuit. There has been atremendous amount of activity in the field of metama-terials the last two decades, the results of which havebeen summarized in a number of review articles [1–8]and books [9–16]. One of metamaterial’s most remark-able properties is that of the negative refraction index ,which results from simultaneously negative dielectric per-mittivity and diamagnetic permeability.An important subclass of metamaterials is that of su-perconducting ones [17, 18], in which the elementaryunits (i. e., the SRRs) are made by a superconduct-ing material, typically Niobium (
N b ) [19] or NiobiumNitride (
N bN ) [20], as well as perovskite superconduc-tors such as yttrium barium copper oxide (
Y BCO ) [21].In superconductors, the dc resistance vanishes below acritical temperature T c ; thus, below T c , superconduct-ing metamaterials have the advantage of ultra-low losses,a highly desirable feature for prospective applications.Moreover, when they are in the superconducting state,these metamaterials exhibit extreme sensitivity in exter-nal stimuli such as the temperature and magnetic fields,which makes their thermal and magnetic tunability pos-sible [22]. Going a step beyond, the superconductingSRRs can be replaced by SQUIDs [23, 24], where theacronym stands for Superconducting QUantum Interfer-ence Devices . The simplest version of such a device con-sists of a superconducting ring interrupted by a Joseph- son junction (JJ) [25], as shown schematically in Fig.1(a); the most common type of a JJ is formed when-ever two superconductors are separated by a thin insu-lating layer (superconductor / insulator / superconduc-tor JJ). The current through the insulating layer and thevoltage across the JJ are then determined by the cele-brated Josephson relations. Through these relations, theJJ provides a strong and well-studied nonlinearity to theSQUID, which makes the latter a unique nonlinear oscil-lator that can be actually manipulated through multipleexternal means.
FIG. 1: (a) Schematic of a SQUID (superconducting quan-tum interference device) in a magnetic field. (b) Equivalentelectrical circuit. (c) Schematic top view of a one-dimensionalperiodic array of SQUIDs in a magnetic field H SQUID metamaterials are extended systems contain-ing a large number of SQUIDs arranged in various con-figurations which, from the dynamical systems point ofview, can be viewed theoretically as an assembly ofweakly coupled nonlinear oscillators that inherit the flex- a r X i v : . [ n li n . PS ] F e b ibility of their constituting elements (i.e, the SQUIDs).They present a nonlinear dynamics laboratory in whichnumerous classical as well as quantum complex spatio-temporal phenomena can be explored. Recent experi-ments on SQUID metamaterials have revealed several ex-traordinary properties such as negative permeability [26],broad-band tunability [26, 27], self-induced broad-bandtransparency [28], dynamic multistability and switch-ing [29], as well as coherent oscillations [30]. More-over, nonlinear effects such as localization of the discretebreather type [31] and nonlinear band-opening (nonlineartransmission) [32], as well as the emergence of counter-intuitive dynamic states referred to as chimera states incurrent literature [33–35], have been demonstrated nu-merically in SQUID metamaterial models [36].The chimera states, in particular, which were first dis-covered in rings of non-locally and symmetrically coupledidentical phase oscillators [37], have been reviewed thor-oughly in recent articles [38–40], are characterized by thecoexistence of synchronous and asynchronous clusters ofoscillators; their discovery was followed by intense the-oretical [41–61] and experimental [62–76] activities, inwhich chimera states have been observed experimentallyor demonstrated numerically in a huge variety of physicaland chemical systems.Here, the possibility for generating chimera states inSQUID metamaterials driven by a time-dependent mag-netic flux by proper initialization or by the applicationof a dc flux gradient is demonstrated numerically. TheSQUIDs in such a metamaterial are coupled togetherthrough magnetic dipole-dipole forces due to their mu-tual inductance. This kind of coupling between SQUIDsfalls-off approximately as the inverse cube of their center-to-center distance, and thus it is clearly non-local. How-ever, due to the magnetic nature of the coupling, itsstrength is weak [27, 30], and thus a nearest-neighborcoupling approach (i. e., a local coupling approach) isoften sufficient in capturing the essentials of the dynam-ics of SQUID metamaterials. Chimera states emerge inSQUID metamaterials with either non-local [33, 35] orlocal [34] coupling between SQUIDs. They can be gener-ated from a large variety of initial conditions, and theyare characterized using well-established measures. It isalso demonstrated that chimera states emerge in SQUIDmetamaterials with zero initial conditions using a dc fluxgradient; in that case, control over the obtained chimerastates can be achieved.In the next section, a model for a single SQUID thatrelies on the equivalent electrical circuit of Fig. 1(b) is de-scribed, and the dynamic equation for the flux throughthe ring of the SQUID is derived and normalized. InSection 3, the dynamic equations for a one-dimensional(1D) SQUID metamaterial with non-local coupling arederived, and subsequently they are reduced to the localcoupling limit. In Section 4, various types of chimerastates are presented and characterized using appropri-ate measures. In Section 5, the possibility to generateand control chimera states with a dc flux gradient, is ex- plored. A brief discussion and conclusions are presentedin Section 6. II. THE SQUID OSCILLATOR
The simplest version of a SQUID consists of a super-conducting ring interrupted by a JJ (Fig. 1(a)), whichcan be modeled by the equivalent electrical circuit of Fig.1(b); according to that model, the SQUID features a self-inductance L , a capacitance C , a resistance R , and a crit-ical current I c which characterizes an ideal JJ. A “real”JJ (brown-dashed square in Fig. 1(b)) is however mod-eled as a parallel combination of an ideal JJ, the resis-tance R , and the capacitance C . When a time-dependentmagnetic field is applied to the SQUID in a directiontransverse to its ring, the flux threading the SQUID ringinduces two types of currents; the supercurrent, whichis lossless, and the so-called quasiparticle current whichis subject to Ohmic losses. The latter roughly corre-sponds to the current through the branch containing theresistor R in Fig. 1(b). The (generally time-dependent)flux threading the ring of the SQUID is described in themodel as a flux source, Φ ext . Many variants of SQUIDshave been studied for several decades (since 1964) andthey have found numerous applications in magnetic fieldsensors, biomagnetism, non-destructive evaluation, andgradiometers, among others [77, 78]. SQUIDs exhibitvery rich dynamics including multistability, complex bi-furcation structure, and chaotic behavior [79].The magnetic flux Φ threading the ring of the SQUIDis given by Φ = Φ ext + L I, (1)where Φ ext is the external flux applied to the SQUID,and I = − C d Φ dt − R d Φ dt − I c sin (cid:18) π ΦΦ (cid:19) , (2)is the total current induced in the SQUID as provided bythe resistively and capacitively shunted junction (RCSJ)model of the JJ [80] (the part of the circuit in Fig. 1(b)contained in the brown-dashed square), Φ is the fluxquantum, and t is the temporal variable. The three termsin the right-hand-side of Eq. (1) correspond to the cur-rent through the capacitor C , the current through theresistor R , and the supercurrent through the ideal JJ,respectively. The combination of Eqs. (1) and (2) gives C d Φ dt + 1 R d Φ dt + I c sin (cid:18) π ΦΦ (cid:19) + Φ − Φ ext L = 0 . (3)Note that losses decrease with increasing Ohmic resis-tance R , which is a peculiarity of the SQUID device.The external flux usually consists of a constant (dc) termΦ dc and a sinusoidal (ac) term of amplitude Φ ac and fre-quency Ω, i. e., it is of the formΦ ext = Φ dc + Φ ac cos( ωt ) . (4)The normalized form of Eq. (3) be obtained by using therelations φ = ΦΦ , φ ac,dc = Φ ac,dc Φ , τ = tω − LC , Ω = ωω LC , (5)where ω LC = 1 / √ LC is the inductive-capacitive ( L C )SQUID frequency (geometrical frequency), and the defi-nitions β = I c L Φ = β L π , γ = 1 R (cid:114) LC . (6)for the rescaled SQUID parameter and the loss coeffi-cient, respectively. Thus we get¨ φ + γ ˙ φ + φ + β sin (2 πφ ) = φ dc + φ ac cos(Ω τ ) . (7)By substituting γ = 0 and φ ext = 0 and β sin (2 πφ ) (cid:39) FIG. 2: SQUID potential curves u SQ ( φ ) for β L = 0 . φ ac =0, and (a) φ dc = 0; (b) φ dc = 0 .
25; (c) φ dc = 0 .
5; (d) φ dc =0 .
75; (e) φ dc = 1 . β L φ into Eq. (7), we get ¨ φ + Ω SQ φ = 0, with Ω SQ = √ β L being the linear eigenfrequency (resonance fre-quency) of the SQUID. Eq. (7) can be also written as¨ φ + γ ˙ φ = − du SQ dφ , (8)where u SQ = − φ ext ( τ ) φ + 12 (cid:20) φ − βπ cos(2 πφ ) (cid:21) , (9)is the normalized SQUID potential, and φ ext ( τ ) = φ dc + φ ac cos(Ω τ ) , (10)is the normalized external flux. The SQUID potential u SQ given by Eq. (9) is time-dependent for φ ac (cid:54) = 0 andΩ (cid:54) = 0. Here, parameter values of β L less than unity( β L <
1) are considered, in accordance with recent ex-periments; in that case, u SQ is a single-well, althoughnonlinear potential. For φ ext = φ dc , there is no time-dependence; however, the shape of u SQ varies with vary-ing φ dc , as it can be seen in Fig. 2. The potential u SQ is symmetric around a particular φ for integer and half-integer values of φ dc . (In Figs. 2(a), (c), and (e), thepotential u SQ is symmetric around φ = 0, 0 .
5, and 1,respectively.) For all the other values of φ dc , the poten-tial u SQ is asymmetric; this asymmetry of u SQ allowsfor chaotic behavior to appear in an ac and dc driven single SQUID through period-doubling bifurcation cas-cades. Such cascades and the subsequent transition tochaos are prevented by a symmetric u SQ which rendersthe SQUID a symmetric system in which period-doublingbifurcations are suppressed [81]. Actually, suppression ofperiod-doubling bifurcation cascades due to symmetryoccurs in a large class of systems, including the sinu-soidally driven-damped pendulum. FIG. 3: Flux amplitude - driving frequency ( φ max − Ω)curves for a SQUID with β L = 0 . γ = 0 . φ dc = 0, and(a) φ ac = 10 − , (b) φ ac = 2 × − , (c) φ ac = 10 − , (d) φ ac = 10 − . For zero dc flux, the strength of the SQUID nonlinear-ity increases with increasing ac flux amplitude φ ac . Thiseffect is illustrated in Fig. 3 in which the flux amplitude- driving frequency ( φ max − Ω) curves, i. e., the reso-nance curves, for four values of φ ac spanning four ordersof magnitude are shown (for φ dc = 0). In Fig. 3(a), for φ ac = 0 . φ max − Ω curve is apparently symmetric around thelinear SQUID eigenfrequency, Ω SQ = √ β L (cid:39) . φ ac = 0 . φ ac = 0 .
01, the nonlinear effects arealready strong enough to generate a multistable φ max − Ωcurve. In Fig. 3(d), for φ ac = 0 .
1, the SQUID is in thestrongly nonlinear regime and the φ max − Ω curve hasacquired a snake-like form . Indeed, the curve “snakes”back and forth within a narrow frequency region via suc-cessive saddle-node bifurcations [79]. Note that in Figs.3(c) and (d), the frequency region with the highest mul-tistability is located around the geometrical frequency ofthe SQUID, i. e., at Ω (cid:39)
L C frequency in normal-ized units). Inasmuch the frequency at which φ max ishighest can be identified with the “resonance” frequencyof the SQUID, it can be observed that this resonance fre-quency lowers with increasing φ ac from the linear SQUIDeigenfrequency Ω SQ to the inductive-capacitive (geomet-rical) frequency Ω (cid:39)
1. Thus, the resonance frequencyof the SQUID, where its multistability is highest, can beactually tuned by nonlinearity, i. e., by varying the acflux amplitude φ ac . Note that the multistability of theSQUID is a purely dynamic effect, which is not relatedto any local minima of the SQUID potential (which is ac-tually single-welled for the values of β L considered here,i. e., for β L < FIG. 4: (a) Bifurcation diagram of φ ( nT ) as a function ofthe driving frequency Ω, for β L = 0 . γ = 0 . φ dc = 0 . φ ac = 0 .
18. (b) A typical chaotic attractor on the φ − ˙ φ phase-plane for Ω = 0 .
6. The other parameters are as in (a).
For φ dc (cid:54) = 0, chaotic behavior appears in wide fre-quency intervals below the geometrical frequency (Ω = 1)for relatively high φ ac . As it was mentioned above, theSQUID potential u SQ is asymmetric for φ dc (cid:54) = 0, andthus the SQUID can make transitions to chaos throughperiod-doubling cascades [79]. In the bifurcation diagramshown in Fig. 4(a), the flux φ is plotted at the end ofeach driving period T = 2 π/ Ω for several tenths of driv-ing periods (transients have been rejected) as a functionof the driving frequency Ω. This bifurcation diagram re-veals multistability as well as a reverse period-doublingcascade leading to chaos. That reverse cascade, specifi-cally, begins at Ω = 0 .
64 with a stable period-2 solution(i. e., whose period is two times that of the driving pe-riod T ). A period-doubling occurs at Ω = 0 .
638 resultingin a stable period-4 solution. The next period-doubling,at Ω = 0 .
62, results in a stable period-8 solution. Thelast period-doubling bifurcation which is visible in thisscale occurs at Ω = 0 .
614 and results in a stable period-16 solution. More and more period-doubling bifurca-tions very close to each other lead eventually to chaosat Ω = 0 . φ − ˙ φ phase plane in Fig. 4(b) for Ω = 0 . III. SQUID METAMATERIALS: MODELLINGA. Flux dynamics equations
Consider a one-dimensional periodic arrangement of N identical SQUIDs in a transverse magnetic field H as inFig. 1(c), which center-to-center distance is d and they are coupled through (non-local) magnetic dipole-dipoleforces [33]. The magnetic flux Φ n threading the ring ofthe n − th SQUID isΦ n = Φ ext + L I n + L (cid:88) m (cid:54) = n λ | m − n | I m , (11)where n, m = 1 , ..., N , Φ ext is the external flux in eachSQUID, λ | m − n | = M | m − n | /L is the dimensionless cou-pling coefficient between the SQUIDs at the sites m and n , with M | m − n | being their mutual inductance, and − I n = C d Φ n dt + 1 R d Φ n dt + I c sin (cid:18) π Φ n Φ (cid:19) (12)is the current in the n − th SQUID as given by the RCSJmodel [80]. The combination of Eqs. (11) and (12) gives C d Φ n dt + 1 R d Φ n dt + I c sin (cid:18) π Φ n Φ (cid:19) + 1 L N (cid:88) m =1 (cid:16) ˆΛ − (cid:17) nm (Φ m − Φ ext ) = 0 , (13)where ˆΛ − is the inverse of the symmetric N × N couplingmatrix with elements ˆΛ nm = (cid:26) , if m = n ; λ | m − n | = λ | m − n | − , if m (cid:54) = n, (14)with λ being the coupling coefficient betwen nearestneighboring SQUIDs. In normalized form Eq. (13) reads( n = 1 , ..., N )¨ φ n + γ ˙ φ n + β sin (2 πφ n ) = N (cid:88) m =1 (cid:16) ˆΛ − (cid:17) nm ( φ ext − φ m ) , (15)where Eq. (5) and the definitions Eq. (6) have beenused. When nearest-neighbor coupling is only taken intoaccount, Eq. (15) reduces to the simpler form¨ φ n + γ ˙ φ n + φ n + β sin (2 πφ n ) = λ ( φ n − + φ n +1 )+(1 − λ ) φ ext , (16)where λ = λ . B. Local and nonlocal linear frequency dispersion
Equation (11) with Φ ext = 0 can be written in matrixform as L ˆΛ (cid:126)I = (cid:126) Φ , (17)where the elements of the coupling matrix ˆΛ are givenin Eq. (14), and (cid:126)I , (cid:126) Φ are N − dimensional vectors withcomponents I n , Φ n , respectively. The linearized equationfor the current in the n − th SQUID, in the lossless case( R → ∞ ), is given from Eq. (12) as − (cid:126)I = C d dt (cid:126) Φ + 2 π I c Φ (cid:126) Φ , (18)where the approximation sin( x ) (cid:39) x has been employed.By substituting Eq. (18) into Eq. (17), we get ˆΛ (cid:18) ω LC d dt (cid:126) Φ + β L (cid:126) Φ (cid:19) + (cid:126) Φ = 0 . (19)In component form, the corresponding equation reads (cid:88) m ˆΛ nm (cid:18) ω LC d dt Φ m + β L Φ m (cid:19) + Φ n = 0 , (20)or, in normalized form (cid:88) m ˆΛ nm (cid:18) ω LC ¨ φ m + β L φ m (cid:19) + φ n = 0 , (21)where the overdots denote derivation with respect to thenormalized time τ .Substitute the trial (plane wave) solution φ n = exp i ( κn − Ω τ ) , (22)where κ is the dimensionless wavenumber (in units of d − ), into Eq. (21) to obtainΩ = 1 S (1 + β L S ) , (23)where S = (cid:88) m ˆΛ nm exp iκ ( m − n ) . (24)It can be shown that, for the infinite system, the funcion S is S = 1 + 2 λ ∞ (cid:88) s =1 cos( κs ) | s | = 1 + 2 λCi ( κ ) , (25)where s = m − n , and Ci ( κ ) is a Clausen function.Putting Eq. (25) into Eq. (23), we obtain the nonlocalfrequency dispersion for the 1D SQUID metamaterial asΩ κ = (cid:115) Ω SQ + 2 λβ L Ci ( κ )1 + 2 λCi ( κ ) , (26)where Ω SQ = 1 + β L . In the case of local (nearest-neighbor) coupling the Clausen function Ci ( κ ) is re-placed by cos( κ ). Then, by neglecting terms of order λ or higher, the local frequency dispersion Ω κ (cid:39) (cid:113) Ω SQ − λ cos( κ ) (27)is obtained. FIG. 5: Linear frequency dispersion Ω = Ω κ for nonlocal(red) and local (blue) coupling, for β L = 0 .
86, and (a) λ = − .
02, (b) λ = − .
04; (c) λ = − . The linear frequency dispersion Ω = Ω κ , calculated fornonlocal and local coupling from Eq. (26) and (27), re-spectively, is plotted in Fig. 5 for three values of the cou-pling coefficient λ . The differences between the nonlocaland local dispersion are rather small, especially for lowvalues of λ , i. e., for λ = − .
02 (Fig. 5(a)), which aremostly considered here. Although the linear frequencybands are narrow, the bandwidth ∆Ω = Ω max − Ω min in-creases with increasing λ . For simplicity, the bandwidth∆Ω can be estimated from Eq. (27); from that equa-tion the minimum and maximum frequencies of the bandcan be approximated by Ω min,max (cid:39) Ω SQ (cid:16) ± λ Ω SQ (cid:17) , sothat ∆Ω (cid:39) | λ | Ω SQ . (28)That is, the bandwidth is roughly proportional to themagnitude of λ . Note that for physically relevant pa-rameters, the minimum frequency of the linear band iswell above the geometrical (i. e., inductive-capacitive)frequency of the SQUIDs in the metamaterial. Thus, forstrong nonlinearity, for which the resonance frequency ofthe SQUIDs is close to the geometrical one (Ω = 1), noplane waves can be excited. It is this frequency regionwhere localized and other spatially inhomogeneous statessuch as chimera states are expected to emerge (given alsothe extreme multistability of individual SQUIDs there). IV. CHIMERAS AND OTHER SPATIALLYINHOMOGENEOUS STATES
Eqs. (15) are integrated numerically in time withfree-end boundary conditions ( φ N +1 = φ = 0) usinga fourth-order Runge-Kutta algorithm with time-step h = 0 .
02. The initial conditions have been chosen sothat they lead to chimera states. It should be noted thatchimera states can be obtained from a huge variety ofinitial conditions. Here we choose φ n ( τ = 0) = (cid:26) , for n (cid:96) < n ≤ n r ;0 , otherwise , (29)˙ φ n ( τ = 0) = 0 , (30) FIG. 6: Maps of (cid:104) ˙ φ n (cid:105) ( τ ) on the n − τ plane for β L = 0 . γ = 0 . λ = − .
02, Ω = 1 . N = 54, φ dc = 0, and (a) φ ac = 0 . φ ac = 0 .
04, (c) φ ac = 0 .
06, (d) φ ac = 0 .
08, (e) φ ac = 0 .
10, (f) φ ac = 0 . with n (cid:96) = 18 and n r = 36. Eqs. (15) are first inte-grated in time for a relatively long time-interval, 10 T time-units, where T = 2 π/ Ω is the driving period, sothat the system has reached a steady-state. While theSQUID metamaterial is in the steady-state, Eqs. (15)are integrated for τ sst = 1000 T more time-units. Then,the profiles of the time-derivatives of the fluxes, averagedover the driving period T , i. e., (cid:104) ˙ φ n (cid:105) T = 1 T (cid:90) T ˙ φ n dτ, n = 1 , ..., N, (31)are mapped as a function of τ . Such maps are shown inFig. 6, for several values of the ac flux amplitude, φ ac .In these maps, areas with uniform colorization indicatethat the SQUID oscillators there are synchronized, whileareas with non-uniform colorization indicate that theyare desynchronized.In Figs. 6(a) and (b), i. e., for low values of φ ac ,chimera states are not excited since the (cid:104) ˙ φ n (cid:105) T are practi-cally zero during the steady-state integration time. How-ever, this does not mean that the state of the SQUIDmetamaterial is spatially homogeneous, as we shall seebelow. For higher values of φ ac , chimera states beginto appear, in which one or more desynchronized clus-ters of SQUID oscillators roughly in the middle of theSQUID metamaterial are visible (Figs. 6(c)-(e)). Foreven higher values of φ ac , as can be seen in Fig. 6(f), thewhole SQUID metamaterial is desynchronized. In or-der to quantify the degree of synchronization for SQUIDmetamaterials at a particular time-instant τ , the magni-tude of the complex synchronization (Kuramoto) param-eter r is calculated, where r ( τ ) = | Ψ( τ ) | = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) n e πiφ n ( τ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (32)Below, two averages of r ( τ ) are used for the characteriza-tion of a particular state of SQUID metamaterials, i. e.,the average of r ( τ ) over the driving period T , (cid:104) r (cid:105) T ( τ ),and the average of r ( τ ) over the steady-state integration time (cid:104) r (cid:105) sst . These are defined, respectively, as (cid:104) r ( τ ) (cid:105) T = 1 T (cid:90) T r ( τ ) dτ, (cid:104) r (cid:105) sst = 1 τ sst (cid:90) τ sst r ( τ ) dτ. (33)The calculated (cid:104) r (cid:105) T ( τ ) for the states shown in Fig. 6,clarify further their nature. In Fig. 7(a), (cid:104) r (cid:105) T ( τ ) is plot-ted as a function of time τ for all the six states presentedin Fig. 6. It can be seen that for φ ac = 0 .
02 and 0 . (cid:104) r (cid:105) T ( τ ) is constant in time, although less than unity. Forsuch states, (cid:104) r (cid:105) T ( τ ) = (cid:104) r (cid:105) sst , where (cid:104) r (cid:105) sst can be in-ferred from Fig. 7(b) for the curves of interest to be (cid:104) r (cid:105) sst (cid:39) .
972 and (cid:104) r (cid:105) sst (cid:39) .
894 for φ ac = 0 .
02 and 0 . (cid:104) r (cid:105) T ( τ ) <
1. The exact natureof these partially synchronized states can be clarified byplotting the flux profiles φ n at the end of the steady-stateintegration time as shown in Figs. 7(c) and (d). In thesefigures, it can be observed that all but a few SQUID os-cillators are synchronized; in addition, those few SQUIDsexecute high-amplitude flux oscillations. Moreover, ithas been verified that the frequency of all the flux os-cillations is that of the driving, Ω. Such states can beclassified as discrete breathers/multi-breathers, i. e., spa-tially localized and time-periodic excitations which havebeen proved to emerge generically in nonlinear networksof weakly coupled oscillators [82]. In the present case,the multibreathers shown in Figs. 7(c) and (d) can befurther characterized as dissipative ones [83], since theyemerge through a delicate balance of input power andintrinsic losses. They have been investigated in some de-tail in SQUID metamaterials in one and two dimensions[31, 84–86], as well as in SQUID metamaterials on two-dimensional Lieb lattices [87].The corresponding (cid:104) r (cid:105) T ( τ ) for the states shown in Figs.6(c), (d), (e), and (f), are shown in 7(a) as green, blue,orange, and brown curves, respectively. In these curves FIG. 7: (a) The magnitude of the synchronization parameteraveraged over the driving period, (cid:104) r (cid:105) T , as a function of time τ for β L = 0 . γ = 0 . λ = − .
02, Ω = 1 . N = 54, φ dc = 0, and φ ac = 0 .
02 (black), φ ac = 0 .
04 (red), φ ac = 0 . φ ac = 0 .
08 (blue), φ ac = 0 .
10 (orange), φ ac = 0 . τ sst , (cid:104) r (cid:105) sst ,as a function of the ac flux amplitude φ ac . The other param-eters are as in (a). (c) The flux profile φ n for φ ac = 0 .
02 andthe other parameters as in (a). (d) The flux profile φ n for φ ac = 0 .
04 and the other parameters as in (a). there are apparently fluctuations around their temporalaverage over the steady-state integration time (shown inFig. 7(b)). These fluctuations are typically associatedwith the level of metastability of the chimera states [88,89]; an appropriate measure of metastability for SQUIDmetamaterials is the full-width half-maximum (FWHM)of the distribution of (cid:104) r (cid:105) T [33]. The FWHM can be usedto compare the metastability levels of different chimerastates. For synchronized (spatially homogeneous) andpartially synchronized states such as those in Figs. 6(a)and (b), the FWHM of the corresponding distribution ofthe values of (cid:104) r (cid:105) T is practically zero.Another set of initial conditions which gives rise tochimera states is of the form [34] φ n ( τ = 0) = 12 cos (cid:18) jπnN (cid:19) , ˙ φ n ( τ = 0) = 0 , (34)where n = 1 , ..., N . The initial conditions in Eq. (34) al-low for generating multiclustered chimera states, in whichthe number of clusters depends on j . In Figs. 8(a) and(b), maps of (cid:104) ˙ φ n (cid:105) T on the n − τ plane for j = 1 and j = 2, respectively, are shown. In Fig. 8(a), three largeclusters can be distinguished; in the two of them, theSQUID oscillators are synchronized, while in the thirdone, in between the two sychronized clusters, the SQUIDoscillators are desynchronized. The flux profile φ n ofthat state at the end of the steady-state integration time τ sst = 6000, is shown in Fig. 8(c) as blue circles (theblack curve is a guide to the eye) along with the ini-tial condition (red curve). It can be seen that two moredesynchronized clusters at the ends of the metamaterial,which are rather small (they consist of only a few SQUIDs each), are visible. Obviously, the synchronized clusterscorrespond to the spatial interval indicated by the almosthorizontal segments in the φ n profile. The corresponding (cid:104) ˙ φ n (cid:105) T map and flux profile φ n for j = 2 is shown in Figs.8(b) and (d), respectively. In this case, a number of six(6) synchronized clusters and seven (7) desynchronizedclusters are visible in both Figs. 8(b) and (d). In Figs.8(d), the red curve is the initial condition from Eq. (34)with j = 2. Chimera states with even more “heads” canbe generated from the initial condition Eq. (34) for j > N = 54).Similar chimera states can be generated with local(nearest-neighbor) coupling between the SQUIDs of themetamaterial. For that purpose, Eq. (16) is integratedin time using a fourth order Runge-Kutta algorithm withfree-end boundary conditions and the initial conditionsof Eq. (29). As above, in order to eliminate transientsand reach a steady-state, Eq. (16) is integrated for 10 T time units and the results are discarded. Then, (16) isintegrated for τ sst = 10 T more time units (steady-stateintegration time), and (cid:104) ˙ φ n (cid:105) T is mapped on the n − τ plane (Fig. 9). The emerged states are very similar tothose shown in Fig. 6, which is the case of non-localcoupling between the SQUIDs. In particular, the statesshown in Fig. 9(a), (b), and (c), have been generated forexactly the same parameters and initial-boundary con-ditions as those in Fig. 6(c), (e), and (f), respectively,i.e, for φ ac = 0 .
06, 0 .
1, and 0 .
12. Note that the stateof the SQUID metamaterial for φ ac = 0 .
12 is completelydesynchronized both in Fig. 6(f) and 9(c). One mayalso compare the plots of the corresponding (cid:104) r (cid:105) T as afunction of τ , which are shown in Fig. 9(d) for the localcoupling case. The averages of r over the steady-stateintegration time τ sst for φ ac = 0 .
06, 0 .
1, 0 .
12 are, re-spectively, (cid:104) r (cid:105) sst = 0 . . .
136 for the nonlocalcoupling case and (cid:104) r (cid:105) sst = 0 . . .
146 for thelocal coupling case. The probability distribution func-tion of the values of (cid:104) r (cid:105) T , pdf ( (cid:104) r (cid:105) T ), for the three statesin Figs. 9(a)-(c) are shown in Figs. 9(e)-(g), respec-tively. As it was mentioned above, the FWHM of sucha distribution is a measure of the metastability of thecorresponding chimera state. The FWHM for the distri-butions in Fig. 9(e) and (f), calculated for the chimerastates shown in Fig. 9(a) and (b), are respectively 0 . . values of (cid:104) r (cid:105) T have been used to obtain each of the threedistributions. Also, these distributions are normalizedsuch that their area sums to unity.The chimera states do not result from destabilizationof the synchronized state of the SQUID metamaterial; in-stead, they coexist with the latter, which can be reachedsimply by integrating the relevant flux dynamics equa-tions with zero initial conditions, i. e., with φ n ( τ = 0) = 0 FIG. 8: (a) Map of (cid:104) ˙ φ n (cid:105) T on the n − τ plane for β L = 0 . γ = 0 . λ = − .
02, Ω = 1 . N = 54, φ dc = 0, φ ac = 0 . j = 1. (b) Same as in (a) with initial conditions given by Eq. (34) with j = 2.(c) Flux profile φ n at the end of the steady-state integration time (blue circles, the black line is a guided to the eye), obtainedwith the initial conditions Eq. (34) with j = 1 (red curve). (d) Flux profile φ n at the end of the steady-state integration time(blue circles, the black line is a guided to the eye), obtained with the initial conditions Eq. (34) with j = 2 (red curve). and ˙ φ n ( τ = 0) = 0 for any n . In order to reach a chimerastate, on the other hand, appropriately chosen initialconditions such as those in Eq. (29) or Eq. (34) haveto be used. However, one cannot expect that the syn-chronized state is stable over the whole external param-eter space, i. e., the ac flux amplitude φ ac , the frequencyof the ac flux field Ω, and the dc flux bias φ dc . In or-der to explore the stability of the synchronized state ofthe SQUID metamaterial, the magnitude of the synchro-nization parameter averaged over the steady-state inte-gration time, (cid:104) r (cid:105) sst , is calculated and then mapped onthe φ dc − φ ac parameter plane. For each pair of φ ac and φ dc values, the SQUID metamaterial is initializedwith zeros (it is at “rest”). Once again, the frequencyΩ is chosen to be very close to the geometrical reso-nance Ω LC (Ω (cid:39) (cid:104) r (cid:105) sst on the φ dc − φ ac plane are shown for four driving frequencies Ωaround unity. These maps are a kind of “synchroniza-tion phase diagrams” , in which (cid:104) r (cid:105) sst = 1 indicates asynchronized state while (cid:104) r (cid:105) sst < . synchronization-desynchronization transitions do not go through a stagein which chimera states are generated; instead, the desta-bilization of a synchronized state results either in a com-pletely desynchronized state (light blue areas) or a clus-tered state (green areas). Thus, it seems that chimerastates cannot be generated when the SQUID metama-terial is initially at “rest”, i. e., with zero initial condi-tions. As we shall see in the next Section (Section 5),this is not true for a position-dependent external flux φ ext = φ ext ( n ). V. CHIMERA GENERATION BY DC FLUXGRADIENTSA. Modified flux dynamics equations
In obtaining the results of Fig. 10, a spatially homo-geneous dc flux φ dc over the whole SQUID metamaterialis considered. Although, all the chimera states presentedhere are generated at φ dc = 0, such states can be also gen-erated in the presence of a spatially constant, non-zero φ dc , by using appropriate initial conditions (not shownhere). In this Section, the generation of chimera states in FIG. 9: (a) Map of (cid:104) ˙ φ n (cid:105) T on the n − τ plane for β L = 0 . γ = 0 . λ = − .
02, Ω = 1 . N = 54, φ dc = 0, φ ac = 0 .
06, andinitial conditions given by Eq. (29). (b) Same as in (a) but with φ ac = 0 .
10. (c) Same as in (a) and (b) but with φ ac = 0 .
12. (d)The magnitude of the synchronization parameter averaged over the driving period, (cid:104) r (cid:105) T , as a function of time τ for φ ac = 0 . φ ac = 0 . φ ac = 0 .
12 (green). The other parameters are as in (a). (e) The distribution of 10 values of (cid:104) r (cid:105) T , pdf ( (cid:104) r (cid:105) T ), for the chimera state shown in (a). (f) Same as in (e) for the chimera state shown in (b). (g) Same as in (e) and (f)for the completely desynchronized state shown in (c). SQUID metamaterials driven by an ac flux and biased bya dc flux gradient is demonstrated, for the SQUID meta-material being initially at “rest”. The application of adc flux gradient along the SQUID metamaterial is exper-imentally feasible with the set-up of Ref. [28]. Considerthe SQUID metamaterial model in Section 3 . φ dc = φ dcn . Then,Eqs. (16) can be easily modified to become¨ φ n + γ ˙ φ n + φ n + β sin(2 πφ n ) = φ effn ( τ )+ λ ( φ n − + φ n +1 ) , (35)where φ effn = φ extn − λ ( φ extn − + φ extn +1 ) , (36)with φ extn = φ dcn + φ ac cos(Ω τ ) . (37)In the following, the dc flux function φ dcn is assumed tobe of the form φ dcn = n − N − φ dcmax , n = 1 , ..., N, (38)so that the dc flux bias increases linearly from zero (forthe SQUID at n = 1) to φ dcmax (for the SQUID at the n = N ). B. Generation and control of chimera states
Equations (35) are integrated numerically in time withfree-end boundary conditions (Eqs. (35)) using a fourth-order Runge-Kutta algorithm with time-step h = 0 . FIG. 10: Map of the magnitude of the synchronizationparameter averaged over the steady-state integration time, (cid:104) r (cid:105) sst , on the dc flux bias - ac flux amplitude ( φ dc − φ ac ) pa-rameter plane, for β L = 0 . γ = 0 . λ = − . N = 54,and (a) Ω = 1 .
03, (b) Ω = 1 .
02, (c) Ω = 1 .
01, (d) Ω = 0 . The SQUID metamaterial is initially at “rest”, i. e., φ n ( τ = 0) = 0 , ˙ φ n ( τ = 0) = 0 , n = 1 , ..., N. (39)This system is integrated for 10 T time units to elim-inate the transients and then for more τ sst = 10 T time units during which the temporal averages (cid:104) r (cid:105) sst and (cid:104) r (cid:105) T ( τ ) are calculated. Note that the transients die-outfaster in this case since the SQUID metamaterial is ini-tialized with zeros. Typical flux profiles φ n , plotted atthe end of the steady-state integration time are shownin Figs. 11(a)-(i). The varying parameter in this case is φ dcmax , which actually determines the gradient of the dcflux. The state of the SQUID metamaterial remains al-most homogeneous in space for φ dcmax increasing from zero0 FIG. 11: Flux profiles φ n as a function of n for β L = 0 . γ = 0 . λ = − . N = 54, φ ac = 0 .
04, Ω = 1 .
01, and(a) φ dcmax = 0 .
25; (b) 0 .
30; (c) 0 .
35; (d) 0 .
40; (e) 0 .
45; (f)0 .
50; (g) 0 .
55; (h) 0 .
60; (i) 0 .
65. (j) The magnitude of thesynchronization parameter averaged over the steady-state in-tegration time (cid:104) r (cid:105) sst as a function of φ dcmax for the parametersof (a)-(i) but with φ ac = 0 .
02 (black), 0 .
04 (red), 0 .
06 (green),0 .
08 (blue), 0 .
10 (magenta), 0 .
12 (brown). (k) Distributionsof the values of (cid:104) r (cid:105) T for φ ac = 0 .
04, and φ dcmax = 0 .
30 (black),0 .
40 (red), 0 .
50 (green), 0 .
60 (blue). The other parametersas in (a)-(i). The numbers next to the distributions are thecorresponding full-width half-maximums. to φ dcmax = 0 .
22. At that critical value of φ dcmax , the spa-tially homogeneous (almost synchronized) state breaksdown, for several SQUIDs close to n = N become desyn-chronized with the rest (because the dc flux is higher atthis end). The number of desynchronized SQUIDs for φ dcmax = 0 .
25 is about 6 − φ dcmax , more and more SQUIDs become desyn-chronized, until they form a well-defined desynchronizedcluster (Fig. 11(b) for φ dcmax = 0 . φ dcmax continuesto increase, the desynchronized cluster clearly shifts tothe left, i. e., towards n = 1 (Fig. 11(c)-(e)). Furtherincrease of φ dcmax generates a second desynchronized clus-ter around n = N for φ dcmax = 0 .
50 (Fig. 11(f)), whichpersists for values of φ dcmax at least up to 0 .
65. With theformation of the second desynchronized cluster, the firstone clearly becomes smaller and smaller with increasing φ dcmax (see Figs. 11(f)-(i)). Above, the expression “al-most homogeneous” was used instead of simply “homo-geneous”, because complete homogeneity is not possibledue to the dc flux gradient. However, for φ dcmax < . (cid:104) r (cid:105) sst are higher than 0 .
99 ( (cid:104) r (cid:105) sst > . (cid:104) r (cid:105) sst on φ dcmax for several values of the ac fluxamplitude φ ac is shown in Fig. 11(j). The SQUID meta-material remains in an almost synchronized state (with (cid:104) r (cid:105) sst > .
96 below a critical value of φ dcmax , which de-pends on the ac flux amplitude φ ac . That critical valueof φ dcmax is lower for higher φ ac . For values of φ dcmax higherthan the critical one, (cid:104) r (cid:105) sst gradually decreases until itsaturates at (cid:104) r (cid:105) sst (cid:39) .
12. For φ ac = 0 .
12, the SQUIDmetamaterial is in a completely desynchronized state for any value of φ dcmax (brown curve). The distributions ofthe values of (cid:104) r (cid:105) T , obtained during the steady-state in-tegration time, are shown in Fig. 11(k) for φ dcmax = 0 . .
40 (red), 0 .
50 (green), and 0 .
60 (blue). As ex-pected, the maximum of the distributions shifts to lower (cid:104) r (cid:105) T with increasing φ dcmax . These distributions have beendivided by their maximum value for easiness of presen-tation, and the number next to each distribution is itsfull-width half-maximum (FWHM). FIG. 12: The magnitude of the synchronization parameteraveraged over the steady-state integration time (cid:104) r (cid:105) sst mappedas a function of the ac flux amplitude and the maximum dcflux bias ( φ ac − φ dcmax plane), for β L = 0 . γ = 0 . N = 54,Ω = 1 .
01, and (a) λ = − .
02, (b) λ = − . Two typical “synchronization phase diagrams”, inwhich (cid:104) r (cid:105) sst is mapped on the φ ac − φ dcmax parameterplane, are shown in Figs. 12(a) and (b) for λ = − .
02 and λ = − .
06, respectively. The frequency of the driving acfield has been chosen once again to be very close to thegeometrical resonance of a single SQUID oscillator, i. e.,at Ω = 1 .
01. For each point on the φ ac − φ dcmax plane, Eqs.(35) are integrated in time with a standard fourth orderRunge-Kutta algorithm using the initial conditions of Eq.(39), with a time-step h = 0 .
02. First, Eqs. (35) areintegrated for 10 T time-units to eliminate transients,and then they are integrated for τ sst = 10 T more time-units during which (cid:104) r (cid:105) sst is calculated. A comparison be-tween Fig. 12(a) and (b) reveals that the increase of thecoupling strength between nearest-neighboring SQUIDsfrom λ = − .
02 to λ = − .
06 results in relatively mod-erate, quantitative differences only. In both Figs. 12(a)and (b), for values of φ ac and φ dcmax in the red areas,the state of the SQUID metamaterial is synchronized.For values of φ ac and φ dcmax in the dark-green, light-greenand light-blue areas, the state of the SQUID metama-terial is either completely desynchronized, or a chimerastate with one or more desynchronized clusters. In orderto obtain more information about these states, additionalmeasures should be used, such as the incoherence index S and the chimera index η [90, 91]. These are defined asfollows: First, define v n ( τ ) ≡ (cid:104) ˙ φ n (cid:105) T ( τ ) , (40)where the angular brackets denote averaging over T , and¯ v n ( τ ) ≡ n + 1 + n / (cid:88) n = − n / v n ( τ ) , (41)1the local spatial average of v n ( τ ) in a region of length n + 1 around the site n at time τ ( n < N is an integer).Then, the local standard deviation of v n ( τ ) is defined as σ n ( τ ) ≡ (cid:42)(cid:118)(cid:117)(cid:117)(cid:117)(cid:116) n + 1 + n / (cid:88) n = − n / ( v n − ¯ v n ) (cid:43) sst , (42)where the large angular brackets denote averaging overthe steady-state integration time. The index of incoher-ence is then defined as S = 1 − N N (cid:88) n =1 s n , (43)where s n = Θ( δ − σ n ) with Θ being the Theta function.The index S takes its values in [0 , chimera index is defined as η = N (cid:88) n =1 | s n − s n +1 | / , (44)and takes positive integer values. The chimera index η gives the number of desynchronized clusters of a (multi-)chimera state, except in the case of a completely desyn-chronized state where it gives zero. In Fig. 13, the in-coherence index S and the chimera index η are mappedon the φ ac − φ dcmax plane for the same parameters as inFig. 12(a). Figs. 13(a) and (b) provide more information FIG. 13: The index of incoherence S and the chimera index η are mapped on the φ ac − φ dcmax plane, for the same parametersas in Fig. 12(a) and n = 4, δ = 10 − . about the state of the SQUID metamaterial at a partic-ular point on the φ ac − φ dcmax plane. In Fig. 13(a), forvalues of φ ac and φ dcmax in the light-green area ( S = 0)the SQUID metamaterial is in a synchronized state (seethe corresponding area in Fig. 13(b) in which η = 0).For values of φ ac and φ dcmax in the red area ( S = 1),the SQUID metamaterial is completely desynchronized(the corresponding area in Fig. 13(b) has η = 0 due totechnical reasons). For values of φ ac and φ dcmax in one ofthe other areas, the SQUID metamaterial is in a chimerastate with one, two, or three desynchronized clusters, asit can be inferred from Fig. 13(b).Using the combined information from Figs. 12 and 13,the form of the steady-state of a SQUID metamaterial FIG. 14: Flux and voltage profiles φ n (blue) and v n = ˙ φ n (red), respectively, as a function of n for β L = 0 . γ = 0 . . φ ac = 0 .
04, and (a) φ dcmax = 0 .
2, (b) φ dcmax = 0 . φ dcmax = 0 . can be predicted for any physically relevant value of φ ac and φ dcmax . In Fig. 14, three flux profiles φ n are shownas a function of n , along with the corresponding profilesof their time-derivatives, ˙ φ n . The profiles in Figs. 14(a),(b), and (c), are obtained for φ ac = 0 .
04 and φ dcmax = 0 . .
4, and 0 .
6, respectively, which are located in the light-green, light-blue, and dark-green area of Fig. 14(b). Asit is expected, the state in Fig. 14(a) is an almost syn-chronized one, in Fig. 14(b) is a chimera state with onedesynchronized cluster, while in Fig. 14(c) is a chimerastate with two desynchronized clusters. At this point,the use of the expression “almost synchronized” shouldbe explained. In the presence of a dc flux gradient, itis impossible for a SQUID metamaterial to reach a com-pletely synchronized state. This is because each SQUIDis subject to a different dc flux, which modifies accord-ingly its resonance (eigen-)frequency. As a result, theflux oscillation amplitudes of the SQUIDs, whose oscilla-tions are driven by the ac flux field of amplitude φ ac andfrequency Ω, are slightly different. On the other hand,the maximum of the flux oscillations for all the SQUIDsis attained at the same time. Indeed, as can be observedin Fig. 14(a). the flux profile φ n is not horizontal, as itshould be in the case of complete synchronization. In-stead, that profile increases almost linearly from n = 1to n = N (that increase is related to the dc flux gradi-ent). However, the voltage profile ˙ φ n is zero for any n ,indicating that all the SQUID oscillators are in phase.Since, in such a state of the SQUID metamaterial thereis phase synchronization but no amplitude synchroniza-tion, the synchronization is not complete. However, thevalue of (cid:104) r (cid:105) sst in such a state is in the worst case higherthan 0 .
96 for moderately high values of φ ac = 0 . − . φ n profiles, also exhibit a very high degree of globalsynchronization ( (cid:104) r (cid:105) sst > . VI. DISCUSSION AND CONCLUSIONS
The emergence of chimera and multi-chimera states ina 1D SQUID metamaterial driven by an ac flux field isdemonstrated numerically, using a well-established modelthat relies on equivalent electrical circuits. Chimerastates may emerge both with local coupling betweenSQUID (nearest-neighbor coupling) and nonlocal cou-pling between SQUIDs which falls-off as the inverse cubeof their center-to-center distance. A large variety of ini-tial conditions can generate chimera states which persistfor very long times. In the previous Sections, the expres-sion “steady-state integration time” is used repeatedly;however, in some cases this may not be very accurate,since chimera states are generally metastable and sud-den changes may occur at any instant of time-integrationwhich results in sudden jumps the synchronization pa-rameter (cid:104) r (cid:105) T [33]. For the chimera states presented here,however, no such sudden changes have been observed.Along with the ac flux field, a dc flux bias, the same atany SQUID, can be also applied to the 1D SQUID meta-material. Chimera states can be generated in that caseas well, although not shown here.The emergence of those counter-intuitive states, theirform and their global degree of synchronization de-pends crucially on the initial conditions. If theSQUID metamaterial is initialized with zeros, the gen-eration of chimera states does not seem to be pos-sible for spatially constant dc flux bias φ dc . Inthat case, synchronization-desynchronization and reversesynchronization-desynchronization transitions may occurby varying the ac flux amplitude φ ac or the dc flux bias φ dc . In the former transition, a completely synchronizedstate suddenly becomes a completely desynchronized one.The replacement of the spatially constant dc flux bias bya position-dependent one, φ dcn , makes possible the genera-tion of chimera states from zero initial conditions. Here,a dc flux gradient is applied to the SQUID metamate-rial, which provides the possibility to control the chimerastate. Indeed, it is demonstrated that the position of thedesynchronized cluster(s) and the global degree of syn-chronization can be controlled to some extent by varyingthe dc flux gradient. Moreover, in the presence of a dcflux gradient, the ac flux amplitude controls the size ofthe desynchronized cluster.Here, the driving frequency is always chosen to bevery close to the geometrical frequency of the individualSQUIDs. In the case of relatively strong nonlinearity,considered here, the resonance frequency of individualSQUIDs is shifted to practically around the geometri-cal frequency. That is, for relatively strong nonlinearity,the driving frequency was chosen so that the SQUIDsare at resonance. For a single SQUID driven close toits resonance, the relatively strong nonlinearity makes ithighly multistable; then, several stable and unstable sin-gle SQUID states may coexist (see the snake-like curvespresented in Section 2). This dynamic multistability ef-fect is of major importance for the emergence of chimera states in SQUID metamaterials, as it is explained below.The dynamic complexity of N SQUIDs which are cou-pled together increases with increasing N ; this effect hasbeen described in the past for certain arrays of couplednonlinear oscillators as attractor crowding [92, 93]. Thiscomplexity is visible already for two coupled SQUIDs,where the number of stable states close to the geomet-rical resonance increases more than two times comparedto that of a single SQUID [34]; some of these states caneven be chaotic. Interestingly, the existence of homoclinicchaos in a pair of coupled SQUIDs has been proved byanalytical means [94, 95]. It has been argued that thenumber of stable limit cycles (i. e., periodic solutions) insuch systems scales with the number of oscillators N as( N − N .The multistability of individual SQUIDs around the res-onance frequency enhances the attractor crowding effectin SQUID metamaterial. Apart from the large numberof periodic solutions (limit cycles), a number of coex-isting chaotic solutions may also appear as in the two-SQUID system. All these states are available for eachSQUID to occupy. Then, with appropriate initializationof the SQUID metamaterial, or by applying a dc fluxgradient to it, a number of SQUIDs that belong to thesame cluster may occupy a chaotic state. The flux os-cillations of these SQUIDs then generally differ in boththeir amplitude and phase, resulting for that cluster tobe desynchronized. Alternatingly, a number of SQUIDsthat belong to the same cluster may find themselves ina region of phase-space with a high density of periodicsolutions. Then, the flux in these SQUID oscillators mayjump irregularly from one periodic state to another re-sulting in effectively random dynamics and in effect forthat cluster to be desynchronized. At the same time, theother cluster(s) of SQUIDs can remain synchronized and,as a result, a chimera state emerges. Conflict of Interest Statement
The authors declare that the research was conducted inthe absence of any commercial or financial relationshipsthat could be construed as a potential conflict of interest.
Author Contributions
NL conceived the structure of the manuscript, JH andNL performed the simulations, and all authors listed per-formed data analysis and have made intellectual contri-bution to the work, and approved it for publication.
Acknowledgments
This research has been financially supported by theGeneral Secretariat for Research and Technology (GSRT)3and the Hellenic Foundation for Research and Innovation (HFRI) (Code: 203). [1] Smith DR, Pendry JB, Wiltshire. Metamaterials andnegative refractive index.
Science
Vol. 305 no. 5685 (2004) 788–792.[2] Linden S, Enkrich C, Dolling G, Klein MW, Zhou J,Koschny T, et al. Photonic metamaterials: magnetism atoptical frequencies.
IEEE J. Selec. Top. Quant. Electron. (2006) 1097–1105.[3] Padilla WJ, Basov DN, Smith DR. Negative refractiveindex metamaterials. Materials Today (2006)28–35.[4] Shalaev VM. Optical negative-index metamaterials.
Na-ture Photon. (2007) 41–48.[5] Litchinitser NM, Shalaev VM. Photonic metamaterials. Laser Phys. Lett. (2008) 411–420.[6] Soukoulis CM, Wegener M. Past achievements and fu-ture challenges in the development of three-dimensionalphotonic metamaterials. Nature Photon. (2011) 523–530.[7] Liu Y, Zhang X. Metamaterials: a new frontier of scienceand technology. Chem. Soc. Rev. (2011) 2494–2507.[8] Simovski CR, Belov PA, Atrashchenko AV, Kivshar YS.Wire metamaterials: Physics and applications. Adv.Mater. (2012) 4229–4248.[9] Engheta N, Ziokowski RW. Metamaterials: Physics andengineering explorations (Wiley-IEEE Press) (2006).[10] Pendry JB.
Fundamentals and applications of negativerefraction in metamaterials (Princeton University Press)(2007).[11] Ramakrishna SA, Grzegorczyk T.
Physics and applica-tions of negative refractive index materials (SPIE andCRC Press) (2009).[12] Cui TJ, Smith DR, Liu R.
Metamaterials theory, designand applications (Springer: Springer) (2010).[13] Cai W, Shalaev V.
Optical metamaterials, fundamentalsand applications (Heidelberg: Springer) (2010).[14] Solymar L, Shamonina E.
Waves in metamaterials (NewYork: Oxford University Press) (2009).[15] Noginov MA, Podolskiy VA.
Tutorials in metamaterials (Taylor & Francis) (2012).[16] Tong XC.
Functional Metamaterials and Metadevices,Springer Series in Materials Science Vol.262 (SpringerInternational Publishing AG 2018) (2018).[17] Anlage SM. The physics and applications of supercon-ducting metamaterials.
J. Opt. (2011) 024001–10.[18] Jung P, Ustinov AV, Anlage SM. Progress in super-conducting metamaterials. Supercond. Sci. Technol. (2014) 073001 (13pp).[19] Jin B, Zhang C, Engelbrecht S, Pimenov A, Wu J, XuQ, et al. Low loss and magnetic field-tunable supercon-ducting terahertz metamaterials. Opt. Express (2010)17504–17509.[20] Zhang CH, Wu JB, Jin BB, Ji ZM, Kang L, Xu WW,et al. Low-loss terahertz metamaterial from supercon-ducting niobium nitride films. Opt. Express
20 (1) (2012) 42–47.[21] Gu J, Singh R, Tian Z, Cao W, Xing Q, He MX, et al.Terahertz superconductor metamaterial.
Appl. Phys.Lett. (2010) 071102 (3pp). [22] Zhang X, Gu J, Han J, Zhang W. Tailoring electromag-netic responses in terahertz superconducting metamate-rials. Frontiers of Optoelectronics (2014) 44–56.[23] Du C, Chen H, Li S. Quantum left-handed metamate-rial from superconducting quantum-interference devices.
Phys. Rev. B (2006) 113105–4.[24] Lazarides N, Tsironis GP. rf superconducting quantuminterference device metamaterials. Appl. Phys. Lett. (2007) 163501 (3pp).[25] Josephson B. Possible new effects in superconductivetunnelling.
Phys. Lett. A (1962) 251–255.[26] Butz S, Jung P, Filippenko LV, Koshelets VP, UstinovAV. A one-dimensional tunable magnetic metamaterial. Opt. Express
21 (19) (2013) 22540–22548.[27] Trepanier M, Zhang D, Mukhanov O, Anlage SM. Real-ization and modeling of rf superconducting quantum in-terference device metamaterials.
Phys. Rev. X (2013)041029.[28] Zhang D, Trepanier M, Mukhanov O, Anlage SM. Broad-band transparency of macroscopic quantum supercon-ducting metamaterials. Phys. Rev. X (2015) 041045[10 pages].[29] Jung P, Butz S, Marthaler M, Fistul MV, Lepp¨akangasJ, Koshelets VP, et al. Multistability and switching ina superconducting metamaterial. Nat. Comms. (2014)3730.[30] Trepanier M, Zhang D, Mukhanov O, Koshelets VP, JungP, Butz S, et al. Coherent oscillations of driven rf squidmetamaterials. Phys. Rev. E (2017) 050201(R).[31] Lazarides N, Tsironis GP, Eleftheriou M. Dissipative dis-crete breathers in rf squid metamaterials. Nonlinear Phe-nom. Complex Syst. (2008) 250–258.[32] Tsironis GP, Lazarides N, Margaris I. Wide-band tune-ability, nonlinear transmission, and dynamic multistabil-ity in squid metamaterials. Appl. Phys. A (2014)579–588.[33] Lazarides N, Neofotistos G, Tsironis GP. Chimeras insquid metamaterials.
Phys. Rev. B
91 (05) (2015)054303 [8 pages].[34] Hizanidis J, Lazarides N, Tsironis GP. Robust chimerastates in squid metamaterials with local interactions.
Phys. Rev. E (2016) 032219.[35] Hizanidis J, Lazarides N, Neofotistos G, Tsironis G.Chimera states and synchronization in magneticallydriven squid metamaterials. Eur. Phys. J.-Spec. Top. (2016) 1231–1243.[36] Lazarides N, Tsironis GP. Superconducting metamateri-als.
Phys. Rep. (2018) 1–67.[37] Kuramoto Y, Battogtokh D. Coexistence of coherenceand incoherence in nonlocally coupled phase oscillators.
Nonlinear Phenom. Complex Syst. (2002) 380–385.[38] Panaggio MJ, Abrams DM. Chimera states: Coexistenceof coherence and incoherence in network of coulped os-cillators.
Nonlinearity
28 (3) (2015) R67–R87.[39] Sch¨oll E. Synchronization patterns and chimera states incomplex networks: Interplay of topology and dynamics.
Eur. Phys. J.-Spec. Top. (2016) 891–919.[40] Yao N, Zheng Z. Chimera states in spatiotemporal sys- tems: Theory and applications. Int. J. Mod. Phys. B (2016) 1630002 (44 pages).[41] Abrams DM, Strogatz SH. Chimera states for coupledoscillators. Phys. Rev. Lett. (2004) 174102 [4pages].[42] Kuramoto Y, Shima SI, Battogtokh D, Shiogai Y. Mean-field theory revives in self-oscillatory fields with non-localcoupling.
Prog. Theor. Phys., Suppl. (2006) 127–143.[43] Omel’chenko OE, Maistrenko YL, Tass PA. Chimerastates: The natural link between coherence and incoher-ence.
Phys. Rev. Lett. (2008) 044105 [4 pages].[44] Abrams DM, Mirollo R, Strogatz SH, Wiley DA. Solvablemodel for chimera states of coupled oscillators.
Phys.Rev. Lett. (2008) 084103.[45] Pikovsky A, Rosenblum M. Partially integrable dynamicsof hierarchical populations of coupled oscillators.
Phys.Rev. Lett. (2008) 264103.[46] Ott E, Antonsen TM. Long time evolution of phase os-cillator systems.
Chaos (2009) 023117 [6 pages].[47] Martens EA, Laing CR, Strogatz SH. Solvable model ofspiral wave chimeras. Phys. Rev. Lett. (2010) 044101[4 pages].[48] Omelchenko I, Maistrenko Y, H¨ovel P, Sch¨oll E. Lossof coherence in dynamical networks: spatial chaos andchimera states.
Phys. Rev. Lett. (2011) 234102.[49] Yao N, Huang ZG, Lai YC, Zheng ZG. Robustness ofchimera states in complex dynamical systems.
Sci. Rep. (2013) 3522 [8 pages].[50] Omelchenko I, Omel’chenko OE, H¨ovel P, Sch¨oll E. Whennonlocal coupling between oscillators becomes stronger:Matched synchrony or multichimera states. Phys. Rev.Lett. (2013) 224101 [5 pages].[51] Hizanidis J, Kanas V, Bezerianos A, Bountis T. Chimerastates in networks of nonlocally coupled hindmarsh-roseneuron models.
Int. J. Bifurcation Chaos
24 (3) (2014)1450030 [9 pages].[52] Zakharova A, Kapeller M, Sch¨oll E. Chimera death:Symmetry breaking in dynamical networks.
Phys. Rev.Lett. (2014) 154101.[53] Bountis T, Kanas V, Hizanidis J, Bezerianos A. Chimerastates in a two-population network of coupled pendulum-like elements.
Eur. Phys. J.-Spec. Top. (2014) 721–728.[54] Yeldesbay A, Pikovsky A, Rosenblum M. Chimeralikestates in an ensemble of globally coupled oscillators.
Phys. Rev. Lett. (2014) 144103.[55] Haugland SW, Schmidt L, Krischer K. Self-organizedalternating chimera states in oscillatory media.
Sci. Rep. (2015) 9883.[56] Bera BK, Ghosh D, Lakshmanan M. Chimera states inbursting neurons. Phys. Rev. E (2016) 012205.[57] Shena J, Hizanidis J, Kovanis V, Tsironis GP. Turbulentchimeras in large semiconductor laser arrays. Sci. Rep. (2017) 42116.[58] Sawicki J, Omelchenko I, Zakharova A, Sch¨oll E.Chimera states in complex networks: interplay of frac-tal topology and delay. Eur. Phys. J.-Spec. Top. (2017) 1883–1892.[59] Ghosh S, Jalan S. Engineering chimera patterns in net-works using heterogeneous delays.
Chaos (2018)071103.[60] Shepelev IA, Vadivasova TE. Inducing and destruction ofchimeras and chimera-like states by an external harmonicforce. Phys. Lett. A (2018) 690–696. [61] Banerjee A, Sikder D. Transient chaos generates smallchimeras.
Phys. Rev. E (2018) 032220.[62] Tinsley MR, Nkomo S, Showalter K. Chimera and phase-cluster states in populations of coupled chemical oscilla-tors. Nature Phys. (2012) 662–665.[63] Hagerstrom AM, Murphy TE, Roy R, H¨ovel P,Omelchenko I, Sch¨oll E. Experimental observation ofchimeras coulped-map lattices. Nature Phys. (2012)658–661.[64] Wickramasinghe M, Kiss IZ. Spatially organized dynam-ical states in chemical oscillator networks: Synchroniza-tion, dynamical differentiation, and chimera patterns. Plos One (2013) e80586 [12 pages].[65] Nkomo S, Tinsley MR, Showalter K. Chimera statesin populations of nonlocally coupled chemical oscillators.
Phys. Rev. Lett. (2013) 244102.[66] Martens EA, Thutupalli S, Fourri´ere A, Hallatschek O.Chimera states in mechanical oscillator networks.
Proc.Natl. Acad. Sci.
110 (26) (2013) 10563–10567.[67] Sch¨onleber K, Zensen C, Heinrich A, Krischer K. Paternformation during the oscillatory photoelectrodissolutionof n-type silicon: Turbulence, clusters and chimeras.
NewJ. Phys. (2014) 063024 [10 pages].[68] Viktorov EA, Habruseva T, Hegarty SP, Huyet G, Kelle-her B. Coherence and incoherence in an optical comb. Phys. Rev. Lett. (2014) 224101 [5 pages].[69] Rosin DP, Rontani D, Haynes ND, Sch¨oll E, GauthierDJ. Transient scaling and resurgence of chimera statesin coupled boolean phase oscillators.
Phys. Rev. E (2014) 030902(R).[70] Schmidt L, Sch¨onleber K, Krischer K, Garc´ıa-Morales V.Coexistence of synchrony and incoherence in oscillatorymedia under nonlinear global coupling. Chaos (2014)013102.[71] Gambuzza LV, Buscarino A, Chessari S, Fortuna L,Meucci R, Frasca M. Experimental investigation ofchimera states with quiescent and synchronous domainsin coupled electronic oscillators. Phys. Rev. E (2014)032905.[72] Kapitaniak T, Kuzma P, Wojewoda J, Czolczynski K,Maistrenko Y. Imperfect chimera states for coupled pen-dula. Sci. Rep. (2014) 6379.[73] Larger L, Penkovsky B, Maistrenko Y. Laser chimeras asa paradigm for multistable patterns in complex systems. Nat. Comms. (2015) 7752.[74] Hart JD, Bansal K, Murphy TE, Roy R. Experimentalobservation of chimera and cluster states in a minimalglobally coupled network. Chaos (2016) 094801.[75] English LQ, Zampetaki A, Kevrekidis PG, SkowronskiK, Fritz CB, Abdoulkary S. Analysis and observationof moving domain fronts in a ring of coupled electronicself-oscillators. Chaos (2017) 103125.[76] Totz JF, Rode J, Tinsley MR, Showalter K, Engel H.Spiral wave chimera states in large populations of coupledchemical oscillators. Nature Phys. (2018) 282–285.[77] Clarke J, Braginski AI. The SQUID Handbook Vol. I:Fundamentals and Technology of SQUIDs and SQUIDSystems (Weinheim, Germany: Wiley-VCH) (2004).[78] Clarke J, Braginski AI.
The SQUID Handbook Vol.II: Applications of SQUIDs and SQUID Systems (Wein-heim, Germany: Wiley-VCH) (2004).[79] Hizanidis J, Lazarides N, Tsironis GP. Flux bias-controlled chaos and extreme multistability in squid os-cillators.
Chaos (2018) 063117 [8 pages]. [80] Likharev KK. Dynamics of Josephson Junctions and Cir-cuits. (Philadelphia: Gordon and Breach) (1986).[81] Swift JW, Wiesenfeld K. Suppression of period doublingin symmetric systems.
Phys. Rev. Lett. (1984) 705.[82] Flach S, Gorbach AV. Discrete breathers - advances intheory and applications. Phys. Rep. (2008) 1–116.[83] Flach S, Gorbach A. Discrete breathers with dissipa-tion.
Lect. Notes Phys. (Berlin Heidelberg: Springer-Verlag) (2008), pp. 289320.[84] Tsironis GP, Lazarides N, Eleftheriou M. Dissipativebreathers in rf squid metamaterials.
PIERS Online (2009) 26–30.[85] Lazarides N, Tsironis GP. Intrinsic localization in non-linear and superconducting metamaterials. Proc. SPIE (2012) 84231K.[86] Lazarides N, Tsironis GP. Nonlinear localization in meta-materials. Shadrivov I, Lapine M, Kivshar YS, editors,
Nonlinear, Tunable and Active Metamaterials (Switzer-land: Springer International Publishing) (2015), 281–301.[87] Lazarides N, Tsironis GP. Multistable dissipativebreathers and collective states in squid lieb metamate-rials.
Phys. Rev. E (2018) 012207. [88] Shanahan M. Metastable chimera states in community-structured oscillator networks. Chaos (2010) 013108.[89] Wildie M, Shanahan M. Metastability and chimera statesin modular delay and pulse-coupled oscillator networks. Chaos (2012) 043131.[90] Gopal R, Chandrasekar VK, Venkatesan A, LakshmananM. Observation and characterization of chimera states incoupled dynamical systems with nonlocal coupling. Phys.Rev. E (2014) 052914.[91] Gopal R, Chandrasekar VK, Venkatesan A, LakshmananM. Chimera at the phase-flip transition of an ensemble ofidentical nonlinear oscillators. Commun. Nonlinear Sci.Numer. Simulat. (2018) 30-46.[92] Wiesenfeld K, Hadley P. Attractor crowding in oscillatorarrays. Phys. Rev. Lett. (1989) 1335–1338.[93] Tsang KY, Wiesenfeld K. Attractor crowding in joseph-son junction arrays. Appl. Phys. Lett. (1990) 495–497.[94] Agaoglou M, Rothos VM, Susanto H. Homoclinic chaosin a pair of parametrically-driven coupled squids. Journalof Physics: Conference Series (2015) 012027.[95] Agaoglou M, Rothos VM, Susanto H. Homoclinic chaosin coupled squids.
Chaos, Solitons and Fractals99