## Pattern formation and chimera states in 2D SQUID metamaterials

1 1 1

J. Hizanidis, N. Lazarides, and G. P. Tsironis Department of Physics, University of Crete, Herakleio, 71003, Greece National University of Science and Technology MISiS, Leninsky Prospect 4, Moscow, 119049, Russia a)

(Dated: 15 January 2020)

The Superconducting QUantum Interference Device (SQUID) is a highly nonlinear oscillator with rich dynamical behavior, including chaos. When driven by a time-periodic magnetic flux, the SQUID exhibits extreme multistability at frequencies around the geometric resonance which is manifested by a 'snake-like' form of the resonance curve. Repeating motifs of SQUIDs form metamaterials, i. e. artificially structured media of weakly coupled discrete elements that exhibit extraordinary properties, e. g. negative diamagnetic permeability. We report on the emergent collective dynamics in two-dimensional lattices of coupled SQUID oscillators, which involves a rich menagerie of spatio-temporal dynamics, including Turing-like patterns and chimera states. Using Fourier analysis we characterize these patterns and identify characteristic spatial and temporal periods. In the low coupling limit, the Turing-like patterns occur near the synchronization-desynchronization transition which can be related to the bifurcation scenarios of the single SQUID. Chimeras emerge due to the multistability near the geometric resonance, and by varying the dc component of the external force we can make them appear and reappear and, also, control their location. A detailed analysis of the parameter space reveals the coexistence of Turing-like patterns and chimera states in our model, as well as the ability to transform between these states by varying the system parameters.

The Superconducting QUantum Interference Device (SQUID), hereafter referred to as 'SQUID', is a highly nonlinear oscillator that exhibits strong resonant response to external magnetic fields. Its dynamics shows a wealth of phenomena such as hysteresis, multistability, subharmonic resonances, saddle-node and period-doubling bifurcations and chaos, which can be revealed through its complex bifurcation structure. Its resonance curve, in particular, acquires a 'snake-like' shape around its geometric resonance.

When many SQUIDs are arranged on a periodic array, they form magnetic metamaterials with extraordinary electromagnetic properties such as negative permeability, broad-band tunability, self-induced broad-band transparency, dynamic multistability and switching, as well as coherent oscillations. Besides their appeal as superconducting devices, SQUID metamaterials provide a unique testbed for exploring complex spatiotemporal dynamics. Here we demonstrate numerically that two-dimensional SQUID metamaterials (SQUID metasurfaces) support the emergence of certain spatially non-homogeneous dynamic states such as chimera states and patterned states of the Turing type.

Chimera states in SQUID metasurfaces make themselves apparent as domains of SQUIDs with synchronized (coherent) dynamics that coexist with domains of SQUIDs with desynchronized (incoherent) dynamics. Our system is an excellent physical and technologically relevant example of a driven system, where studies on chimeras are limited. Turing patterns, on the other hand, typically hexagons, stripes, rhombi, or labyrinths, emerge in reaction -diffusion systems. Since the SQUIDs in a metasurface are diffusively coupled, the emergence of Turing-like

a) Electronic mail: hizanidis@physics.uoc.gr

patterns are shown, for the first time, to be possible. Our investigation relies on a well-established model for SQUID metasurfaces, whose parameters lie in the experimentally accessible ranges of the applied constant (dc) flux as well as the amplitude and the frequency of the applied timedependent (ac) flux. We present numerically generated chimera states using appropriately selected initial conditions, as well as several Turing-like patterns whose characteristic length is determined from their corresponding two-dimensional spatial Fourier transforms. Moreover, in the low coupling limit, the region of stability for Turinglike patterns is related to the saddle-node bifurcation lines of the reduced equations for the SQUID metasurface.

The interplay between chimera states and Turing-like patterns for such a system of driven nonlinear oscillators is being addressed for the first time. Through a detailed analysis of the parameter space we reveal the coexistence and transformation between these states in certain regions.

## I. INTRODUCTION

Superconducting metamaterials comprising Superconducting QUantum Interference Devices (SQUIDs), are artificial materials that exhibit exceptional properties not found in nature, such as negative magnetic permeability, dynamic multistability 1,2 , broadband tunability, and self-induced broadband transparency 2 . Some of these extraordinary properties have been predicted theoretically both for the quantum 3 and the classical regime 4-6 . They can be implemented and studied in various designs and arrangements, both in one and two dimensions 1,2,7-12 . Recently, the degree of spatiotemporal coherence of SQUID metamaterials was examined experimentally and numerically using microwave transmission measurements 13 . Moreover, its quantum counterpart, the qubit, has been proposed as an important 'building block' of

quantum computers 14,15 .

Apart from their technological applications, SQUID metamaterials provide a unique testbed for exploring complex spatio-temporal dynamics. A SQUID metamaterial is in essence a system of nonlinear coupled oscillators with inertia, which are driven, damped, and are characterized by a nonlinear term which enters through the Josephson effect 16 . A very prominent dynamical feature of SQUID metamaterials are dissipative breathers 17 , which emerge as a result of their discreteness, the weak coupling between their elements and the Josephson nonlinearity.

In this work, we will study the collective behavior beyond these localized states, in a dynamical regime where the whole lattice is oscillating and is capable of creating Turing-like patterns 18 . Pattern forming systems have been in the center of scientific research for decades in a rich variety of natural and laboratory scenarios 19 . These include oscillatory chemical reactions governed by reaction-diffusion dynamics 20 , static and excitable biological media 21 , dissipative spatio-temporal solitons in nonlinear optics 22 , and many more. Here, we will use synchronization measurements in order to explore pattern formation in two-dimensional (2D) SQUID metamaterials.

Another phenomenon related to synchronization which has been observed in SQUID metamaterials are chimera states, where domains of coherent and incoherent motion coexist in an otherwise symmetric network of identical oscillators 23-27 . SQUID chimeras have mainly been studied in onedimensional (1D) arrays 28-31 . Here we will explore this phenomenon for locally coupled SQUIDs on a tetragonal lattice. Higher-dimensional chimeras have been the subject of recent works involving networks of Kuramoto and neuronal oscillators 32-35 . The interplay, however, between chimeras and Turing-like patterns has not been addressed sufficiently and this a new element that our work focuses on. Note that the system under study is a physical, technologically relevant example of a forced system, where such dynamics is still to be investigated 36,37 .

## II. SINGLE SQUID DYNAMICS

ASQUID consists of a superconducting ring interrupted by a Josephson junction (JJ) as shown schematically inside the dashed box of Fig. 1(a). When placed in a perpendicular, spatially uniform magnetic field H , a current I is induced which is the sum of the supercurrent Is flowing through the JJ and the quasiparticle current. Then, the magnetic flux F threading the loop of the SQUID is given by:

$$\Phi = \Phi _ { e x t } + L I, \quad \quad \ \ ( 1 ) \quad \text{hy}$$

where L is the self-inductance of the SQUID ring and F ext = F dc + F ac cos ( w t ) is the external flux applied to the SQUID, containing both a constant (dc) flux bias F dc and an alternating (ac) flux of amplitude F ac and frequency w .

The current I in the SQUID is given by the resistively and capacitively shunted junction (RCSJ) model of the JJ 38 ,

FIG. 1. (a) Schematic of a two-dimensional SQUID metamaterial in a magnetic field H t ( ) and (b) equivalent electrical circuit of the single SQUID (marked by the dashed box in (a)) in the RCSJ framework.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-2551a461-1e90-4fd5-8f3c-f6896c8af073.jpg)

schematically shown in Fig. 1(b):

$$I = - C \frac { d ^ { 2 } \Phi } { d t ^ { 2 } } - \frac { 1 } { R } \frac { d \Phi } { d t } - I _ { c } \sin \left ( 2 \pi \frac { \Phi } { \Phi _ { 0 } } \right ), \quad \ \ ( 2 )$$

where C is the capacitance of the JJ of the SQUID, R is the resistance, Ic is the critical current which characterizes the JJ, F 0 is the flux quantum, and t is the temporal variable. Combining Eqs. (1) and (2) we get:

$$\ddot { \phi } + \gamma \dot { \phi } + \phi + \beta \sin { ( 2 \pi \phi ) } = \phi _ { d c } + \phi _ { a c } \cos { ( \Omega \tau ) }, \quad ( 3 )$$

where all fluxes have been normalized to the flux quantum, f = F F / 0, f ac dc , = F ac dc , / F 0, while the frequency and time variable have been normalized to the inductive-capacitive SQUID frequency, w LC = 1 / √ LC and its inverse, respectively, i. e. W = w w / LC and t = t / w -1 LC . The parameter b = LIc / F 0 = b L / 2 p is the so-called rescaled SQUID parameter and g = w LCL / R corresponds to the loss coefficient.

Typical values of the design parameters of a SQUID 2,7 provide the dimensionless coefficients b glyph[similarequal] 0 1369 ( . b L glyph[similarequal] 0 86) . and g glyph[similarequal] 0 024 which appear in the normalized Eq. . (3) for the flux f = F F / 0 through the loop of the SQUID. They also provide experimentally plausible values f LC = w LC / ( 2 p ) glyph[similarequal] 13 9 . GHz ( W glyph[similarequal] 1) and f SQ = w SQ / ( 2 p ) glyph[similarequal] 18 9 . GHz ( W = W SQ glyph[similarequal] 1 364) for the geometric and the linear resonance fre-. quency of the SQUID, respectively 2,7,9 . The values of the externally controlled parameters f dc , f ac , and W used here, are within the range of the experimentally accessible values, i. e., f dc in the interval [ -1 2 , ] 7 , f ac in the interval [ 0 001 0 18 . , . ] 2 , and W in the interval 2 p w LC [ 10 22 5 , . ] GHz 7 .

By expanding the sine nonlinearity in Eq. 3 in a Taylor series and keeping the cubic term only, the SQUID model reduces to the famous driven Duffing oscillator. The latter is known to exhibit a nonlinear frequency response, bistability, hysteresis phenomena, and chaotic behavior. Similarly, the SQUID is capable of demonstrating complex dynamics, but with additional features owing to its higher-order nonlinear term. For a certain range of parameters the SQUID exhibits a 'snake-like' resonance curve in which multiple stable and unstable periodic orbits coexist and vanish through saddle-node bifurcations of limit cycles 29,39,40 . The detailed bifurcation structure for zero and finite dc flux was first reported in 29,39 .

FIG. 2. Resonance curves of a single SQUID for various dc flux values: (a) f dc = 0 0, (b) . f dc = 0 3, and (c) . f dc = 0 5. The vertical line marks . the value of the geometric resonance frequency and the insets show the corresponding stable periodic solutions at that value. Other parameters are: f ac = 0 06, . g = 0 024 and . b = 0 1369. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-059e0fde-8bf1-474b-a339-fb8935c925cc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-87059a00-a21a-4c61-b983-a7fe3e63d25f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-11da6f8e-d0b1-4131-93da-cb4cc65257ba.jpg)

Here we significantly further this analysis and explore the stability of solutions in the full ( f dc , W ) parameter space.

and so on.

Figure 2 shows the resonance curves of the single SQUID as the dc flux increases from 0 to 0.5. A saddle node bifurcation of limit cycles occurs at each turning point of the curve where stable and unstable branches merge 29,39 . The vertical line marks the geometric resonance frequency and the insets show the phase portraits of the corresponding stable periodic solutions at that particular value of W . As the dc flux increases, the 'center' of these solutions shifts to the right and the number of coexisting limit cycles at W LC changes. For f dc = 0 0 . (Fig. 2(a)) we have five coexisting periodic solutions of different amplitudes, centered around the origin. As the dc flux increases, the number of coexisting orbits gradually shrinks to one, at f dc = 0 30, while new subresonances 39,42 . make an appearance at lower frequencies (Fig. 2(b)). At f dc = 0 50 the . SQUID is again multistable with four coexisting periodic solutions (Fig. 2(c)) and, at the same time, the occurring saddlenode bifurcations have transformed from subcritical to supercritical and vice versa 39 . This scenario repeats itself periodically with the resonance curve moving back and forth, with respect to W LC , as we vary f dc . This 'rocking' of the snake-like resonance curve and, consequently, periodic appearance and disappearance of solutions around the geometric resonance frequency, is better visualized in the video SM1 of the Supplementary Material.

All of the aforementioned features are reflected in Fig. 3(a), where the co-dimension 2 bifurcation diagram in the ( f dc , W ) plane is depicted. The bifurcation lines have been obtained using a very powerful software tool that executes a root-finding algorithm for continuation of steady state solutions and bifurcation problems 43 . Cyan and red lines denote saddle-node bifurcations of limit cycles and period-doubling bifurcations, respectively. The bifurcation structure is extremely delicate and periodic in f dc with a period of unity. This periodicity can be proven as follows: Assuming that f is a solution of the single SQUID equation and by plugging f ± 1 into Eq. 3 we get: d 2 ( f ± 1 ) / dt 2 + g d ( f ± 1 ) / dt +( f ± 1 )+ b sin ( 2 pf ± 2 p ) = f dc + f ac cos ( W t ) . After simple manipulations we obtain: ¨ f + gf ˙ + f + b sin ( 2 pf ) = ( f dc ± 1 ) + f ac cos ( W t ) . Therefore, f satisfies Eq. 3 also for a dc flux f dc ± 1, and f dc ± 2,

Looking at Fig. 3(a) again, for fixed f dc values and moving in the W direction, we can recreate the resonance curves shown in Fig. 2 and the video SM1 of the Supplementary Material: The multiple and interwoven cyan lines correspond to the multiplicity of solutions around the geometric resonance frequency (better visible in the inset), while the red lines around f dc = 0 5 (and its symmetric . f dc = -0 5) are related . to the subresonances that make their appearance for those dc flux values. The bifurcation diagram of Fig. 3(a) presents additional, long period-doubling bifurcation branches extending to higher W values, which are not captured in the resonance curves of Fig. 2. The period doubling lines are symmetrical around f dc = ± 0 5 and for higher ac flux values are associated . with corresponding chaotic regions, as shown in 39 , where the maximum Lyapunov exponent was calculated in the ( f dc , W ) plane.

## III. TWO-DIMENSIONAL SQUID LATTICES

In this work, we will focus on two regimes of the driving frequency: The vicinity of W LC , and at lower values around W = 0 3. . Through the single SQUID complex dynamics, we aim at interpreting the collective behavior of the 2D SQUID lattice. We consider a planar N × N SQUID array consisting of identical units as shown in Fig. 1(a), arranged in an orthogonal lattice with a constant distance d in both x and y directions. The induced current Inm produces a magnetic field which couples each SQUID with all the others due to magnetic dipole-dipole interactions through their mutual inductance. To a good approximation, we may assume that the SQUIDS are coupled only to their nearest neighbors, neglecting further-neighbor interactions. The dynamic equations for the normalized flux through the ring of the ( n m , ) -th SQUID, f nm , are given by 17 :

$$& \ddot { \phi } _ { n m } + \gamma \dot { \phi } _ { n m } + \phi _ { n m } + \beta \sin { ( 2 \pi \phi _ { n m } ) } \\ = & \ \lambda ( \phi _ { n - 1, m } + \phi _ { n + 1, m } + \phi _ { n, m - 1 } + \phi _ { n, m - 1 } ) \\ + & \ ( 1 - 4 \lambda ) ( \phi _ { d c } + \phi _ { a c } \cos ( \Omega \tau ) ), \quad n, m = 1 \dots N, \quad ( 4 )$$

FIG. 3. Left: Bifurcation diagram in the ( f dc , W ) plane of the single SQUID oscillator. Blue and red lines correspond to saddle-node bifurcations of limit cycles and period doubling bifurcations, respectively. The inset shows a blowup around W = 1. Right: Value of the synchronization measure 〈 g 0 〉 in the ( f dc , W ) parameter space for a 30 × 30 SQUID lattice with coupling strength l = -0 01. The black and . cyan curves correspond to saddle-node bifurcation lines of the reduced system for l = 0 0 and . l = -0 01, respectively. Other parameters are: . f ac = 0 06, . g = 0 024 and . b = 0 1369. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-507aa765-2c25-4313-a107-aaf346a22f03.jpg)

where l ≡ M L / is the coupling constant between any two neighboring SQUIDs, coupled through their mutual inductance M . The value of M is negative due to the fact that the magnetic field generated by one SQUID crosses the neighboring SQUID in the opposite direction. In the following, we will study the nature of the synchronization-desynchronization transitions and will identify the collective states that emerge in relevant regimes of the parameter space. The latter involves the two parameters which can be easily tuned in an experiment, namely the dc flux and the frequency of the ac flux, with the other parameters (ac flux amplitude f ac , g and b ) kept constant.

Equations (4) are integrated numerically in time using a standard fourth-order Runge-Kutta algorithm with a time-step equal to 0 02 and periodic boundary conditions, i. e., . f n ( t ) = f N + ( n t ) for all n . This particular choice of boundary conditions does not affect the dynamics significantly. By using instead, e. g., free-end boundary conditions, only slight differences would have been observed which account for 1 -2% deviations of the corresponding g 0 value. Moreover, if nonlocal interaction between SQUIDs were assumed, a slight suppression of global synchronization would have been observed with both periodic and free-end boundary conditions.

The initial conditions for the f nm values follow a Gaussian random distribution in the interval [ -3 3 , ] and ˙ f nm = 0 0. We . will employ a quantification measure, originally introduced for the classification for chimera states 44 , the local curvature , which is calculated at each time instance by applying the absolute value of the discrete Laplacian on the spatial data of the magnetic flux:

$$\overset { \cdot } { e } \quad \hat { L } \phi _ { n m } ( t ) & = L _ { n m } ( t ) = 4 \phi _ { n, m } ( t ) - \phi _ { n + 1, m } ( t ) - \phi _ { n - 1, m } ( t ) \\ \overset { \cdot } { \cdot } \quad & - \phi _ { n, m + 1 } ( t ) - \phi _ { n, m - 1 } ( t ), \quad n, m = 1 \dots N. \quad ( 5 )$$

The local curvature is a measure for amplitude synchronization and in the synchronization regime it is close to zero while in the asynchronous regime it is finite and fluctuating. If g is the normalized probability density function of | ˆ L | , then g ( | ˆ L | = ) 0 measures the relative size of spatially coherent regions in each temporal realization and characterizes the entire lattice. For a fully synchronized system g ( | ˆ L | = ) = 0 1, while for a totally incoherent system it holds that g ( | ˆ L | = 0 ) = 0. An intermediate value of g ( | ˆ L | = ) 0 indicates the coexistence of synchronous and asynchronous SQUIDs and, therefore, potentially interesting collective behavior. From g , which is time-dependent, we calculate the spatial extent occupied by the coherent SQUIDs which is defined by the integral: g 0 ( ) = t ∫ d 0 g t ( , | ˆ L | ) d L | ˆ | , where d = 0 01 . L max is a threshold value distinguishing between coherence and incoherence and is related to the maximum local curvature ( L max ).

In order to correspond one single value to each realization, we calculate the time-average 〈 g 0 ( ) t 〉 , and we plot it in the ( f dc , W ) parameter space. The result is shown in Fig. 3(b) for a coupled lattice with l = -0 025. Yellow (bright) and purple . (dark) regions denote a synchronized and desynchronized lattice, respectively. The cyan lines mark two (for visualization simplicity) of the saddle-node bifurcation lines of the single SQUID system of Fig. 3(a). By comparing the two plots, it

is evident that the bifurcation lines of the single SQUID almost mark the borders between synchronization and desynchronization of the coupled system. For relatively weak coupling, which is the case in Fig. 3(b), this is plausible: When the single SQUID has one stable solution we may claim that the whole lattice acts like one SQUID and therefore the solution for the coupled system is the fully synchronized state. However, when the single SQUID loses its stability through the bifurcations shown in Fig. 3(a), each node of the lattice may behave differently resulting, thus, in a desynchronized state. For stronger coupling strengths (not shown here), the regions of incoherence around f dc = ± 0 5 broaden, while the . one around f dc = 0 shrinks. The structure of the parameter space (symmetry and periodicity), however, is maintained. In terms of mean field theory, the weak coupling case essentially corresponds to the single SQUID limit, where the collective and single SQUID behavior are seen to be in close correspondence. In the other extreme of strong coupling (not addressed in this work) one may possibly introduce an order parameter and proceed in analyzing the system from the point of view of collective phenomena.

## A. Pattern formation

Based on Fig. 3 we will study the collective dynamics emerging near the synchronization-desynchronization transition. Additionally, based on the resonance curve of the single SQUID (Fig. 2), we select the two W regimes that show the most interesting behavior: One around the geometric resonance frequency where the single SQUID is extremely multistable through successive saddle-node bifurcations of periodic solutions, and one at lower frequencies where period-doubling takes place. We prepare the lattice such that the initial conditions for the f nm values follow a Gaussian random distribution in the interval [ -3 3 , ] and ˙ f nm = 0 0. . As a control parameter we consider the coupling strength l , which, in principle, can be tuned in an experiment by increasing or decreasing the distance between the SQUIDs. The results are shown in Fig. 4: The top panels (a-d) show snapshots of the spatial distribution of the magnetic fluxes for the points (I-IV) marked in the parameter space in Fig. 3 (b). Figures 4 (a) and (b) correspond to a low driving frequency value ( W = 0 345) i. e. . far from the geometric resonance, where the single SQUID obtains low-amplitude periodic solutions and undergoes period doubling. For a coupling strength l = -0 025 the SQUID lat-. tice self-organizes into a labyrinthine-like pattern (Fig. 4(a)), while for a stronger coupling (Fig. 4 (b)) the collective state is a striped pattern, where smaller 'zigzag' patterns exist within each stripe. As mentioned previously, these patterns emerge from a completely random magnetic flux initialization and are, therefore, a result of the nonlinearity of the single SQUID and the collective dynamics of the coupled system.

In the middle panels, the two-dimensional Fourier power spectra | ˜ f k | 2 are plotted in the inverse space domain. The maximum values of the power spectra correspond to the characteristic wavenumber k = √ k x 2 + k y 2 of each pattern. From

| ˜ f k | 2 we obtain the 1D Radially-Averaged Power Spectrum (RAPS) 45 in terms of the wavelength l k = 2 p / k , shown in the lower panels of Fig. 4. From the peaks of these RAPSs we can extract the characteristic wavelength of each pattern. For example, for the pattern in Fig. 4(a) this value is glyph[similarequal] 3 53, . which is roughly the distance between two stripes in the m -direction of the corresponding plot in the top panel, in other words, the spatial period of the pattern. In the case of Fig. 4 (b), on the other hand, the RAPS obtains two maxima: The first one reflects the distance within the 'zigzag' patterns inside the stripes ( glyph[similarequal] 2) and the second one, the distance between the stripes themselves ( glyph[similarequal] 3 53). .

Similarly, Figs. 4(c) and (d) show the patterns, and the corresponding Fourier power spectra in k -space and RAPS, obtained near the geometric resonance, where the single SQUID may achieve high magnetic flux values through saddle-node bifurcations of limit cycles. The pattern in Fig. 4 (c), similar to Fig. 4 (a), is labyrinthine-striped, but with, evidently, a higher characteristic wavelength glyph[similarequal] 5 45. . For stronger coupling and a smaller f dc value, the emerging pattern consists of spots, with a characteristic wavelength equal to 4 (Fig. 4 (d)).

The patterns of Fig. 4 are spatio-temporal, and apart from a spatial period they also have a temporal period and corresponding frequency. These characteristic frequencies are given by the peaks of the Fourier power spectra in the inverse time domain, shown in the upper panels of Fig. 5. Figures 5(a) and (b) refer to Figs. 4 (a) and (c) (we have omitted the spectral analysis of cases Fig. 4 (b) and (d) because they are identical to Fig. 4 (a) and (c), respectively). We have plotted the spectra of all the SQUIDs in the lattice as well as their average (thick blue line). For the case of the pattern Fig. 4 (a), the spectra are very similar and the lattice is highly synchronized in frequency. As expected, the dominant frequency is that of the driving force, marked with the vertical dashed line. Moreover, the spectra are rather 'noisy' and they possess multiple secondary frequencies. This is typical for quasiperiodic motion as demonstrated by the phase diagrams in the lower panels of Fig. 5(a). The situation is similar, but 'cleaner' for the pattern in Fig. 4 (c). As seen in the power spectra of Fig. 5(b), the SQUIDs in the lattice are almost perfectly frequency-locked with the dominant frequency again being that of the driving force. The corresponding phase diagrams show, again, quasiperiodicity but the motion now is closer to harmonic since we are very close to the geometric resonance.

In correspondence with the traditionally discussed Turing patterns in Reaction-Diffusion systems, the 2D SQUID metamaterials can be characterized as partially cross-diffusive systems whose two components are the magnetic fluxes threading the loops of the SQUIDs and their time-derivatives. This can be readily inferred by taking the continuous limit of Eqs. 4. Although our system exhibits similarities with classical Reaction-Diffusion systems, it also exhibits differences, with the most important being the presence of the driving force. The emergence of Turing-like patterns in forced, discrete, 2D systems such as the one considered here has not been addressed very often in the literature, see e. g. 46 . Due to the forcing term in the dynamic equations, the simple pro-

FIG. 4. Top: Snapshots of the spatio-temporal patterns of the magnetic flux in a n × m SQUID lattice ( n = m = 30), corresponding to the cases I-IV marked on Fig. 3(b): (a) W = 0 345, . f dc = 0 5 and . l = -0 025 (case I), (b) . W = 0 345, . f dc = 0 5 and . l = -0 039 (case II), (c) . W = 1 06, . f dc = 0 3 and . l = -0 032 (case III), and (d) . W = 1 06, . f dc = 0 23 and . l = -0 05 (case IV). Other parameters are: . f ac = 0 06, . g = 0 024 and . b = 0 1369. . (See videos SM2(a)-(d) of the Supplementary Material for the corresponding videos). Middle: Corresponding Fourier Power Spectra in the 2D k -space. Bottom: Radially Averaged Power Spectrum (RAPS) in l k .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-e5a02a88-62d7-4f95-9a8f-2b71b64e36e0.jpg)

cedure to identify Turing instabilities cannot be applied. The main reason is the multistability of the individual SQUIDs, that results in a large number of periodic solutions for the SQUID metamaterial, even in the uncoupled case. Some of these solutions may be synchronized for l = 0 and may not be destroyed as the coupling is switched on. As a result, they may coexist with Turing-like patterns, indicating a complexity which cannot be handled by the simple analysis applied in classical Reaction-Diffusion systems.

## B. Chimera states

Chimeras are known to coexist with the fully synchronized state and, therefore, in many cases, they can be very sensitive to initial conditions. This holds for our system too, where chimeras can be achieved only for certain spatial distributions of the initial values of f . For example, in the locally coupled 1D SQUID array 29 , a 'sine wave' magnetic flux distribution was used for the initial conditions. It was shown that the SQUIDs that were prepared at lower values formed the coherent clusters of the chimera state, while those that were initially set at higher magnetic flux values, oscillated incoherently. Moreover, as the 'wavelength' of the initial magnetic flux distribution increased, so did the chimera state multiplicity (number of (in)coherent clusters). Note that in our system, since the frequency of the SQUID oscillators is imposed by the external driving, we are dealing with amplitude chimera states 41 . Here we will employ a set of different initial conditions, inspired by experimental feasibility. In particular, we will use a spatial gradient for the magnetic fluxes f nm = n -1 N -1 f max , where f max = 1 5 is the slope of the gradi-. ent, and zero values ( ˙ f nm = ) 0 for their derivatives. Another important factor for achieving robust chimeras in our system is the choice of the driving frequency. As reported in 29 , it is crucial to be near the geometric resonance where the phenomenon of 'attractor crwoding' 29 favors the emergence of such states. From section I, however, we know that by varying the dc flux, the snake-like form of the resonance curve shifts, resulting in loss of the SQUID multistability. It is interesting, therefore, to see what the effect of f dc will be on the creation of chimera states.

Figure 6 shows 3D snapshots of the magnetic flux (left), and their corresponding normalized local curvature values, when

FIG. 5. Top: Fourier power spectra in the frequency domain of all the SQUIDs and their average (thick blue line) for the spatio-temporal patterns of (a) Fig. 4 I, (b) Fig. 4 I, and (c) the chimera state of Fig. 6(a). Bottom: Phase diagrams of some typical timeseries in the SQUID lattice.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-2f2c12c5-f64d-404b-8545-f745bdd062d3.jpg)

the SQUID lattice is prepared with gradient initial conditions. In Fig. 6(a) the dc flux is zero and the incoherent cluster forms at the part of the lattice which is initially set at high magnetic flux values. On the other hand, the left half of the lattice is coherent and performs low-amplitude oscillations (better vizualized in the video SM2(e) of the Supplementary Material). By changing the f dc value to 0.3, the chimera state is destroyed and the collective state exhibits no spatio-temporal structure (Fig. 6(b)). This is due to the fact that for this dc flux value, the single SQUID is no longer multistable and chimera states are not possible. By further increasing f dc to unity, where the single SQUID is again multistable, the chimera reappears. Interestingly, comparing to Fig. 6(a), we observe a 'swap' in the position of the (in)coherent clusters, although the initial conditions are unchanged. This is due to the fact that for f dc = 1 the 'center' of the periodic solutions has shifted by 1 and the situation is reversed compared to Fig. 6(a) where f dc = 0 0. .

Apart from the single chimeras of Fig. 6, we can also achieve multichimera states (with more than one (in)coherent clusters), simply by increasing the slope of the initial conditions gradient. For instance, for a slope of 3 5, a multichimera . state with two (in)coherent domains is formed (not shown here). Recently, this mechanism for the generation of chimera states was reported, for non-identical coupled SQUIDs, where the gradient was in the dc flux distribution rather than in the initial conditions 47 . Such chimeras are similar to the equivalent one-dimensional structures 29 , extended in the second spatial dimension. For different special initial conditions, other types of chimeras are also possible, which are specific to the 2D geometry and are not present in the one-dimensional array. Here, however, we chose to focus on the 'stripe' chimeras of

Fig. 6, since the gradient flux initialization is easy to achieve experimentally.

Finally, we take a look at the Fourier power spectrum of the chimera state, namely that of Fig. 6(a), in the frequency domain. As we can see in Fig. 5(c), the sharpest peak is located at the value of the driving frequency, and there are secondary broader peaks at higher frequencies too. The corresponding phase diagram in the panel below shows some typical solutions of SQUIDs in the lattice. It is clear that the coexistence of smaller and bigger amplitude attractors (which is absent in the patterns discussed in Subsection III A) is the key to the emergence of chimera states in our system.

It should be noted that the emergence and form of chimera states does not depend crucially on the particular lattice geometry, as long as the dimensionality of the SQUID metamaterial is the same and the same initial conditions are used. For example, chimera states similar to those demonstrated in 48 for a 2D SQUID metamaterial on a Lieb lattice can be obtained for the 2D tetragonal lattice considered here, when the same initial conditions are used.

## C. Interplay of patterns

As demonstrated in the previous sections III A and III B, our system is capable of exhibiting Turing-like patterns associated with the single SQUID bifurcation structure in the low coupling limit, as well as chimera states when the driving frequency is chosen close to the geometric resonance frequency. Chimeras emerge through special initial conditions and may disappear and reappear as the dc flux varies; Turing-like pat-

FIG. 6. 3D snapshots of the magnetic flux (left) and the normalized discrete Laplacian (right) of a 30 × 30 SQUID lattice, for gradient initial conditions, l = -0 025, . W = 1 007 and different dc flux val-. ues: (a) f dc = 0 0 (see Supplementary Material for the corresponding . video SM2(e)), (b) f dc = 0 3, and (c) . f dc = 1 0. . Other parameters are: f ac = 0 06, . g = 0 024 and . b = 0 1369. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-7626dd9f-9880-4c68-8a83-56b15c22bd75.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-4685f801-2142-4550-aa9e-70065fd4bfc9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-21d40c35-e49d-471e-8e5e-eb5cbe389cad.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-2932b834-2e53-4795-9a85-6facc916253d.jpg)

terns, on the other hand, can be obtained for a random lattice initialization and a wider parameter range. Naturally, the question arises, under which circumstances do these different patterns coexist and how do they interact with each other as the system parameters change.

Figure 7 shows a map of the possible patterns observed in our system, in the ( f dc , l ) parameter space, for a driving frequency close to the resonance ( W = 1 03). . We distinguish areas of synchronized states (SYNCH), chimera states (CH), Turing-like states (TL), and states with no clear pattern structure which we will refer to as 'formless' states (FL). There are regions where different patterns may coexist, depending on the choice of initial conditions. We highlight two such examples, marked by points a and b. The corresponding coexisting states are shown in Fig. 8(a) and (b): In Fig. 8(a) a Turing-like patterns (left) coexists with a chimera state (right). Note that the unsychronized part of the chimera has an evident

FIG. 7. Classification of 2D patterns in the ( f dc , l ) parameter space for W = 1 03. . 'SYNCH' stands for synchronized state, 'TL' for Turing-like state, 'CH' for chimera state, and 'FL' for formless state. Points a and b correspond to ( f dc , l ) = ( 0 2 . , -0 03 . ) and ( 0 17 . , -0 025 , . ) respectively. Other parameters are: f ac = 0 06, . g = 0 024 and . b = 0 1369. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-1c878088-f864-47d2-8a98-56d92b486500.jpg)

spatial structure, resembling a 'half' Turing-like pattern. On the other hand, in Fig. 8(b), the chimera state (right) has a rather formless desynchronized domain and coexists with an unstructured pattern shown in the left.

FIG. 8. (a) Coexisting 'TL' state (left) and 'CH' state (right) for point a of Fig. 7. (b) Coexisting 'FL' state (left) and 'CH' state (right) for point b of Fig. 7. Plots (c)-(d) show the evolution of the 'CH' state of (b), as l increases along the arrow in Fig. 7. Specifically: (c) l = -0 024, (d) . l = -0 022, and (e) . l = -0 02. . Other parameters are: W = 1 03, . f ac = 0 06, . g = 0 024 and . b = 0 1369. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00004v2-figure-e9c99114-f528-4351-8e73-f355679a5799.jpg)

These coexisitng patterns are robust and do not transform between each other, as our long simulations (of the order of

10 4 periods) can confirm. The transformation between states, however, can be achieved by varying the system parameters, namely the coupling strength and the dc flux, as shown in Fig. 7. Specifically, a chimera state may evolve into a Turinglike pattern, but not the reverse since chimeras require special initial conditions in order to occur. Additionally, a formless state may also change into a Turing-like pattern and vice versa. By performing a continuation of states while varying the coupling strength l in the direction of the arrow in Fig. 7, we can see that the chimera state of Fig. 8(b) loses its structure and becomes a formless state in Fig. 8 (c). In turn, this state, by further increase of l , evolves into a Turing-like pattern as depicted in Figs. 8 (d) and (e). Interestingly, the coexistence of chimeras and Turing-like patterns has been reported before in two-dimensional networks of nonlocally coupled neurons in the low coupling limit 33 but, in general, is a question yet to be explored.

## IV. CONCLUSIONS

In conclusion, we have shown that a 2D SQUID lattice with nearest neighbor interactions is capable of exhibiting a rich menagerie of Turing-like pattern forming states. In the low coupling limit, this collective behavior emerges near the transition from synchronization to desynchronization where the single SQUID undergoes complex bifurcations. Moreover, near the geometric resonance, we observe 2D chimera states, as a result of the extreme multistability of the single SQUID. What is interesting is that by proper choice of initial conditions and tuning of the dc flux of the driving force, we are able to control the multiplicity and position of the chimera states, respectively. Additionally, in certain regions of the parameter space, chimeras may coexist with Turing-like patterns and also evolve into such, by proper tuning of the relevant parameters. Recent experiments 49 on the imaging of collective states in SQUID metamaterials through laser scanning microscopy (LSM technique), are very promising in terms of verifying our theoretical findings in the lab.

## V. SUPPLEMENTARY MATERIAL

See Supplementary Material for the videos corresponding to Fig. 2, the patterns of Fig. 4, and the chimera state of Fig. 6(a).

## ACKNOWLEDGEMENT

This work was financially supported by the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of NUST 'MISiS' (grant No. K3-2018-027). JH and NL acknowledge support by the General Secretariat for Research and Technology (GSRT) and the Hellenic Foundation for Research and Innovation (HFRI) (Code: 203). JH would also like to thank Jan Sieber for helping with the continuation tool.

- 1 P. Jung, S. Butz, M. Marthaler, M. V. Fistul, J. Leppäkangas, V. P. Koshelets, and A. V. Ustinov, 'Multistability and switching in a superconducting metamaterial', Nat. Commun. 5 , 3730 (2014).
- 2 D. Zhang, M. Trepanier, O. Mukhanov, and S. M. Anlage, 'Tunable Broadband Transparency of Macroscopic Quantum Superconducting Metamaterials', Phys. Rev. X 5 , 041045 (2015).
- 3 C. Du, H. Chen, and S. Li, 'Quantum left-handed metamaterial from superconducting quantum-interference devices', Phys. Rev. B 74 , 113105 (2006).
- 4 N. Lazarides and G. P. Tsironis, 'rf superconducting quantum interference device metamaterials', Appl. Phys. Lett. 90 , 163501 (2007).
- 5 N. Lazarides and G. P. Tsironis, 'Multistability and self-organization in disordered SQUID metamaterials', Supercond. Sci. Technol. 26 , 084006 (2013).
- 6 N Lazarides and G. P. Tsironis, 'Superconducting metamaterials', Phys. Rep. 752 , 1 (2018).
- 7 M. Trepanier, D. Zhang, O. Mukhanov, and S. M. Anlage, 'Realization and modeling of rf superconducting quantum interference devices', Phys. Rev. X 3 , 041029 (2013).
- 8 D. Zhang, M. Trepanier, T. Antonsen, E. Ott, and S. M. Anlage, 'Intermodulation in nonlinear SQUID metamaterials: Experiment and theory', Phys. Rev. B 94 , 174507 (2016).
- 9 S. Butz, P. Jung, L. V. Filippenko, V. P. Koshelets, and A. V. Ustinov, 'A one-dimensional tunable magnetic metamaterial', Opt. Express 21 , 22540 (2013).
- 10 P. Jung, A. V. Ustinov, and S. M. Anlage, 'Progress in superconducting metamaterials', Supercond. Sci. Technol. 27 , 073001 (2014).
- 11 A. V. Ustinov, 'Experiments With Tunable Superconducting Metamaterials', IEEE Trans. Terahertz Sci. Technol. 5 , 22 (2015).
- 12 E. I. Kiselev, A. S. Averkin, M. V. Fistul, V. P. Koshelets, and A. V. Ustinov, 'Two-tone spectroscopy of a SQUID metamaterial in the nonlinear regime', Phys. Rev. Research 1 , 033096 (2019).
- 13 M. Trepanier, D. Zhang, O. Mukhanov, V. P. Koshelets, P. Jung, S. Butz, E. Ott, T. M. Antonsen, A. V. Ustinov, and S. M. Anlage, 'Coherent oscillations of driven rf SQUID metamaterials', Phys. Rev. E 95 , 050201(R) (2017).
- 14 S. Saito, X. Zhu, R. Amsüss, Y. Matsuzaki, K. Kakuyanagi, T. Shimo-Oka, N. Mizuochi, K. Nemoto, W. J. Munro, and K. Semba, 'Towards Realizing a Quantum Memory for a Superconducting Qubit: Storage and Retrieval of Quantum States', Phys. Rev. Lett. 111 , 107008 (2013).
- 15 K. V. Shulga, E. Il'chev, M. V. Fistul, I. S. Besedin, S. Butz, O. V . Astafiev, U. Hübner, and A. V. Ustinov, 'Magnetically induced transparency of a quantum metamaterial composed of twin flux qubits', Nat. Comms. 9 , 150 (2018).
- 16 B. Josephson, 'Possible new effects in superconductive tunnelling', Phys. Lett. A 1 , 251 (1962).
- 17 N. Lazarides, G. P. Tsironis, and M. Eleftheriou, 'Dissipative discrete breathers in rf SQUID metamaterials', Nonlinear Phenom. Complex Syst. 11 , 250 (2008).
- 18 A. M. Turing, 'The Chemical Basis of Morphogenesis', Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, Vol. 237, No. 641. (Aug. 14, 1952), pp. 37-72.
- 19 M. C. Cross and P. C. Hohenberg, 'Pattern formation outside of equilibrium', Rev. Mod. Phys. 65 , 851 (1993).
- 20 K. Showalter and I. R. Epstein, 'From chemical systems to systems chemistry: Patterns in space and time', Chaos 25 , 097613 (2015).
- 21 A. J. Koch and H. Meinhardt, 'Biological Pattern Formation : from Basic Mechanisms to Complex Structures', Rev. Mod. Phys. 66 , 1481 (1994).
- 22 F. T. Arecchi, S. Boccaletti, and P. Ramazza, 'Pattern formation and competition in nonlinear optics', Phys. Rep. 318 , 1-83 (1999).
- 23 Y. Kuramoto and D. Battogtokh, 'Coexistence of coherence and incoherence in nonlocally coupled phase oscillators', Nonlinear Phenom. Complex Syst. 5 , 380 (2002).
- 24 D. M. Abrams and S. H. Strogatz, 'Chimera states for coupled oscillators', Phys. Rev. Lett. 93 , 174102 (2004).
- 25 M. J. Pannagio and D. Abrams, 'Chimera states: Coexistence of coherence and incoherence in network of coulped oscillators', Nonlinearity 28 , R67 (2015).
- 26 N. Yao and Z. Zheng, 'Chimera states in spatiotemporal systems: Theory and applications', Int. J. Mod. Phys. B 30 1630002 (2016).

- 27 O. E. Omel'chenko, 'The mathematics behind chimera states', Nonlinearity 31 , R121 (2018).
- 28 N. Lazarides, G. Neofotistos, and G. P. Tsironis, 'Chimeras in SQUID Metamaterials', Phys. Rev. B 91 , 054303 (2015).
- 29 J. Hizanidis, N. Lazarides, and G. P. Tsironis, 'Robust chimera states in SQUID metamaterials with local interactions', Phys. Rev. E 94 , 032219 (2016).
- 30 J. Hizanidis, N. Lazarides, G. Neofotistos, and G. P. Tsironis, 'Chimera states and synchronization in magnetically driven SQUID metamaterials', Eur. Phys. J. Special Topics, Springer 225 , 1231 (2016).
- 31 A. Banerjee and D. Sikder, 'Transient chaos generates small chimeras', Phys. Rev. E 98 , 032220 (2018).
- 32 Y. Maistrenko, O. Sudakov, O. Osiv, and V. Maistrenko, 'Chimera states in three dimensions', New J. Phys. 17 , 073037 (2015).
- 33 A. Schmidt, T. Kasimatis, J. Hizanidis, A. Provata, P. Hövel, 'Chimera patterns in two-dimensional networks of coupled neurons', Phys. Rev. E 95 , 032224 (2017).
- 34 T. Kasimatis, J. Hizanidis, and A. Provata, 'Three-dimensional chimera patterns in networks of spiking neuron oscillators', Phys. Rev. E 97 , 052213 (2018).
- 35 G. Argyropoulos, T. Kasimatis, and A. Provata, 'Chimera patterns and subthreshold oscillations in two-dimensional networks of fractally coupled leaky integrate-and-fire neurons', Phys. Rev. E 99 , 022208 (2019).
- 36 D. Dudkowski, Y. Maistrenko, and T. Kapitaniak, 'Occurrence and stability of chimera states in coupled externally excited oscillators', Chaos 26 , 116306 (2016).
- 37 M. G. Clerc, S. Coulibaly, M. A. Ferré, and R. G. Rojas, 'Chimera states in a Duffing oscillators chain coupled to nearest neighbors', Chaos 28 , 083126 (2018).
- 38 K. K. Likharev, Dynamics of Josephson Junctions and Circuits , Gordon and Breach, Philadelphia (1986).
- 39 J. Hizanidis, N. Lazarides, and G. P. Tsironis, 'Flux bias-controlled chaos and extreme multistability in SQUID oscillators', Chaos 28 , 063117 (2018).
- 40 I. M. Dmitrenko, G. M. Tsoi, V. I. Shnyrkov, and V. V. Kartsovnik, 'RF SQUID in the Nonhysteretic Regime with k 2 Q glyph[lscript] &gt; 1', J. Low Temp. Phys. 49 , 417 (1982).
- 41 T. Banerjee, D. Biswas, D. Ghosh, E. Schöll, and A Zakharova, 'Networks of coupled oscillators: From phase to amplitude chimeras', Chaos 28 , 113124 (2018).
- 42 A. Marchionne, P. Ditlevsen, and S. Wieczorek, 'Synchronisation vs. resonance: Isolated resonances in damped nonlinear oscillators', Physica D 380 , 8 (2018) .
- 43 K. Engelborghs, T. Luzyanina, and D. Roose, 'Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL', ACM Trans. Math. Softw. 28 , 1 (2002).
- 44 F. P. Kemeth, S. W. Haugland, L. Schmidt, I. G. Kevrekidis, and K. Krische, 'A classification scheme for chimera states', Chaos 26 , 094815 (2016).
- 45 Weiming Wang, Quan-Xing Liu, and Zhen Jin, 'Spatiotemporal complexity of a ratio-dependent predator-prey system', Phys. Rev. E 75 , 051913 (2007).
- 46 P. Muruganandam, K. Murali and M. Lakshmanan, 'Spatiotemporal dynamics of coupled array of Murali Lakshmanan Chua circuits', International Journal of Bifurcation and Chaos 9 (5), 805-830 (1999).
- 47 N Lazarides, J Hizanidis, and G. P. Tsironis, 'Controlled Generation of Chimera States in SQUID Metasurfaces using DC Flux Gradients', Chaos, Solitons and Fractals 130 , 109413 (2020).
- 48 N. Lazarides and G. P. Tsironis, 'Multistable dissipative breathers and collective states in SQUID Lieb metamaterials', Phys. Rev. E 98 , 012207 (2018).
- 49 A. P. Zhuravel, S. Bae, A. V. Lukashenko, A. S. Averkin, A. V. Ustinov, and S. M. Anlage, 'Imaging collective behavior in an rf-SQUID metamaterial tuned by DC and RF magnetic fields', Appl. Phys. Lett. 114 , 082601 (2019).
Johanne Hizanidis Nikos Lazarides Giorgos P. Tsironis
Pattern formation and chimera states in 2D SQUID metamaterials
The Superconducting QUantum Interference Device (SQUID) is a highly nonlinear
oscillator with rich dynamical behavior, including chaos. When driven by a
time-periodic magnetic flux, the SQUID exhibits extreme multistability at
frequencies around the geometric resonance which is manifested by a
"snake-like" form of the resonance curve. Repeating motifs of SQUIDs form
metamaterials, i. e. artificially structured media of weakly coupled discrete
elements that exhibit extraordinary properties, e. g. negative diamagnetic
permeability. We report on the emergent collective dynamics in two-dimensional
lattices of coupled SQUID oscillators, which involves a rich menagerie of
spatio-temporal dynamics, including Turing-like patterns and chimera states.
Using Fourier analysis we characterize these patterns and identify
characteristic spatial and temporal periods. In the low coupling limit, the
Turing-like patterns occur near the synchronization-desynchronization
transition which can be related to the bifurcation scenarios of the single
SQUID. Chimeras emerge due to the multistability near the geometric resonance,
and by varying the dc component of the external force we can make them appear
and reappear and, also, control their location. A detailed analysis of the
parameter space reveals the coexistence of Turing-like patterns and chimera
states in our model, as well as the ability to transform between these states
by varying the system parameters.