Multi-branched resonances, chaos through quasiperiodicity, and asymmetric states in a superconducting dimer
MMulti-branched resonances, chaos through quasiperiodicity, andasymmetric states in a superconducting dimer
J. Shena, N. Lazarides, and J. Hizanidis a) National University of Science and Technology “MISiS”, Leninsky Prospekt 4, Moscow 119049,Russia Department of Physics, University of Crete, 71003 Herakleio, Greece (Dated: 15 June 2020)
A system of two identical SQUIDs (superconducting quantum interference devices) symmetrically coupled throughtheir mutual inductance and driven by a sinusoidal field is investigated numerically with respect to dynamical propertiessuch as its multibranched resonance curve, its bifurcation structure, as well as its synchronization behavior. The
SQUIDdimer is found to exhibit a hysteretic resonance curve with a bubble connected to it through Neimark-Sacker (torus)bifurcations, along with coexisting chaotic branches in their vicinity. Interestingly, the transition of the SQUID dimerto chaos occurs through a period-doubling cascade of a two-dimensional torus (quasiperiodicity-to-chaos transition).The chaotic states are identified through the calculated Lyapunov spectrum, and their basins of attraction have beendetermined. Bifurcation diagrams have been constructed on the parameter plane of the coupling strength and thedriving frequency of the applied field, and they are superposed to maps of the maximum Lyapunov exponent on thesame plane. In this way, a clear connection between chaotic behavior and torus bifurcations is revealed. Moreover,asymmetric states that resemble localized synchronization have been detected using the correlation function betweenthe fluxes threading the loop of the SQUIDs. The effect of intermittent chaotic synchronization, which seems to bepresent in the SQUID dimer, is only slightly touched.
Networks of coupled limit cycle oscillators represent aclass of systems with special interest in physics, chemistry,and biological sciences. Arrays of coupled nonlinear os-cillators, in particular, which often exhibit complex dy-namical behavior, have attracted large amounts of theo-retical, numerical, and experimental efforts. However, theprominent features of a finite large number of such os-cillators can be sometimes understood by analysing justtwo coupled oscillators. A unique in many aspects oscil-lator is the superconducting quantum interference device(SQUID) that has been investigated extensively for manyyears. The SQUID is a low-loss, highly nonlinear reso-nant element that responds strongly to applied magneticfield(s) and exhibits rich dynamical behavior. Moreover, ithas recently been the elementary unit for the constructionof metamaterials in one and two dimensions. SQUIDs andSQUID arrays are technologically important devices, andthey also serve as a testbed for exploring complex dynam-ics. Two SQUIDs in close proximity are coupled togetherthrough magnetic dipole-dipole forces, and their dynam-ical complexity significantly increases. Specifically, morecomplex bifurcation scenarios appear, along with transi-tions to chaos through quasiperiodicity, the emergence oflocalized synchronization, and intermittent chaotic syn-chronization. These dynamical properties are explorednumerically with a well-established model whose parame-ters are acquired from recent relevant experiments. a) Electronic mail: [email protected]
I. INTRODUCTION
There has been great interest on the behavior of systemswith interacting and forced nonlinear oscillations with ap-plications in physics, biology, and chemistry, which in ad-dition are in abundance in the natural world . The in-teraction introduces qualitatively new behavior to the net-work of oscillators, as compared to that of a single os-cillator. For example, synchronization of the oscillat-ing units of the system , extreme multistability , multi-ple resonance and anti-resonance effects , complex bifur-cation structure , localization , emergence of exceptionalpoints , the intriguing effect of amplitude death , chaossynchronization , localized synchronization , chaos to hy-perchaos transitions , and quasi-periodicity with subse-quent transition to chaos , are some of the experimentallyor numerically observed behaviors.Here we consider a pair of coupled nonlinear oscillators(i.e., a nonlinear dimer) that belong to the large class of ex-ternally driven and dissipative systems. More specifically,we consider a pair of identical mutually coupled SQUIDoscillators , where the acronym stands for “superconduct-ing quantum interference device”. SQUIDs are mesoscopicsuperconducting devices which are modeled efficiently byequivalent electrical circuits, while their dynamics is governedby a second order nonlinear ordinary differential equation.The simplest variant of a SQUID consists of a superconduct-ing loop interrupted by a Josephson junction (JJ), called rfSQUID, which was found to exhibit very rich dynamic be-havior, including complex bifurcation structure and chaos .Moreover, rf SQUIDs are employed as elementary units forthe fabrication of superconducting metamaterials, whose in-vestigation has revealed several extraordinary properties .Recent theoretical works on SQUID metamaterials have alsoreported the existence of chimera states and pattern forma- a r X i v : . [ n li n . C D ] J un ition . FIG. 1. Schematic diagram of a SQUID dimer in a magnetic fieldwhere “Mf” is the Magnetic field, “Sr” is the Superconducting ring,“JJ” is the Josephson Junction, and I , I are the induced currents.Panels (a) and (b) correspond to negative and positive magnetic cou-pling strength, respectively, appropriate for the planar and the axialgeometry. Two SQUIDs in relatively close proximity are coupled to-gether by magnetic dipole-dipole forces due to their mutualinductance M , while the sign and the magnitude of their cou-pling strength coefficient depends strongly on their relativepositions. Specifically, the two SQUIDs may be arranged ei-ther in the planar geometry (Fig. 1(a)), where they lie in thesame plane with their axes in parallel, or in the axial geome-try, they lie the one on top of the other with their axes lyingon the same line (Fig. 1(b)). The SQUIDs are driven by ex-ternally applied periodic and constant (bias) magnetic fieldsof appropriate polarization, to which they respond resonantly.Although its potential for technological applications, the dy-namics of a SQUID dimer are not properly investigated. Re-cently however, the existence of homoclinic chaos in a pair ofparametrically driven SQUIDs has been exposed .In the present work, particular dynamical properties of aSQUID dimer are addressed, such as its highly unusual res-onance curve that consists of a bistable part with a bubbleconnected to it through Neimark-Sacker (NS) and pitchforkbifurcations (PBs). These bubbles are of very different char-acter compared to those reported long ago by Bier and Boun-tis and they deserve to be investigated on their own. More-over, in addition to those solutions, there are at least twomore coexisting solutions in the frequency interval spanned bythe resonance curve, which end up becoming chaotic throughquasiperiodicity, i.e., through period doubling bifurcations oftwo-dimensional (2D) tori. The basins of attraction lead-ing to chaos are calculated for two parameter sets of partic-ular importance. Furthermore, asymmetric solutions eithersynchronized (leading to localized synchronization) or anti-synchronized (i.e., having opposite phases) are investigatedby calculating the maximum Lyapunov exponent and the cor-relation function over the entire control parameter space.In the next section (II), the electrical circuit equivalentmodel for the SQUID dimer is described in detail, and thedynamic equations are obtained. In section III, the bifurca-tion structure for the SQUID dimer is discussed along withthe emergence of quasiperiodicity and chaos. Section IV is devoted to the transition to chaos through quasiperiodicity,through a Neimark-Sacker (NS) bifurcation. The existenceof highly asymmetric solutions related to localized synchro-nization is revealed using an appropriate correlation functionin section V. Conclusions are given in section VI. II. MODELING TWO COUPLED SQUIDS
Consider two identical rf SQUIDs in close proximity so thatthey are coupled by magnetic dipole-dipole forces throughtheir mutual inductance M (Fig. 1). Each of the SQUIDsis modeled by an equivalent electrical circuit that features aself-inductance L due to the superconducting ring, which isconnected in series with a “real” Josephson junction charac-terized by a critical current I c , capacitance C , and Ohmic re-sistance R . The SQUIDs are subject to an externally appliedspatially constant and time-periodic (ac) magnetic field and aconstant in time and space (dc) magnetic field. Those fieldsadd an electromotive force in series to the SQUID’s equivalentcircuits. The magnetic flux threading the loops of the SQUIDsincludes supercurrents as well as normal (i.e., quasiparticle)currents around the SQUID’s rings through Faraday’s law. Inturn, the induced currents produce their own magnetics fieldwhich counter-acts the applied ones. Then, the flux Φ and Φ through the loop of the SQUID number 1 and 2, respectively,is given by the following flux-balance relations Φ = Φ ext + L I + M I , (1) Φ = Φ ext + L I + M I , (2)where I and I are the currents induced by the external mag-netic fields, and the external flux Φ ext = Φ dc + Φ ac cos ( ω t ) , (3)with Φ dc the dc flux bias, Φ ac the amplitude of the ac flux, ω its frequency, and t the temporal variable. Eqs. (1) and (2) canbe solved for the currents and written in matrix form as (cid:20) I I (cid:21) = L (cid:20) λλ (cid:21) − (cid:20) Φ − Φ ext Φ − Φ ext (cid:21) , (4)where λ = M L is the magnetic coupling strength between theSQUIDs. Note that the sign of λ depends on the mutualposition of the two SQUIDs. In the planar geometry (Fig.1(a)), that sign is negative, while in the axial geometry (Fig.1(b)), it is positive. The currents flowing in the SQUIDs aregiven within the framework of the resistively and capacitivelyshunted junction (RCSJ) model , as I n = − C d Φ n dt − R d Φ n dt − I c sin (cid:18) π Φ n Φ (cid:19) , (5)where Φ is the flux quantum, and n = ,
2. Combining Eqs.(4) and (5) we get the coupled equations for the fluxes throughii
FIG. 2. (a) The potential U SQ ( φ , φ ) of the SQUID dimer for τ = T = π / Ω where Ω = . U SQ ( φ , φ ) on the φ = λ = . φ ac = . γ = . β = . the loops of the SQUIDs in natural units, as LC d Φ dt + LR d Φ dt + L I c sin (cid:18) π Φ Φ (cid:19) + − λ Φ − λ − λ Φ = + λ Φ ext , (6) LC d Φ dt + LR d Φ dt + L I c sin (cid:18) π Φ Φ (cid:19) + − λ Φ − λ − λ Φ = + λ Φ ext , (7)In normalized form, Eqs. (6) and (7) can be written as¨ φ + γ ˙ φ + β sin ( πφ ) + − λ φ − λ − λ φ = + λ φ ext , (8)¨ φ + γ ˙ φ + β sin ( πφ ) + − λ φ − λ − λ φ = + λ φ ext , (9)with φ ext = φ dc + φ ac cos ( Ω τ ) , (10)where all the fluxes Φ , Φ , Φ ac , Φ dc are in units of the fluxquantum Φ and the new temporal variable τ is defined as τ = t ω LC with ω LC = / √ LC being the inductive-capacitive( LC ) SQUID frequency. Note that the normalized driving fre-quency Ω is in units of ω LC , i.e., Ω = ω / ω LC . The overdots inEqs. (8) and (9) denote differentiation with respect to the nor-malized temporal variable τ . The rescaled SQUID parameterand the loss coefficient are given respectively by β = L I c Φ , γ = R (cid:114) LC . (11)In the lossless case, i.e., when γ = R → ∞ ), the (Newtonian)dynamics of the SQUID dimer is governed by the equations¨ φ = − ∂ U SQ ( φ , φ ) ∂ φ , ¨ φ = − ∂ U SQ ( φ , φ ) ∂ φ (12) in terms of the SQUID dimer potential U SQ ( φ , φ ) = − + λ ( φ + φ ) φ ext ( τ ) − β π [ cos ( πφ ) + cos ( πφ )]+ ( − λ ) (cid:0) φ − λ φ φ + φ (cid:1) . (13)Since the SQUID dimer considered here consists of twoidentical non-hysteretic SQUIDs (2 πβ < U SQ ( φ , φ ) is a two-dimensional corrugated parabola with asingle minimum that depends periodically on time τ . A snap-shot of the potential at τ = T , with T = π / Ω being the periodof the (driving) external periodic field, is shown in Fig. 2(a)while its one-dimensional projection on the φ = φ dc =
0. For obtaining the values of the SQUID-dependent normal-ized parameters β and γ that go into the model equations (8)and (9), we choose L = pH , C = . pF , R = Ω ,and I c = . µ A for their equivalent lumped-circuit elements.These values are typical for non-hysteretic SQUIDs investi-gated recently in the context of SQUID metamaterials .By substituting them into Eq. (11), we get β L = .
86 ( β (cid:39) . γ = . LC or geometrical fre-quency is f LC = ω LC / ( π ) (cid:39) . GHz . As mentioned earlier,the value of the coupling strength λ depends significantly onthe relative positions of the two SQUIDs, and may assumeits value from the relatively broad range of [ − . , + . ] that corresponds to technologically feasible SQUID dimerdesigns . Besides, there are also the parameters φ ac and Ω ( φ dc = φ ac = .
02, which is well into therange of the experimentally accessible values, and use λ and Ω as control parameters. III. BIFURCATIONS, QUASI-PERIODICITY, AND CHAOS
It is known that for relatively high φ ac , the single SQUIDexhibits a multistable response which is reflected in its cor-responding “snake-like” resonance curve , i.e., the curveformed by the maximum flux within one driving period T = π / Ω , φ max , as a function of the driving frequency Ω . Forthe value of φ ac = .
02 considered here, the single SQUID isbistable. When two SQUIDs are coupled positively (i.e., when λ >
0) the bistability is maintained but the dynamics presentsadditional complexity. The resonance curve for the SQUIDdimer is plotted in Fig. 3(a) in terms of the magnetic flux φ threading the loop of the SQUID number 1. Specifically,the maximum flux threading the loop of the SQUID number1 within one driving period T , φ , max , is plotted against thedriving frequency Ω . We see the typical tilted curve asso-ciated with bistability and hysteresis, where two stable solu-tions (black curves) coexist with an unstable one (red curve)below the geometric resonance frequency Ω =
1. Starting atlower driving frequencies and following the periodic solutionv
FIG. 3. (a) The resonance curve of the SQUID dimer. Black lines correspond to stable branches and red and blue lines correspond to unstableones. The latter are born through saddle-node bifurcations of limit cycles (SN), Neimark-Sacker bifurcations (NS), or pitchfork bifurcationsof limit cycles (PB). (b) Maximum values of the magnetic flux φ (top) and all five Lyapunov exponents (bottom), as a function of the drivingfrequency in the interval corresponding to the left blue branches of (a). For each value of Ω , n i =
50 different initial conditions were used. (c)The same as in (b) for the right blue branches of (a). Other parameters: λ = . φ ac = . γ = . β = . in Ω there is a saddle-node bifurcation of limit cycles (SN) oc-curring at Ω = .
2. The solution then becomes unstable andwhen it reaches the peak of the resonance curve it turns stableagain in a second SN bifurcation at Ω = . Ω = . Ω = . Ω = . Ω = .
12 (right PB)and this family of solutions emanating from the two “oppositeto each other” PB bifurcations, undergoes NS bifurcations (at Ω = .
083 and Ω = . .We focus now on the two Ω intervals corresponding to theblue unstable branches of Fig. 3(a) which lie within the twoPB bifurcations on both right and left side. In the upper panelsof Figs. 3(b) and (c) the maxima of the solution for φ ( τ ) areplotted over Ω for the left and right interval, respectively. Notethat at each frequency Ω , a large number of different initialconditions n i were used (here n i =
50) in order to obtain allthe coexisting solutions within the Ω interval shown, and alsoto obtain a good representation of possible chaotic behavior.Apparently, chaotic behavior appears in both the upper panelsof Figs. 3(b) and (c). These solution branches coexist withthe blue unstable branches of Fig. 3(a) mentioned above; thereason for which they are not superposed to them is merely clarity.In the corresponding lower panels, the Lyapunov exponentsof the obtained solutions that are plotted for the same intervalsof Ω , reveal chaotic behavior in the regions where the maxi-mum Lyapunov exponent L is positive (positive Lyapunovexponents are shown as blue points). Careful inspection of theLyapunov exponents indicate that transitions to chaos throughquasiperiodicity and reversely take place in Figs. 3(b) and (c)(see next section). The Lyapunov exponents were calculatedemploying the algorithm from using the Julia R software li-brary. Similar multibranched resonance curve are sometimesencountered in dissipative-driven oscillators where in ad-dition, isolated branches (“isolas”) may exist. However, tothe best of our knowledge, multibranched resonance curveswhose branches are connected through Neimark-Sacker bifur-cations have never before been observed.Since the system dynamics is very sensitive to initial con-ditions, we have identified the basins of attraction of chaoticand periodic behavior for two Ω values within the intervalsshown in Figs. 3(b) and (c) close to the NS bifurcations oflimit cycles. These are displayed in Figs. 4(a) and (b), re-spectively, both in the ( φ , φ ) plane (left panels) and ( ˙ φ , ˙ φ ) plane (right panels). Blue (dark) regions denote the set of ini-tial conditions leading to chaotic motion, while light yellow(light) regions the ones leading to periodic or quasiperiodicoscillations. The other state variables are initially fixed as˙ φ ( τ = ) = ˙ φ ( τ = ) = φ ( τ = ) = φ ( τ = ) = FIG. 4. Basins of attraction for the SQUID dimer in terms of themaximum Lyapunov exponent L max . Blue and yellow color cor-responds to chaotic and non-chaotic dynamics, respectively. Left: φ versus φ for ( ˙ φ ( ) , ˙ φ ( )) = ( , ) . Right: ˙ φ versus ˙ φ for ( φ ( ) , φ ( )) = ( , ) . (a) Ω = . Ω = .
1. The inset inthe left panel of (a) shows a blow-up of the region inside the red box.Other parameters: λ = . φ ac = . γ = . β = . both subfigures, that resemble riddled basins investigated invarious low-dimensional systems . The inset of Fig. 4(a)reveals this intricate structure of the basin of attraction on asmaller scale. Note that a classification of basins of attrac-tion in several low-dimensional systems has been perormed inRef. . The knowledge of the basins is essential for determin-ing the usefulness of the systems in practical applications.The other control parameter in our system is the couplingstrength λ , which (as described in Section II) can obtainboth negative and positive values within a certain range pro-vided by calculations on feasible SQUID dimer designs . InFig. 6(a) the co-dimension 2 bifurcation diagram is shown inthe ( Ω , λ ) parameter plane. Black lines mark the continua-tion of the saddle-node bifurcations and red lines that of theNeimark-Sacker (torus) bifurcations. In Fig. 6(b) the maxi-mum Lyapunov exponent is plotted in the ( Ω , λ ) parameterspace. Blue and yellow regions correspond to chaotic (posi-tive L max ) and non-chaotic (negative or zero L max ) dynamics,respectively. We can see that the chaotic regimes lie within theboundaries of the NS bifurcation curves. In the following sec-tion we explore in detail the transition from quasiperiodicityto chaos in our system. FIG. 5. (a) Bifurcation diagram in the ( Ω , λ ) parameter plane of thetwo coupled SQUIDs. Black and red lines correspond to saddle-node(SN) bifurcations of limit cycles and Neimark-Sacker (NS) bifurca-tions, respectively. (b) Map of the maximum Lyapunov exponent L max in the ( Ω , λ ) parameter space, on which the NS bifurcationcurve from (a) is superposed (red lines). Blue and yellow color cor-responds to chaotic and non-chaotic dynamics, respectively. Otherparameters: φ ac = . γ = . β = . IV. TORUS-DOUBLING TRANSITION TO CHAOS
From Figs. 3(b)-(c) and 5(b) it is clear that the SQUIDdimer described by Eqs. (8) and (9) becomes chaotic throughNeimark-Sacker bifurcations. This quasiperiodic transitionto chaos is a well-known scenario whereby a torus in low-dimensional dynamical systems loses its stability and devel-ops into chaos. A path from a limit cycle, then a transitioninto 2D torus following by a 3D torus and eventually directtransition to chaos has been presented numerically and exper-imentally is an analog electrical circuit representing the ringof unidirectionally coupled single-well Duffing oscillators .Another possible scenario to chaos is a cascade of period dou-blings of the torus which has been investigated numericallylong ago using low-dimensional maps , and experimen-tally in electrochemical reactions . According to the secondscenario, a finite cascade of period doubling of such invari-ant tori leads eventually to chaos. This scenario has been re-ported numerically for a quintic complex Ginzburg-Landauequation and bimodal laser model as well as experimen-tally, in Rayleigh-Benard convection , in a simple thermoa-coustic system , and near the ferroelectric phase transition of KH PO crystals .We examine this route to chaos through period doublingof a 2D torus in our system by fixing the magnetic couplingstrength to a certain value λ = .
15 (marked by a black hori-zontal line in Fig. 5(b)) and by varying the driving frequency Ω in a narrow interval that encloses chaotic states.For this choice of parameters, the resonance curve(Fig. 6(a)) exhibits the same sequence of bifurcations (PBand consequently NS) as in Fig. 3(a), but presents no secondchaotic band and therefore no “bubble” connected to it. Thetransition into chaos is shown in Fig. 6(b) where the maxi-mum values of the magnetic flux φ are plotted for the same Ω interval.For each value of Ω the whole spectrum of the Lyapunoviexponents is calculated and plotted in Fig. 7 (top). Startingfrom the right end of the frequency interval, it is observedthat the two largest Lyapunov exponents are zero ( L = L = Ω decreases, one more Lyapunov exponent reacheszero at Ω = . Ω , we find three more points in which the three largestLyapunov exponents are zero at Ω = . . . FIG. 6. (a) The resonance curve of the SQUID dimer for λ = . φ of stable solutions as a function of the driv-ing frequency in the same Ω interval. Other parameters: φ ac = . γ = . β = . L i , i = ,..., Ω . The bifurcationanalysis depicts a transition from a stable two-dimensional torus tochaos via a sequence of four consecutive period-doubling bifurca-tions marked by P i ( i = , , , λ = . φ ac = . γ = .
024 and β = . torus break down and the system enters into a chaotic state, asit is indicated by the positivity of the largest Lyapunov expo-nent ( L = L max > d L , which according to Kaplanand York is given by: d L = k + | L k + | k ∑ i = L i , (14)where k is defined by the condition that k ∑ i = L i ≥ k + ∑ i = L i < . (15)As the 2D torus transitions into chaos, the Lyapunov dimen-sion changes from the value d L = d L >
3, as shown inFig. 7 (bottom).Typical Poincaré sections of the period-2, -4, -8, andchaotic 2D tori on the φ , − ˙ φ , phase plane are shown inthe upper panels of Figs. 8(a) - (e). Blue and red lines corre-spond to flux threading the loop of the first (1) and the second(2) SQUID, respectively. The insets show enlargements of asmall part of the Poincaré sections of the first SQUID (SQUIDnumber 1), in order to see clearly the varying number of layersof the tori as the frequency Ω decreases. The Poincaré sectionof the period-16 torus is very similar to that of the period-8 torus and thus it has been omitted. In the lower panels ofFig. 8, the corresponding Fourier power spectral densities ina logarithmic scale are shown. The 2D torus (upper panel ofFig. 8(a)) has a dominant frequency which in the present casecoincides with the external (driving) frequency Ω , as well asfour more frequencies denoted as Ω ± , Ω ± (lower panel ofFig. 8(a)). The exact values of these frequencies are obtainedfrom the locations of the sharp peaks. Then, as the period ofthe 2D torus increases to 2, 4, and 8 through period-doublingbifurcations, a doubling of the sharp peaks occurs (as can beseen in Fig. 8(b), (c), and (d), respectively) that correspondto frequencies with significant spectral content. In Fig. 8(e),in which the SQUID dimer is in a chaotic state, the powerspectrum exhibits a substantial noisy background, along withseveral sharp peaks that correspond to frequencies with sig-nificant spectral content. Note that the driving frequency Ω isdominant in all five subfigures.Thus, we have identified a novel route to chaos for theSQUID dimer, which initiates from a single limit cycle. Thelatter then bifurcates into a 2D torus, which follows a finiteperiod doubling cascade which eventually leads to a chaoticstate. Note that in the case of two coupled Duffing oscillatorswith dissipation and driving force, three period-doubling bi-furcations of a three-torus were identified for the transition tochaos through quasiperiodicity . V. LOCALIZED & CHAOTIC SYNCHRONIZATION
The SQUID dimer considered here consists of two sym-metrically coupled identical nonlinear oscillators. Such sys-tems were once believed to support only totally synchronousii
FIG. 8. Poincaré cross sections on the Φ − ˙ Φ plane (upper panels) and the corresponding frequency spectra plotted in a logarithmic scale(bottom panels), for (a) Ω = . Ω = .
27; (c) Ω = . Ω = . Ω = . Φ = φ , ˙ Φ = ˙ φ ) and second ( Φ = φ , ˙ Φ = ˙ φ ) SQUID, respectively. Insets:
Enlargement of a small part of Poincarésections of the trajectories of the first SQUID to make clear the period of the tori. Other parameters as in Fig. 7. or asynchronous states. However, as it is known today, thereare more possible forms of synchronization such as the socalled localized synchronization. This form of synchroniza-tion has been originally investigated in a slightly dissimilarpair of solid state and semiconductor lasers. For two cou-pled oscillators, localized synchronization means that one ofthe oscillators exhibits strong oscillations while the other oneexhibits weak oscillations. In the other limit, that of asyn-chronous dynamics, extreme case of having the one oscillatorexhibiting periodic motion and the other chaotic motion hasbeen recently reported . We explore below the effect of lo-calized synchronization for the SQUID dimer using the cor-relation function between the time-dependent fluxes throughthe loop of the first and the second SQUID φ ( t ) and φ ( t ) ,respectively, i.e., C = (cid:104) ( φ ( t ) − µ )( φ ( t ) − µ ) (cid:105) σ σ , (16)where µ , and σ , are the mean value and the standard de-viation of the corresponding time series of φ and φ , respec-tively. Using this measure we can quantify the asymmetry inthe amplitudes of the two coupled SQUID magnetic fluxes. If φ and φ are increasing or decreasing simultaneously then thecorrelation function will be positive. In contrast, if the valueof φ is increasing while φ is decreasing and vice versa, then C is negative. The maximum/minimum value of the corre-lation function is ±
1. For C (cid:39) φ and φ are almost identical thereby the two SQUIDs exhibit in-phase synchronization, whereas for C (cid:39) − C is mapped onto the Ω , λ parameter space. As men-tioned earlier, the value of C reflects the dynamical behaviorof the SQUID dimer as long as its synchronization propertiesare concerned. As it can be observed, the value of C is equalto 1 in a large area (yellow) of the parameter space but therealso exist regions where C (cid:54) =
1. For example, a dark-blueregion (I) can be observed in which C approaches −
1, in-dicating anti-synchronized (i.e., phase difference π between φ ( t ) and φ ( t ) ) dynamics. Moreover, we observe a “horn”-shaped region of intermediate to high values of C which co-incides with the corresponding “horn”-shaped chaotic regionof Fig. 6(b). This particular region, which is bounded by theNeimark-Sacker (NS) bifurcation lines, intermittent chaoticsynchronization is observed. In the case of the SQUIDdimer, the effect of intermittent chaotic synchronization ismanifested by the entrainment of the two SQUID time seriesfor φ ( t ) and φ ( t ) in random temporal intervals of finite du-ration . This will be further investigated in future works in-volving SQUID trimers where more interesting synchroniza-tion phenomena can be observed, like small chimeras .Next we focus on the behavior of the correlation func-tion C in the non-chaotic region of the ( Ω , λ ) parameterspace and reveal another interesting behavior. For a fixed Ω we study the dependence of C as the coupling strength λ varies from region I to II of Fig. 9(a). A cross section ofiii C for Ω = . n i =
30 sets ofdifferent, randomly chosen initial conditions (bottom). Thedifference between the maximum values of the two SQUIDmagnetic fluxes φ − φ is also plotted in black color.We make the following observations: first of all, the fullyphase-synchronized state ( C =
1) exists for all λ in this in-terval. Secondly, the correlation function appears to have a(co)sinusoidal dependence on the coupling strength λ as ex-pected, since the correlation function of two out-of-phase co-sine signals is a cosine function of their phase difference . Inthe low-coupling non-chaotic regions I and II we can assumethat the two SQUID solutions are also (co)sine-like. For φ ( t ) and φ ( t + τ ) or vice versa for φ ( t + τ ) and φ ( t ) , τ is thephase difference between the SQUID fluxes and is dependenton the set of initial conditions. From Fig. 9(b) (bottom) wemay claim that the τ = C = C corresponds to the case where τ is a function of λ ) . In the former case the two SQUIDs arein phase and φ − φ is zero, while in the latter case theirphase difference varies from values close to − + φ − φ attains a nonzero positive/negative valuedepending on whether φ is larger/smaller than φ , for all λ and all initial conditions. Therefore this case corresponds toan asymmetric solution as far as the amplitudes of the twocoupled SQUIDs are concerned.This is illustrated in Fig. 9(c) in the phase portraits (top) andtime series (bottom) of the magnetic flux for the first SQUID(green line) and the second SQUID (blue line) for parameterstha belong to region I ( λ = − . , C = − . ) and re-gion II ( λ = . , C = . ) . These plots show thephase portraits of the corresponding stable periodic solutionswhere large asymmetry in the amplitudes of the two coupledSQUIDs occurs. Due to that difference between the ampli-tudes, these two states are illustrative examples of localizedsynchronization. Moreover, through the time series, it can beclearly seen that the fluxes ( φ and φ ) are almost out of phase(I) and almost in phase (II) because the correlation functionhas values close to − VI. CONCLUSIONS
In summary, a dimer comprising two magnetically cou-pled SQUIDs is investigated numerically with respect to itsbifurcation structure and its transitions from quasiperiodic-ity to chaos through a finite period-doubling cascade of tori.The SQUID dimer is a prototype, dissipative-driven dynam-ical system of considerably complex dynamic behavior, onwhich numerous nonlinear dynamic effects can be investi-gated. In the present work, the bifurcation structure of theSQUID dimer was obtained using the driving frequency Ω andthe coupling strength λ as external control parameters thatform a two-dimensional parameter space. Bifurcation (stro-boscopic) diagrams as a function of Ω reveal multibranchedresonance curves resulting through different kinds of bifurca-tions, transitions to chaos after four period-doubling bifurca- tions of a 2D torus, and riddle-like basins. Bifurcation dia-grams for SN and NS bifurcations are shown on the Ω − λ plane. In combination with maps of the maximum Lyapunovexponent Λ max , it is observed that the areas of the Ω − λ planewhere chaotic behavior is expected is bounded by the NS bi-furcation curves. Eventually, localized synchronization is ob-served and quantified in the SQUID dimer using the correla-tion function C . The possibility of intermittent chaotic syn-chronization is briefly mentioned. ACKNOWLEDGMENTS
This work was supported by the General Secretariat for Re-search and Technology (GSRT) and the Hellenic Foundationfor Research and Innovation (HFRI) (Code: 203). JS ac-knowledges support by the Ministry of Education and Scienceof the Russian Federation in the framework of the IncreaseCompetitiveness Program of NUST “MISiS” (Grant numberK4-2018-049). The authors would like to thank S. M. Anlagefor fruitful discussions. Jan Awrejcewicz, “
Bifurcation And Chaos In Coupled Oscillators ,” WorldScientific, Singapore (1991). M. Rosenblum and A. Pikovsky, “
Synchronization: from pendulum clocksto chaotic lasers and chemical oscillators ,” Contemporary Physics
44 (5) ,401-416 (2003). K. Wiesenfeld and P. Hadley, “
Attractor crowding in oscillator arrays ,”Phys. Rev. Lett. , 1335-1338 (1989). Y. Nishio and S. Mori, “
Mutually coupled oscillators with an extremelylarge number of steady states ,” IEEE Xplore, 819-822 (1992). R. Jothimurugan, K. Thamilmaran, S. Rajasekar, and M. A. F. Sanjuán,“
Multiple resonance and anti-resonance in coupled Duffing oscillators ,”Nonlinear Dyn. Y. Kominis, K. D. Choquette, A. Bountis, and V. Kovanis, “
Antiresonancesand Ultrafast Resonances in a Twin Photonic Oscillator ,” IEEE PhotonicsJournal
11 (1) , 1500209 (2019). Prasun Sarkar and Deb Shankar Ray, “
Vibrational antiresonance in nonlin-ear coupled systems ,” Phys. Rev. E , 052221 (2019). Hua-Wei Yin, Jian-Hua Dai, and Hong-Jun Zhang, “
Phase effect of twocoupled periodically driven Duffing oscillators ,” Phys. Rev. E
58 (5) , 5683-5688 (1998). A. Kenfack, “
Bifurcation structure of two coupled periodically drivendouble-well Duffing oscillators ,” Chaos Soliton Fract. , 205-218 (2003). N. Lazarides, M. I. Molina, G. P. Tsironis, and Yu. S. Kivshar, “
Multistabil-ity and localization in coupled nonlinear split-ring resonators ,” Phys. Lett.A , 2095-2097 (2010). Y. Kominis, K. D. Choquette, A. Bountis, and V. Kovanis, “
Exceptionalpoints in two dissimilar coupled diode lasers ,” Appl. Phys. Lett. ,081103 (2018). E. Tafo Wembe and R. Yamapi, “
Chaos synchronization of resistively cou-pled Duffing systems: Numerical and experimental investigations ,” Com-mun. Nonlinear Sci. Numer. Simul. , 1439-1453 (2009). R. Herrero, M. Figueras, J. Rius, F. Pi, and G. Orriols, “
Experimental obser-vation of the amplitude death effect in two coupled nonlinear oscillators ,”Phys. Rev. Lett.
84 (23) , 5312-5315 (2000). A. Hohl, A. Gavrielides, T. Erneux, and V. Kovanis, “
Localized synchro-nization in two coupled nonidentical semiconductor lasers ,” Phys. Rev.Lett.
78 (25) , 4745-4748 (1997). T. Kapitaniak and W.-H. Steeb, “
Transition to hyperchaos in coupled gen-eralized van der Pol equations ,” Phys. Lett. A
152 (1,2) , 33-36 (1991). T. Kapitaniak, K.-E. Thylwe, I. Cohen, and J. Wojewoda, “
Chaos-Hyperchaos Transition ,” Chaos Soliton Fract. , 2003-2011 (1995). T. Kapitaniak, Y. Maistrenko, and S. Popovych, “
Chaos-hyperchaos transi-tion ,” Phys. Rev. E
62 (2) , 1972-1976 (2000). x FIG. 9. (a) Dynamical behavior in the ( Ω , λ ) parametric space of the correlation function C . For a completely in-phase synchronized state C = C =
1. (b) The correlation function (red line) and maxima differences between φ and φ (black line), versus the coupling strength ( λ ), for one random set of initial conditions (top) and for n i =
30 different initial conditions(bottom) where Ω = . φ (green line), and thesecond SQUID, φ (blue line), for (I) λ = − . C = − . λ = . C = . φ ac = . γ = . β = . K. Grygiel and P. Szlachetka, “
Chaos and hyperchaos in coupled Kerr os-cillators ,” Optics Communications , 425-431 (2000). Prof. S. M. Anlage (Personal communication) and also https://snf.ieeecsc.org/sites/ieeecsc.org/files/documents/snf/abstracts/Anlage_STP647.pdf , p. 33. D. Rand, S. Ostlund, J. Sethna, and E. D. Siggia, “
Universal transitionfrom quasiperiodicity to chaos in dissipative systems ,” Phys. Rev. Lett. , 132-135 (1982). A. R. Bishop, M. G. Forest, D. W. McLaughlin, E. A. Overman II, “
A quasi-periodic route to chaos in a near-integrable PDE ,” Physica , 293-328(1986). T. W. Dixon, T. Gherghetta, and B. G. Kenny, “
Universality in thequasiperiodic route to chaos ,” Chaos , 32-42 (1996). Z. Elhadj and J. C. Sprott, “
A minimal 2-D quadratic map with quasiperi-odic route to chaos ,” Int. J. Bifurc. Chaos
18 (5) , 1567-1577 (2008). L. Borkowski, P. Perlikowski, T. Kapitaniak, and A. Stefanski, “
Experi-mental observation of three-frequency quasiperiodic solution in a ring ofunidirectionally coupled oscillators ,” Phys. Rev. E , 062906 (2015). J. Kozlowski, U. Parlitz, and W. Lauterborn, “
Bifurcation analysis of twocoupled periodically driven Duffing oscillators ,” Phys. Rev. E , 1861-1867 (1995). R. Kleiner and D. Koelle and F. Ludwig and J. Clarke, “
Superconductingquantum interference devices: State of the art and applications ,” Proceed-ings of the IEEE , 1534-1548 (2004). R. L. Fagaly, “
Superconducting quantum interference device instrumentsand applications ,” Rev. Sci. Instrum. , 101101 (2006). J. Hizanidis, N. Lazarides, and G.P. Tsironis, “
Flux bias-controlled chaosand extreme multistability in SQUID oscillators ,” Chaos , 063117 (2018). M. Trepanier, D. Zhang, O. Mukhanov, and S. M. Anlage, “
Realization andmodeling of RF superconducting quantum interference device metamateri-als ,” Phys. Rev. X , 041029 (2013). D. Zhang, M. Trepanier, O. Mukhanov, and S. M. Anlage, “
Broadbandtransparency of macroscopic quantum superconducting metamaterials ,”Phys Rev X , 041045 (2015). N. Lazarides, G. Neofotistos, and G. P. Tsironis, “
Chimeras in SQUIDMetamaterials ,”, Phys. Rev. B , 054303 (2015). J. Hizanidis, N. Lazarides, and G. P. Tsironis, “
Robust chimera states inSQUID metamaterials with local interactions ”, Phys. Rev. E , 032219(2016). J. Hizanidis, N. Lazarides, G. Neofotistos, and G. P. Tsironis, “
Chimerastates and synchronization in magnetically driven SQUID metamaterials ”,Eur. Phys. J. Special Topics, Springer , 1231 (2016). J. Hizanidis, N. Lazarides, and G. P. Tsironis, “
Pattern formation andchimera states in 2D SQUID metamaterials ”, Chaos , 013115 (2020). M. Agaoglou, V. M. Rothos, and H. Susanto, “
Homoclinic chaos in a pairof parametrically-driven coupled SQUIDs ,” J. Phys. Conf. Ser. , 012027(2015). M. Agaoglou, V. M. Rothos, and H. Susanto, “
Homoclinic chaos in coupledSQUIDs ,” Chaos Soliton Fract. , 133-140 (2017). M. Bier and T. C. Bountis, “
Remerging Feigenbaum Trees in DynamicalSystems ,” Phys. Lett. , 239 (1984). B. D. Josephson, “
Possible new effects in superconductive tunnelling ,”Phys. Lett. A , 251-253 (1962). K. K. Likharev, “
Dynamics of Josephson Junctions and Circuits ,” Gordonand Breach, Philadelphia (1986). K. Engelborghs, T. Luzyanina, and D. Roose, “
Numerical bifurcation anal-ysis of delay differential equations using DDE-BIFTOOL ,” ACM Trans.Math. Softw. , 1 (2002). K. Geist, U. Parlitz, and W. Lauterborn, “
Comparison of different methodsfor computing Lyapunov exponents ,” Prog. Theor. Phys.
83 (5) , 875-893(1990). J. Warminski, “
Frequency locking in a nonlinear MEMS oscillator drivenby harmonic force and time delay ,” Int. J. Dynam. Control , 122-136(2015). A. Marchionne, P. Ditlevsen, and S. Wieczorek, “
Synchronisation vs. res-onance: Isolated resonances in damped nonlinear oscillators ,” Physica D , 8-16 (2018). Jian Zang and Ye-Wei Zhang, “
Responses and bifurcations of a struc-ture with a lever-type nonlinear energy sink ,” Nonlinear Dyn. , 889-906(2019). Ying-Cheng Lai and R. L. Winslow, “
Riddled Parameter Space in Spa-tiotemporal Chaotic DynamicalSystem ,” Phys. Rev. Lett.
72 (11) , 1640- K. Sathiyadevi, S. Karthiga, V. K. Chandrasekar, D. V. Senthilkumar, andM. Lakshmanan, “
Frustration induced transient chaos, fractal and riddledbasins in coupled limit cycle oscillators ,” Commun. Nonlinear Sci. Numer.Simulat. , 586-599 (2019). J. C. Sprott and Anda Xiong, “
Classifying and quantifying basins of attrac-tion ,” Chaos , 083101 (2015). A. Arneodo, P. H. Coullet, and E. A. Spiegel, “
Cascade of period doublingtori ,” Phys. Lett. A
94 (1) , 1-6 (1983). K. Kaneko, “
Oscillation and doubling of torus ,” Prog. Theor. Phys.
72 (2) ,202-215 (1984). M. R. Bassett and J. L. Hudson, “
Experimental evidence of period doublingof tori during an electrochemical reaction ,” Physica D , 289-298 (1989). Jung-Im Kim, Hwa-Kyun Park, and Hie-Tae Moon, “
Period doubling of atorus: Chaotic breathing of a localized wave ,” Phys. Rev. E
55 (4) , 3948-3951 (1997). C. Letellier, M. Bennoud, and G. Martel, “
Intermittency and period-doubling cascade on tori in a bimode laser model ,” Chaos Soliton Fract. , 782-794 (2007). J.-M. Flesselles, V. Croquette, and S. Jucquois, “
Period doubling of a torusin a chain of oscillators ,” Phys. Rev. Lett.
72 (18) , 2871-2874 (1994). S. Mondal, S. A. Pawar, and R. I. Sujith, “
Synchronous behaviour of twointeracting oscillatory systems undergoing quasiperiodic route to chaos ,”Chaos , 103119 (2017). Jong Cheol Shin, “
Experimental observation of a torus-doubling transitionto chaos near the ferroelectric phase transition of a KH PO crystal ,” Phys.Rev. E
60 (5) , 5394-5401 (1999). H. Kaplan and J. A. Yorke, “
Lecture Note in Mathematics ,” Springer-Verlag, Berlin , 228-237 (1979). R. Kuske and T. Erneux, “
Localized synchronization of two coupled solidstate lasers ,” Optics Communications , 125-131 (1997). N. M. Awal, D. Bullara, and I. R. Epstein, “
The smallest chimera: Periodic-ity and chaos in a pair of coupled chemical oscillators ,” Chaos , 013131(2019). G. L. Baker, J. A. Blackburn, and H. J. T. Smith, “
Intermittent Synchro-nization in a Pair of Coupled Chaotic Pendula ,” Phys. Rev. Lett.
81 (3) ,554-557 (1998). Liang Zhao, Ying-Cheng Lai, and Chih-Wen Shih, “
Transition to intermit-tent chaotic synchronization ,” Phys. Rev. E
72 (3) , 036212 (2005). A. Banerjee and D. Sikder, “
Transient chaos generates small chimeras ,”Phys. Rev. E , 032220 (2018). A. V. Oppenheim and A. S. Willsky, with S. H. Nawab, “