- Email: [email protected]
Homoclinic chaos in coupled SQUIDs
Chaos, Solitons and Fractals 99 (2017) 133–140
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Homoclinic chaos in coupled SQUIDs M. Agaoglou a,∗, V.M. Rothos b, H. Susanto c a
Department of Mathematics and Statistics, University of Massachusetts, Lederle Graduate Research Tower, Amherst, MA 01003, USA Lab of Nonlinear Mathematics and Department of Mechanical Engineering, Faculty of Engineering, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece c Department of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, United Kingdom b
a r t i c l e
i n f o
Article history: Received 19 October 2016 Revised 23 March 2017 Accepted 3 April 2017
Keywords: SQUIDs Homoclinic chaos Melnikov theory
a b s t r a c t An rf superconducting quantum interference device (SQUID) consists of a superconducting ring interrupted by a Josephson junction (JJ). The induced supercurrents around the ring are determined by the JJ through the celebrated Josephson relations. We study the dynamics of a pair of parametrically-driven coupled SQUIDs lying on the same plane with their axes in parallel. The drive is through the alternating critical current of the JJs. This system exhibits rich nonlinear behavior, including chaotic effects. We take advantage of the weak damping that characterizes these systems to perform a multiple-scales analysis and obtain amplitude equations, describing the slow dynamics of the system. This picture allows us to expose the existence of homoclinic orbits in the dynamics of the integrable part of the slow equations of motion. Using high-dimensional Melnikov theory, we are able to obtain explicit parameter values for which these orbits persist in the full system, consisting of both Hamiltonian and non-Hamiltonian perturbations, to form so called Shilnikov orbits, indicating a loss of integrability and the existence of chaos. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction A superconducting quantum interference device (SQUID) is a device consisting of tiny loops of superconductors employing the so-called Josephson junction (i.e. a system made up of two superconductors, separated by a thin insulating layer [1]). SQUIDs can measure extremely weak signals due to their characteristic related to the flux quantization in a superconducting loop [2]. There are two main kinds of SQUIDs, the dc and the rf SQUIDs. The main difference of them is that the first has two Josephson junctions whereas the rf SQUID has only one [3]. When driven by an alternating magnetic field, the induced supercurrents around the ring are determined by the Josephson junctions through the celebrated Josephson relations. This system exhibits rich nonlinear behavior, including chaotic effects [4]. Quantum rf SQUIDs have attracted great attention, among others since they constitute essential elements for quantum computing [5]. In this direction, rf SQUIDs have been constructed in, e.g., [6]. Recently, rf SQUIDs have been shown to possibly serve as constituting elements for the so-called nonlinear magnetic metamaterials [7,8]. Metamaterials are artificial, composite, inherently non-magnetic media with (positive or negative) magnetic response. Metama∗
Corresponding author. E-mail addresses: [email protected], [email protected] (M. Agaoglou), [email protected] (V.M. Rothos), [email protected] (H. Susanto). http://dx.doi.org/10.1016/j.chaos.2017.04.003 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
terials provide access to all quadrants of the real permittivitypermeability plane, exhibiting negative refraction index, optical magnetism, and other fascinating properties. The key element for the construction of metamaterials has customarily been the split-ring resonator, which is a subwavelength resonant “particle” with operating frequencies up to the optical range [9]. Thus, the split-ring resonator is effectively a kind of an artificial “magnetic atom” [10]. A periodic arrangement of split-ring resonators in space forms a magnetic metamaterial that exhibits high frequency magnetism and negative permeability [11]. Proposals of using rf SQUIDs as superconducting split-ring resonators promise reduction of losses, which constrain the evanescent wave amplification in these materials, and intrinsic nonlinearities due to the extreme sensitivity of the superconducting state to externally applied fields [12]. Motivated by superconducting metamaterials that consist of arrays of rf SQUIDs, here we consider two coupled inductively by the mutual inductance M, rf SQUIDs as sketched in Fig. 1(a). An rf SQUID itself is normally modeled by an electrical circuit containing an inductance L in series with an ideal Josephson element Ic shunted by a capacitor C and a resistor R. The analogue for the coupled SQUIDs is shown in Fig. 1(b). The parameter L represents the self-inductance of the ring with f being the flux through the inductor, the resistor R plays the role of the resistance of the ring the Josephson junction and the capacitor C models the capacitance across the Josephson junction. The dynamic equation for the flux
134
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
Fig. 1. (a) A sketch of two rf SQUIDs (i.e. each SQUID has only one Josephson junction) that are coupled by a mutual inductance M. The coupling appears as the magnetic field generated by one SQUID will cross the neighboring SQUID. (b) Analogous electric circuit for the coupled SQUIDs in panel (a). Each SQUID is modeled by a circuit containing an inductor L, a resistor R, a capacitor C and a Josephson junction with critical current Ic .
ical resonators of the kind that is typically encountered in applications involving microelectromechanical systems (MEMS) and nanoelectromechanical systems (NEMS) has been considered in [19]. Homoclinic chaos in the same system as in [19], but with nonlinear coupling between the oscillators was considered in [20]. Compared to the previous works, our system has a so-called softening nonlinearity as opposed to stiffening one, i.e. the nonlinearity has the opposite sign. Chaos in superconducting Josephson Junction (JJ) has been studied by many researchers. Chaotic JJs can be used for many purposes. Is has been shown [21] that the high bandwidth and rapid decorrelation of the laser pulse train could provide an unambiguous autocorrelation signal for measuring round-trip delay time. Compared with commercial Ladar systems which use the time-of-flight measurement, chaotic system use enables better accuracy. Moreover chaotic JJs can be used for the short-distance secure communications. In a chaos-based secure communications, a message is masked in the broadband chaotic output of the transmitter and synchronization between the transmitter and receiver systems is used to recover a transmitted message [22]. The paper is organized as follows. In Section 2, we perform a multi-scale analysis for the coupled system and obtain amplitude equations to describe the slow dynamics of the system of the two SQUIDs. It is in the slow dynamics that the homoclinic orbits are found. In Section 3, we perform the analytical method for the existence of homoclinic orbits in the perturbed system. In particular, we calculate the energy-difference function using the Melnikov integral evaluated on the homoclinic solutions and applying the singular perturbation theory we study the dynamics near resonances. In Section 4, we summarize our analytical and numerical results for the SQUID model. 2. Normal mode amplitude equations Here, we consider small parameters γ = γˆ , βˆ = β + f p cos (ω pt ), and drive amplitude f p = hˆ , where 1. The normal
threading a SQUID ring can be obtained by applying the Kirchhoff laws [7,8], i.e.
mode frequencies for the linear system are
f¨ + γ f˙ + f + βˆ sin(2π f ) = 0
ω12 = 1 + 2π β − λ, ω22 = 1 + 2π β + λ.
where γ ≡
(1.1)
Ic L (damping constant), βˆ ≡ (SQUID parameter), 0 where 0 is the flux quantum, and the derivatives are with √ respect to normalized time t ≡ ω0 τ , with ω0 1/ LC and τ the physical time variable. Now for two SQUIDs arranged in series, let us denote I1 the supercurrent induced in the first SQUID and I2 in the second. Then the differential equations for the flux of the first SQUID f1 and of the second SQUID f2 in the normalized form, are given by [8] 1 R
L C
f¨1 + γ f˙1 + f1 + βˆ sin(2π f1 ) − λ f2 = 0 f¨2 + γ f˙2 + f2 + βˆ sin(2π f2 ) − λ f1 = 0.
(1.2)
Here, λ is the electrostatic coupling coefficient between the two SQUIDs. In terms of the inductance M depicted in Fig. 1 (a), λ = M/L. We wish to study the homoclinic chaos near resonances in the coupled SQUIDs described by (1.2) when the nonlinearity coefficient βˆ , i.e. the critical current of the Josephson junctions, is varying periodically in time with a frequency close to the natural frequency of the system. Such a parametric drive may appear due to a periodic change of the temperature [13–16]. In particular, we show that the unperturbed system has a homoclinic orbit in their collective dynamics by applying the geometrical method of singular perturbation theory near the resonances and Melnikov theory of near integrable Hamiltonian system to predict the chaotic behavior near resonances [17]. Moreover, we state the theorem for the existence of multi-homoclinic orbits near resonance following [17,18]. Previously homoclinic chaos in coupled nonlinear mechan-
(2.1)
We use multiple time scales [23,24] to express the phases
f1,2 (t ) =
√ 3 (A1 (t˜)eiω1 t ± A2 (t˜)eiω2 t + c.c ) + 3/2 f1(,12) (t ) + . . . 2 (2.2)
where f1 is taken with the positive sign and f2 with the negative sign, slow time t˜ = t and the normal mode frequencies are given by (2.1). The drive frequency ωp is related with ω2 as
ω p = 2 ω 2 + , ω 1 = ω 2 + 2 1 . Substituting the multiple time scale expression to the system of f1 , f2 generates equations containing secular terms. By requiring the secular terms to vanish and expressing the complex amplitudes using real amplitudes and phases as
A1 (t˜) = a1 (t˜)ei(x1 (t˜)+( /2−2 1 )t˜) ,
A2 (t˜) = a2 (t˜)ei(x2 (t˜)+ t˜/2) , (2.3)
we obtain the following coupled equations
1 da1 3π 3 β π hˆ = a1 a22 sin(2(x2 − x1 )) + sin(2x1 )a1 − a1 γˆ dt 2ω1 2ω1 2 3π 3 β 2 dx1
= − + 2 1 − (a1 + 2a22 + a22 cos(2(x2 − x1 ))) dt 2 2ω1 π hˆ + cos(2x1 ) 2ω1 1 da2 3π 3 β π hˆ = a2 a21 sin(2(x1 − x2 )) + sin(2x2 )a2 − a2 γˆ dt 2ω1 2ω2 2
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
dx2
3π 3 β 2 =− − (a2 + 2a21 + a21 cos(2(x1 − x2 ))) dt 2 2ω2 π hˆ + cos(2x2 ), 2ω2
dφ = dt (2.4)
where for simplicity we have omitted the tilde in the time variable t˜. In our analysis herein we will again treat the drive hˆ and the damping γˆ in (2.4) as perturbations.
y2 ∂ H0 (x, y, I )
x2 = −δ I − − (3 − δ ) − (1 − δ ) ∂I 4 2 2
which has the Hamiltonian function
δ
x4 H0 (x, y, I, φ ) = − I2 − I + 3− 2 4 4
+ x2 y2 1 −
2.1. Unperturbed equations: setting damping and drive to zero
−
For γˆ = hˆ = 0, we obtain the unperturbed system
da1 dt dx1 dt da2 dt dx2 dt
2
− 2 1 −
2
−
ω1 2 a − a21 (2 + cos(2(x1 − x2 ))), ω2 2
under the scaling
a1 −→ a1
2ω2 , 3π 3 β
a2 −→ a2
(2.5)
2ω1 . 3π 3 β
Introducing the new variables (see [25])
B=
1 2 a , 2 1
1 2
θ = x1 − x2 , I = (a21 + a22 ), φ = x2 , δ =
ω1 , t → t/2 ω2
the unperturbed system above becomes
dB dt dθ dt dI dt dφ dt
= −2B(I − B ) sin(2θ )
= 1 − I (2 − δ + cos(2θ )) + B 4 −
δ2 + 1 + 2 cos(2θ ) δ
=0 = −δ I −
4
− B(2 − δ + cos(2θ ))
(2.6)
which can be written in a Hamiltonian form as
B˙ = −∂θ H˜ 0 ,
θ˙ = ∂B H˜ 0 , I˙ = −∂φ H˜ 0 , φ˙ = ∂I H˜ 0
with the Hamiltonian
(2.7)
2
2
dI ∂ H0 (x, y, I ) =− =0 dt ∂φ
is a compact normally hyperbolic invariant manifold with boundary that admits two homoclinic manifolds W0± . Solutions on are typically periodic, as one can readily verify by substituting x = y = 0 into the φ -components of (2.9). An exception occurs for the action value
I0 = −
4δ
,
which is labeled as C, a circle of fixed points, on the manifold . The fixed points in this circle represent a 1-parameter family of nonlinear normal modes. The 2D invariant manifold , has a 3D stable and unstable manifolds Ws (), Wu (), which coincide to form a 3D homoclinic manifold = W s () ∩ W u (). Trajectories on the homoclinic manifold are homoclinic orbits that connect the origin of the (x, y) plane to itself. Moreover, the constant value of the Hamiltonian along such an orbit is equal to its value at the origin, namely H0 (0, 0, I ) = H0 (x, y, I ) = H˜0 (B, θ , I ) which yields an equation for the homoclinic orbits in terms of the action-angle variables, i.e.
Bh (t, I ) =
2δ a2 q cosh(2at ) + p
θ (t, I ) = tan (2.8)
δ2 + 1 2 δ2 + 1 −x y 2− 2δ 2δ
dy ∂ H0 (x, y, I ) δ + 1 δ2 + 1 = = x3 1 − + xy2 2 − dt ∂x 2δ 2δ − x ( I ( 3 − δ ) − 1 )
(2.10)
1 = (x, y, I, φ ) : x = y = 0, 1 < δ < 3, I > 3−δ
h
At this point we have to emphasize that I and H˜0 are constants of the motion in the unperturbed system. Note that (B, θ , I, φ ) ∈ R+ × S × R+ × S, where S is the unit circle, and R+ are the non-negative reals. √ √ Now, introducing the Cartesian variables x = 2Bcosθ , y = 2Bsinθ , the system (2.7) becomes
dx ∂ H0 (x, y, I ) =− = −y3 1 − dt ∂y + y ( I ( 1 − δ ) − 1 )
x2 y2 ( I ( 3 − δ ) − 1 ) − ( I ( 1 − δ ) − 1 ) 2 2
δ2 + 1 H˜ 0 (B, θ , I ) = − I − I + B 2 − 2 4 δ − B(I (2 − δ ) − 1 ) − B(I − B )cos(2θ ) δ
δ 2 + 1 y4 δ2 + 1 + 1− 2δ 4 2δ
The point x = y = 0 is a fixed point of our unperturbed system, as given by Eqs. (2.9a) and (2.9b), and it is a saddle point for 1 < δ < 3 and I > 1 /3 − δ . This result comes from the fact that this fixed point is saddle for values of the positive constant of motion I, that satisfy the inequality [I (1 − δ ) − 1 ][I (3 − δ ) − 1 ] < 0. The set
ω2 2 a − a22 (2 + cos(2(x2 − x1 ))) ω1 1
= a21 a2 sin(2(x1 − x2 )) =−
δ2 + 1 4δ
(2.9d)
2.2. Analytical expressions for homoclinic orbits
= −a1 a22 sin(2(x1 − x2 )) =−
135
(2.9a)
−1
I ( δ − 3 ) + 1 tanh(at ) I ( 1 − δ ) − 1
−a(δ 2 − 1 ) χ1h (t, I ) =
tan−1 − +
p2 − q2
δI −
(2.11)
4
p−q tanh(at ) p+q
t + x1 ( 0 )
φ h (t, I ) = χ1h (t, I ) − θ h (t, I )
(2.12)
(2.13) (2.14)
where
q ≡ I ( δ 2 − 1 ) + 2 δ 1 > 0 p = − 1 (δ 2 − 4δ + 1 ) − I (δ 3 − 6δ 2 + 7δ − 2 )
2
a2 = 21 − 2I 1 (δ − 2 ) − I2 (δ − 3 )(δ − 1 ). (2.9b) (2.9c)
In Fig. 2(a,b), we sketch the scenario of the dynamics explained here. In Fig. 3 we plot the formation of the homoclinic orbit to . The homoclinic manifolds that we obtained in terms of the original amplitude equations (2.4) with hˆ = γˆ = 0 are shown in Fig. 4.
136
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
We assume that every resonant curve C is connected to some other such state for γˆ = 0, hˆ = 0, i.e. there are “single-jump” fast transitions between resonances [17,25]. More precisely, we assume that the stable and unstable directions associated with C give rise to a manifold W0 (C ) of orbit homoclinic to C, see Fig. 2. We shall consider homoclinic manifolds with two components W0+ (C ) and W0− (C ). The orbits in both W0+ (C ), W0− (C ) connect the same pair of points with the phase shift vector
x˜ = lim (xh (t ) − xh (−t )) t→∞
= (x, y, I, φ ) = (0, 0, 0, lim (φ h (t ) − φ h (−t ))) t→∞
being the same for any solution xh (t) in W0+ (C ) ∪ W0− (C ), where φ = x1 − θ . The orbits given by Eqs. (2.11)–(2.14) are homoclinic to (see Fig. 3). In Fig. 2(b) we sketch some of them. At resonance, for I = I0 , the orbits are heteroclinic, connecting fixed points that are φ apart, where
φ = x 1 − θ Fig. 2. (a) The homoclinic manifold of C. (b) A sketch of the assumption that the orbits in both W0+ (C ) and W0− (C ) connect the same pairs of points, i.e. the phase shift vector. (c) A 3− chain of unperturbed orbits. See Definition 2.1 for the meaning of xj (t), j = 1, 2, . . . . The pictures are reprinted from [17].
and
2a ( δ 2 − 1 ) −1 x 1 =
tanh
p2 − q2
θ = 2 tan
−1
p−q p+q
I0 (δ − 3 ) + 1 I0 (1 − δ ) − 1
Therefore
p − q 2a ( δ 2 − 1 ) I0 (δ − 3 ) + 1 −1 φ =
tanh − 2 tan−1 p+q I0 (1 − δ ) − 1 p2 − q2
B
10
5
0 -4
28.9 28.95 -2
29 29.05
0 2
29.1 29.15
4
I
Fig. 3. Orbits homoclinic to . For I = I0 (the orbit in the middle), dφ /dt = 0 on , and the orbit is heteroclinic, connecting fixed points on that are φ apart. For IࣤI0 , dφ /dtࣤ0 on . The parameters are δ = 1.49704, = −110, 1 = 27.6085.
Definition 2.1. Let us consider for our unperturbed system (2.9) a set of solutions {x j (t )}Nj=1 each homoclinic to the torus (circle) C with the property
I=50
lim x1 (t ) = b0 ,
t→−∞
10
a2
lim x j−1 (t ) = lim x j (t ),
t→∞
t→−∞
j = 2, . . . , N
and let {Pj }Nj=1 be a sequence of +1’s and −1’s. The solutions xj (t) are said to form an N− chain with base point b0 and jump sequence {Pj }Nj=1 if ∀t ∈ R we have
9 8 7
x (t ) ∈ j
I=20 0 a1 sin(x1-x 2)
5
0 a cos(x -x ) 1
1
5
I=30
6 10
Remark 1. In resonance problems the fixed points in C typically represent same solution with different initial conditions in a moving frame. In such a case the phase shift (φ ) on an orbit homoclinic to the underlying solution will be independent of the initial phase. One can built “chains” out of the solutions xh (t ) = (x, y, I, θ ) to connect fixed points that are Nφ apart for some integer N (see Fig. 2(b,c)). The members of such a chain can be taken from both components of the homoclinic manifold, and they connect N + 1 different resonant states altogether. We shall seek actual solutions staying close to these chains, i.e., N− chains of unperturbed orbits will serve as “shadowing sets” for real orbits exhibiting N excursions from the set of resonances. We first define an N− chain more precisely.
-5
-10
-5
W0+ (C ) W0− (C )
i f Pj = +1 i f Pj = −1
We denote an N − chain with base point b0 by XN (b0 ). 3. Homoclinic intersections in the damped, driven system
2
Fig. 4. Results of a numerical integration of Equations (2.4) with hˆ = γˆ = 0, ω1 = 1.58114, ω2 = 1.22474, 1 = 35 and = −150.
In this section we want to study how the drive and the damping as perturbations affect the dynamics. Writing the drive
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
hˆ = 2εω1 h¯ /π and the damping γˆ = 2ε γ¯ , ε 1, the perturbed system (2.4) in terms of the action-angle variables takes the form
dB = −2B(I − B ) sin(2θ ) + ε Bh sin(2(φ + θ )) − ε Bγ dt dθ δ2 + 1 = 1 − I (2 − δ + cos(2θ )) + B 4 − + 2 cos(2θ ) dt δ εh + (cos(2(φ + θ ))) − δ cos(2φ ) 2 dI = ε h(B sin(2(φ + θ )) + (I − B )δ sin(2φ )) − εγ I dt dφ
h = −δ I − − B(2 − δ + cos(2θ )) + ε δ cos(2φ ), (3.1) dt 4 2 where we have dropped the bar in h¯ and γ¯ for simplicity. This system can also be written as
dB ∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) =− + ε kB = − dt ∂θ ∂θ ∂ H˜ (B, θ , I, φ ) 1 B +ε − +g
∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) ∂ H˜ 1 (B, θ , I, φ ) + ε kθ = +ε ∂B ∂B ∂B ∂ H˜ (B, θ , I, φ ) ε − 1 + gI ∂φ dφ ∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) ∂ H˜ 1 (B, θ , I, φ ) = + ε kφ = +ε (3.2) dt ∂I ∂I ∂I where
h (B cos(2(φ + θ )) + δ (I − B ) cos(2φ )) 2
and the dissipative perturbations are given by gB = −Bγ and gI = −γ I. In Cartesian variables, the perturbed system takes the form
dx ∂ H0 (x, y, I ) ∂ H1 (x, y, I, φ ) =− +ε − + gx dt ∂y ∂y ∂ H (x, y, I, φ ) dy ∂ H0 (x, y, I ) 1 = +ε + gy dt ∂x ∂x ∂ H (x, y, I, φ ) dI 1 =ε − + gI dt ∂φ
dφ = dt
H1 (x, y, I, φ ) =
h x2 + y2 2
2
+δ I − gx = −
γx 2
,
gy = −
γy 2
x x2 + y2
x +y 2
2
2 ,
cos φ −
cos(2φ )
gI = −γ ,
dI = ε hIδ sin(2φ ) − εγ I, dt
(3.4)
dφ
ε hδ = −δ I − + cos(2φ ). dt 4 2
(3.5)
According [25] and [17] in order to study the slow dynamics, which is induced by the perturbation on Mε near resonance, we √ set a slow variable I = I0 + ε η into Eqs. (3.4) and (3.5) along √ with a slow time scale τ = εt to get
√ dη = (hδ sin(2φ ) − γ0 )(I0 + ε η ) dτ √ dφ εhδ = −δη + cos(2φ ) dτ 2
(3.6) (3.7)
(3.8)
dφ = −δη dτ
(3.9)
which can be obtained from the Hamiltonian
H (η , φ ) = −
y x2 + y2
sin φ
gφ = 0
For 0 < ε 1 the perturbed system (3.3):
p1c (ξ ) = 0,
of Cr -class and ⊂ Mε , while for ε = 0, we have M0 = . • Mε has codimension-one local stable and unstable manifolds s (M ) and W u (M ). Wloc ε ε loc Solutions of (3.3) admit two different time scales near : fast in x, y-variables and slow motion in (I, φ ).
2
+ I0 γ φ +
1 I0 hδ cos(2φ ) 2
(3.10)
1 −1 γ sin 2 hδ 1 2
π − sin−1
(3.11)
γ
(3.12)
hδ
Our goal is to prove the existence of generalized Shilnikov orbits homoclinic or heteroclinic to the slow sinks [17]. We first compute the energy-function N H associated with (3.10) [17]
N H = H ( φ0 + N φ ) − H ( φ0 ) −
N k=1
+∞
−∞
→ DH0 , − g dt
(3.13)
− → where g = (−H1y + gx , H1x + gy , −H1φ + gI , H1I + gφ ) is the perturbed vector field. In order to find the energy-function N H (3.13) we calculate the following
dH˜ 1 ∂ H0 x ∂ H0 y + g + g + dt ∂x ∂y dH˜ 1 dφ dθ =− −γI −γB dt dt dt
→ DH0 , − g=−
• near C there exists a unique, locally invariant manifold
Mε = {(x, y, I, φ ) : x = xε (y, I, φ )}
δη2
The phase portrait of Eqs. (3.6) and (3.7) in the new variables (η, φ ) for different values of ε is given in Fig. 5. The slow local Hamiltonian H admits two saddle-type fixed points p1s , p2s as well as two centers p1c , p2c given by
p1c (ξ ) = 0,
(3.3)
dη = (hδ sin(2φ ) − γ0 )I0 dτ
and
∂ H (x, y, I, φ ) ∂ H0 (x, y, I ) 1 +ε + gφ ∂I ∂I
where
By setting B = 0 in the two latests Equations of (3.1) we obtain the equations that describe the dynamics on Mε near the resonance I = I0 , i.e.
ε we obtain
dθ = dt dI = dt
3.1. Dynamics near resonance
From the leading terms in Eqs. (3.6) and (3.7) independent of
∂θ
H˜ 1 (B, θ , I, φ ) =
137
∂ H0 I g ∂I
where x1 = θ + φ , B(±∞ ) = 0 and dθ /dt = −I (2 − δ − cos(2θ )) −
1 and dx1 /dt = −Iδ − /4 − B(1 − δ /2δ ). Therefore, we obtain
−
N l=1
=−
+∞
−∞ N l=1
− → DH0 , g dt +∞
−∞
dH˜ 1 −
dt
−γI
dφ dθ −γB dt dt
dt
138
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
−
2ξ N hδ
− 2 φ − θ −
L˜ I0
(3.16)
If the dissipation parameter satisfies the condition
ξ<
( N φ / 2 ) sin(Nφ ) 1 + 2 sin sin(φ /2 ) 2N hδ
− 2 φ − θ −
L˜ I0
,
then the energy function admits transverse zeros, where
L˜ =
−2I0 (δ 2 − 1 ) + 2δ 1 tan−1 √ ( δ − 1 ) δ 2 − 6δ + 1
δ−1 δ 2 − 6δ + 1
√
I0 (δ − 3 ) + 1 I0 (1 − δ ) − 1
3.2. A homoclinic connection to the sink pε We are finally in a position to show the existence of an orbit homoclinic to the sink pε . To achieve this, we first show that there exists a homoclinic orbit that approaches pε asymptotically backward in time, and approaches the perturbed annulus Mε asymptotically forward in time. We then estimate the conditions under which the perturbed counterpart of the point, which is reached on forward in time in the unperturbed system, lies within the basin of attraction of the sink pε on Mε . This gives us an estimate for the possibility of obtaining a Shilnikov orbit that connects the sink back to itself. This first step is done by finding the conditions for which the energy function has simple zeros. Theorem 1. System (3.3) depends on a vector κ of system parameters. If (i) N H ( p1c (κ ), κ ) = 0 (ii) D pc N H ( p1c (κ ), κ ) = 0, Dξ N H ( p1c (κ ), κ ) = 0 (iii) The point (0, φc (ζ ) + Nφ ) lies in the domain of attraction of the sink. Then there exist a codimension-one set M+ such that: (a) Our system admits an N−pulse Shilnikov-type orbit x , homoclinic to the point p . The takeoff point for this orbit is p , and the jump sequence of the orbits is given by Fig. 5. (a) Numerical phase portraits of Equations (3.6) and (3.7) with ε = 0. (b) Numerical phase portrait of Equations (3.6) and (3.7) with ε = 1. The parameters are δ = 1.55, = −150, h = 1, ζ = 0.45, γ = hζ δ .
=
hδ I0 2
sin
N φ 2
sin 2φ
2 sin(2φ ) sin(Nφ ) + γ I0 Nφ (3.14)
where I0 = − /4δ is the resonance, N is the number of homoclinic jumping, φ = x1 − θ and
2I0 (δ 2 − 1 ) + 2δ 1 L˜= − tan−1 √ ( δ − 1 ) δ 2 − 6δ + 1
δ−1 √ 2 δ − 6δ + 1
I0 (δ − 3 ) + 1 . I0 (1 − δ − 1 )
Next, consider again and calculate the term H (φ0 + Nφ ) − H (φ0 ) on the slow manifold, i.e.
Remark 2. According to the Theorem 1 we find solutions of N H ( p1c , κ ) = 0 where p1c center point and κ the parameters of system. Some elementary algebra turns this condition (N H ( p1c , κ ) = 0) into
1 − ζ 2 (1 − cos(2N φ ) ) = −ζ sin(2N φ )
N H
H ( φ0 + N φ ) − H ( φ 0 ) =
hI0 δ (cos(2(φ0 + Nφ )) 2 − cos(2φ )) + γ I0 Nφ .
From (3.14) and (3.15), we obtain
N H (φ ) = cos(2(φ0 + Nφ )) − cos(2φ0 ) sin + sin
N φ 2
φ 2
2 sin(2φ0 ) sin(Nφ )
(3.15)
j = 1, ..., N − 1.
(b) There also exists another codimension-one set M− that yields similar homoclinic orbits with jump sequence {−Pj }Nj=1 . (c) Both surfaces M+ and M− have open neighborhoods in the (κ , ) parameter space for which our system admits Smale horseshoes [17].
+ I0 N θ γ + N L˜γ
Pj+1 = σ sign( j H ( p1c , κ )Pj ,
P1 = +1,
sin +2
N φ 2
sin 2φ
1 sin sin−1 (ζ ) sin(Nφ ) 2
− 2 − 2 φ − θ −
1 L˜ Nζ I0
(3.17)
where ζ = ξδ .The Eq. (3.17) gives the perturbation parameter value
ζ=
1 − cos(2Nφ )
( N φ / 2 ) (1 − cos(2Nφ ))2 + −2BˆN + 2 sin sin(Nφ ) − sin(2Nφ ) sin(φ /2 )
2
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
139
(3.18) where
Bˆ = −2φ − θ −
1 L˜. I0
3
Of course, this result is meaningful only for nonzero values
k ∈ Z.
2
(3.19)
The second transversality condition in Theorem 1 gives
D ξ [ H ( N
p1c ,
κ )] = 0
B
φ
kπ = , N
1 0
(3.20)
20
1.5
to be satisfied whenever (3.18) and (3.19) hold. Carrying out the differentiation, we obtain that (3.20) is violated if
1 − ζ 2 −sin(2N φ ) + 2
N φ
sin
1
−1
sin sin 2 sin 2φ 2
(ζ ) sin(Nφ )
1 − 2 − 2φ − θ − L˜ N I0 = ζ (1 − cos(2Nφ ))
(3.21)
However, under condition (3.19), Eqs. (3.17) and (3.21) cannot be satisfied simultaneously. Therefore, the second condition of (ii) of Theorem 1 is satisfied. The first condition in (ii) is satisfied by inequality (2.6). It remains to verify condition (iii) of Theorem 1 to ensure that the landing point of N-pulse orbit taking off from a slow sink lies in the domain of attraction of one of the sinks. According to [17], the landing point lies near the angle φc (ζ ) + Nφ , where φ c (ζ ) is the φ coordinate of the center p1c . For small enough ε , the landing point will fall in the basin of a slow sink if the point
( 0 , φc ( ζ ) + N φ ) falls in D01 . This can be verified by comparing the energy of the approximate landing point to the “energy” H (η, φ ) of the boundary of D01 ∪ D02 . Note, however, that this energy is not well-defined, since H (η, φ ) is only a local Hamiltonian. To obtain a uniquely defined energy, we consider the point that lies in the interval [ π2 , 32π ] and is kπ apart from the approximation landing point. The φ -coordinate of this point can be computed as
φmN (ζ ) = φs (ζ ) + φc (ζ ) + Nφ − φs (ζ ) modπ ,
(3.22)
where φ c (ζ ), φ s (ζ ) are φ -components of center and saddle points. N (ζ ) > φ (ζ ), then redefine φ N (ζ ) by subtracting π from it. If φm s m Now, if N H ( 0 , φm (ζ )) > H (0, φs (ζ ))
cos(2φs (ζ )) − cos(2φ (ζ )) > 2ζ (φs (ζ ) − φ (ζ )) N m
(3.24)
:= ξδ , ξ =
ζ γ /h then (0, φmN (ζ )) ∈ D01 , which implies that (0, φs (ζ ) + Nφ ) ∈ D01 ∪ D02 . Now from the above we want to write an approximate condition which guarantees that this orbit approaches pε as t → −∞. To succeed in this, we find a condition for which the unperturbed heteroclinic orbit, which asymptotes to pc as t → +∞ returns back to the point on the circle of fixed points that is inside the homoclinic separatrix loop connecting the saddle qc to itself. This kind of an orbit is shown in Fig. 6. This condition is presented in terms of the difference φ between the asymptotic values of the angular variable φ as
φs < φc + φ < φ m
0.5 0 30 -0.5
Fig. 6. The heteroclinic orbit given by Equations (2.11) and (2.14) with I = I0 , superimposed with the phase portrait of the unperturbed scaled system on ε near resonance, given by Equations (3.8) and (3.9). The parameters are δ = 1.55, = −150, 1 = 35, ζ = 0.45, h = 1 and γ = hζ δ . This value of ζ sets φc = 0.233383 and we can see from the figure that φs < φc + φ < φm .
Eqs. (3.18) and (3.25) define the conditions which are responsible for the existence of orbits homoclinic to the sink pε . What we want now is to relate these results to the actual physical parameters of the coupled resonators. As we have already defined δ , the value of λ is then given by λ = (−δ 2 − 2π βδ 2 + 1 + 2π β )
/(δ 2 + 1 ). The scaled frequency
1 is 1 = (ω1 − ω2 )/ = ( 1 + 2π β − λ − 1 + 2π β + λ)/2 , so after fixing it is also determined by δ . The ratio ζ = γ /(hδ ), between the damping coefficient and the drive amplitude, has to be positive in order for the damping coefficient γ to be positive and has the standard physical meaning of energy dissipation. The parameter ζ is positive if the inequality,
−2BˆN − sin(2Nφ ) + 2
sin(Nφ /2 ) sin(Nφ ) > 0 sin(φ /2 )
(3.26)
is satisfied. The φ values of the fixed points are shown in Fig. 7(a) along with φc + φ and φ m for a particular value of . The parameter values for which these φ values satisfy the condition (3.25) are displayed in Fig. 7(b), which outlines the values of the coupling and driving frequency, for which orbits homoclinic to the sink pε exist. 4. Conclusion
(3.23)
or N m
25
J
1
(3.25)
where φ m is the maximal value of φ on the homoclinic orbit, connecting the saddle qc to itself.
In this work we studied chaotic dynamics of a pair of parametrically-driven coupled SQUIDs arranged in series and we found under which conditions it can exist by using high dimensional Melnikov theory. This was achieved by applying a method of Kovacic and Wiggins on transformed amplitude equations that were derived from the equations of motion, which model an actual experimental realization of coupled SQUIDs in metamaterials. We considered the amplitude of the drive and the damping to be small and obtained explicit expressions for orbits homoclinic to a twodimensional invariant annulus in the unperturbed equations. At resonance, we were able to calculate the energy function analytically, and provide a primary condition for having homoclinic orbits in the full, perturbed equations. By further studying the effects of perturbations on the invariant annulus near resonance, we found a secondary condition for the existence of orbits homoclinic to a fixed point of a saddle-focus type. We used a numerical scheme to verify our theoretical predictions. Such Shilnikov homoclinic orbits give rise to a particular type of horseshoe chaos, which can be expected in the dynamics of the full system for parameter values
140
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
Fig. 7. (a) The values of φs , φc , φc + φ and φ m as functions of δ , for = −150. For δ > 1.524 the condition is satisfied. (b) Parameter values for which the condition is satisfied are indicated in gray. In both figures = 0.01.
in the vicinity of those presented here. Our work provides parameter regions where homoclinic chaos can be observed that may be exploited for, e.g., chaos communication. As future work we can consider the possibility to have two SQUIDs with one above the other in which case the coupling is stronger and has opposite sign. Acknowledgments M.A. thanked Professor Ron Lifshitz for fruitful discussions. The work of M.A. has been co-financed from resources of the operational program “Education and Lifelong Learning” of the European Social Fund and the National Strategic Reference Framework (NSRF) 2007-2013 within the framework of the Action State Scholarships Foundation’s (IKY) mobility grants programme for the short term training in recognized scientific/research centers abroad for candidate doctoral or postdoctoral researchers in Greek universities or research centers. The work of V.R. has been co-financed by the European Union (European Social Fund ESF) and Greek national funds through the Operational Program “Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF) Research Funding Program D.534 MIS: 379337: THALES Investing in knowledge society through the European Social Fund. V.R and M.A acknowledge support from FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-606096). H.S. is supported by the Engineering and Physical Sciences Research Council (Grants No. EP/M024237/1). References [1] Josephson BD. Possible new effects in superconductive tunnelling. Phys Lett 1962;1:251–3. [2] Clarke J, Braginski AI, editors. The SQUID handbook. Fundamentals and technology of SQUIDS and SQUID systems, vol. I. Weiheim: Wiley-VCH; 2004. [3] Likharev KK. Dynamics of Josephson junctions and circuits. Gordon and Breach, Philadelphia; 1986. [4] Fesser K, Bishop AR, Kumar P. Chaos in rf SQUIDs. Appl Phys Lett 1983;43:123–4. [5] Bocko MF, Herr AM, Feldman MJ, M J. Prospects for quantum coherent computation using superconducting. IEEE Trans Appl Supercond 1997;7:3638–41.
[6] Yamashita T, Takahashi S, Maekawa S. Superconducting pi qubit with three josephson junctions. Appl Phys Lett 2006;88:132501–1-132501-3. [7] Lazarides N, Tsironis GP, G P. rf superconducting quantum interference device metamaterials. Appl Phys Lett 2007;16:163501–163501-3. [8] Lazarides N, Tsironis GP, Eleftheriou N. Dissipative discrete breathers in rf SQUID metamaterials. Nonlinear Phenom. Complex Syst. 2008;11:250–8. [9] Yen TJ, Padilla WJ, Fang N, Vier DC, Smith DR, Pendry JB, et al. Terahertz magnetic response from artificial materials. Science 2004;303:1494–6. [10] Caputo JG, Gabitov I, Maimistov AI. Electrodynamics of a split-ring josephson resonator in a microwave line. Phys Rev B 2012;85:205446. [11] Linden S, Enkrich C, Dolling G, Klein MW, Zhou JF, Koschny T, et al. Photonic metamaterials: magnetism at optical frequencies. IEEE J Selec Top Quant Electron 2006;12:1097–105. [12] Chen HT, Yang H, Singh R, OHara JF, Azad AK, Stuart A, et al. Tuning the resonance in hightemperature superconducting terahertz metamaterials. Phys Rev Lett 2010;105:247402. [13] Overhauser AW. Temperature dependence of the josephson current. Phys Rev Lett 1999;83:180. [14] Han S, Wolf EL. Temperature dependence of the critical josephson current in superconductor insulatornormal-metal proximity junctions near Tc . Phys Rev B 1987;35:4669. [15] Hikino S-i, Mori M, Takahashi S, Maekawa S. Temperature dependence of josephson current in a superconductor/ ferromagnet/superconductor junction. Physica C 2007;460-462:1323–4. [16] Veldhorst M, Snelder M, Hoek M, Gang T, Guduru VK, Wang XL, et al. Josephson supercurrent through a topological insulator surface state. Nat. Mater. 2012;11:417–21. [17] Haller G. Chaos near resonance. Springer, New York; 1999. [18] Menon G, Haller G, Rothos V. Silnikov manifolds in coupled nonlinear schrödinger equations. Phys Lett A 1999;263:175–85. [19] Kenig E, Tsarin YA, Lifshitz R. Homoclinic orbits and chaos in a pair of parametrically driven coupled nonlinear resonators. Phys Rev E 2011;84:016212. [20] Zhang W, Huang YT, Yao MH. Multi-pulse homoclinic orbits and chaotic dynamics of a parametrically excited nonlinear nano-oscillator with coupled cubic nonlinearities. Sci China Phys Mech Astron. 2014;57:1098–110. [21] Dana SK, Sengupta DC, Edoh KD. Chaotic dynamics in josephson junction. IEEE Trans on Circ and Syst I 2001;48:990–6. [22] Myneni K, Barr TA, Reed BR, B R. High-precision ranging using a chaotic laser pulse train. Appl Phys Lett 2001;78:1496–8. [23] Hand LN, Finch JD. Analytical mechanics. Cambridge Univ. Press, Cambridge; 1998. [24] Strogatz SH. Nonlinear dynamics and chaos. Addison-Wesley, Reading MA; 1994. [25] Kovacic G, Wiggins S, S. Orbits homoclinic to resonances, with an application to chaos in a model of the forced and damped sine-gordon equation. Physica D 1992;57:185–225.
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Homoclinic chaos in coupled SQUIDs M. Agaoglou a,∗, V.M. Rothos b, H. Susanto c a
Department of Mathematics and Statistics, University of Massachusetts, Lederle Graduate Research Tower, Amherst, MA 01003, USA Lab of Nonlinear Mathematics and Department of Mechanical Engineering, Faculty of Engineering, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece c Department of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, United Kingdom b
a r t i c l e
i n f o
Article history: Received 19 October 2016 Revised 23 March 2017 Accepted 3 April 2017
Keywords: SQUIDs Homoclinic chaos Melnikov theory
a b s t r a c t An rf superconducting quantum interference device (SQUID) consists of a superconducting ring interrupted by a Josephson junction (JJ). The induced supercurrents around the ring are determined by the JJ through the celebrated Josephson relations. We study the dynamics of a pair of parametrically-driven coupled SQUIDs lying on the same plane with their axes in parallel. The drive is through the alternating critical current of the JJs. This system exhibits rich nonlinear behavior, including chaotic effects. We take advantage of the weak damping that characterizes these systems to perform a multiple-scales analysis and obtain amplitude equations, describing the slow dynamics of the system. This picture allows us to expose the existence of homoclinic orbits in the dynamics of the integrable part of the slow equations of motion. Using high-dimensional Melnikov theory, we are able to obtain explicit parameter values for which these orbits persist in the full system, consisting of both Hamiltonian and non-Hamiltonian perturbations, to form so called Shilnikov orbits, indicating a loss of integrability and the existence of chaos. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction A superconducting quantum interference device (SQUID) is a device consisting of tiny loops of superconductors employing the so-called Josephson junction (i.e. a system made up of two superconductors, separated by a thin insulating layer [1]). SQUIDs can measure extremely weak signals due to their characteristic related to the flux quantization in a superconducting loop [2]. There are two main kinds of SQUIDs, the dc and the rf SQUIDs. The main difference of them is that the first has two Josephson junctions whereas the rf SQUID has only one [3]. When driven by an alternating magnetic field, the induced supercurrents around the ring are determined by the Josephson junctions through the celebrated Josephson relations. This system exhibits rich nonlinear behavior, including chaotic effects [4]. Quantum rf SQUIDs have attracted great attention, among others since they constitute essential elements for quantum computing [5]. In this direction, rf SQUIDs have been constructed in, e.g., [6]. Recently, rf SQUIDs have been shown to possibly serve as constituting elements for the so-called nonlinear magnetic metamaterials [7,8]. Metamaterials are artificial, composite, inherently non-magnetic media with (positive or negative) magnetic response. Metama∗
Corresponding author. E-mail addresses: [email protected], [email protected] (M. Agaoglou), [email protected] (V.M. Rothos), [email protected] (H. Susanto). http://dx.doi.org/10.1016/j.chaos.2017.04.003 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
terials provide access to all quadrants of the real permittivitypermeability plane, exhibiting negative refraction index, optical magnetism, and other fascinating properties. The key element for the construction of metamaterials has customarily been the split-ring resonator, which is a subwavelength resonant “particle” with operating frequencies up to the optical range [9]. Thus, the split-ring resonator is effectively a kind of an artificial “magnetic atom” [10]. A periodic arrangement of split-ring resonators in space forms a magnetic metamaterial that exhibits high frequency magnetism and negative permeability [11]. Proposals of using rf SQUIDs as superconducting split-ring resonators promise reduction of losses, which constrain the evanescent wave amplification in these materials, and intrinsic nonlinearities due to the extreme sensitivity of the superconducting state to externally applied fields [12]. Motivated by superconducting metamaterials that consist of arrays of rf SQUIDs, here we consider two coupled inductively by the mutual inductance M, rf SQUIDs as sketched in Fig. 1(a). An rf SQUID itself is normally modeled by an electrical circuit containing an inductance L in series with an ideal Josephson element Ic shunted by a capacitor C and a resistor R. The analogue for the coupled SQUIDs is shown in Fig. 1(b). The parameter L represents the self-inductance of the ring with f being the flux through the inductor, the resistor R plays the role of the resistance of the ring the Josephson junction and the capacitor C models the capacitance across the Josephson junction. The dynamic equation for the flux
134
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
Fig. 1. (a) A sketch of two rf SQUIDs (i.e. each SQUID has only one Josephson junction) that are coupled by a mutual inductance M. The coupling appears as the magnetic field generated by one SQUID will cross the neighboring SQUID. (b) Analogous electric circuit for the coupled SQUIDs in panel (a). Each SQUID is modeled by a circuit containing an inductor L, a resistor R, a capacitor C and a Josephson junction with critical current Ic .
ical resonators of the kind that is typically encountered in applications involving microelectromechanical systems (MEMS) and nanoelectromechanical systems (NEMS) has been considered in [19]. Homoclinic chaos in the same system as in [19], but with nonlinear coupling between the oscillators was considered in [20]. Compared to the previous works, our system has a so-called softening nonlinearity as opposed to stiffening one, i.e. the nonlinearity has the opposite sign. Chaos in superconducting Josephson Junction (JJ) has been studied by many researchers. Chaotic JJs can be used for many purposes. Is has been shown [21] that the high bandwidth and rapid decorrelation of the laser pulse train could provide an unambiguous autocorrelation signal for measuring round-trip delay time. Compared with commercial Ladar systems which use the time-of-flight measurement, chaotic system use enables better accuracy. Moreover chaotic JJs can be used for the short-distance secure communications. In a chaos-based secure communications, a message is masked in the broadband chaotic output of the transmitter and synchronization between the transmitter and receiver systems is used to recover a transmitted message [22]. The paper is organized as follows. In Section 2, we perform a multi-scale analysis for the coupled system and obtain amplitude equations to describe the slow dynamics of the system of the two SQUIDs. It is in the slow dynamics that the homoclinic orbits are found. In Section 3, we perform the analytical method for the existence of homoclinic orbits in the perturbed system. In particular, we calculate the energy-difference function using the Melnikov integral evaluated on the homoclinic solutions and applying the singular perturbation theory we study the dynamics near resonances. In Section 4, we summarize our analytical and numerical results for the SQUID model. 2. Normal mode amplitude equations Here, we consider small parameters γ = γˆ , βˆ = β + f p cos (ω pt ), and drive amplitude f p = hˆ , where 1. The normal
threading a SQUID ring can be obtained by applying the Kirchhoff laws [7,8], i.e.
mode frequencies for the linear system are
f¨ + γ f˙ + f + βˆ sin(2π f ) = 0
ω12 = 1 + 2π β − λ, ω22 = 1 + 2π β + λ.
where γ ≡
(1.1)
Ic L (damping constant), βˆ ≡ (SQUID parameter), 0 where 0 is the flux quantum, and the derivatives are with √ respect to normalized time t ≡ ω0 τ , with ω0 1/ LC and τ the physical time variable. Now for two SQUIDs arranged in series, let us denote I1 the supercurrent induced in the first SQUID and I2 in the second. Then the differential equations for the flux of the first SQUID f1 and of the second SQUID f2 in the normalized form, are given by [8] 1 R
L C
f¨1 + γ f˙1 + f1 + βˆ sin(2π f1 ) − λ f2 = 0 f¨2 + γ f˙2 + f2 + βˆ sin(2π f2 ) − λ f1 = 0.
(1.2)
Here, λ is the electrostatic coupling coefficient between the two SQUIDs. In terms of the inductance M depicted in Fig. 1 (a), λ = M/L. We wish to study the homoclinic chaos near resonances in the coupled SQUIDs described by (1.2) when the nonlinearity coefficient βˆ , i.e. the critical current of the Josephson junctions, is varying periodically in time with a frequency close to the natural frequency of the system. Such a parametric drive may appear due to a periodic change of the temperature [13–16]. In particular, we show that the unperturbed system has a homoclinic orbit in their collective dynamics by applying the geometrical method of singular perturbation theory near the resonances and Melnikov theory of near integrable Hamiltonian system to predict the chaotic behavior near resonances [17]. Moreover, we state the theorem for the existence of multi-homoclinic orbits near resonance following [17,18]. Previously homoclinic chaos in coupled nonlinear mechan-
(2.1)
We use multiple time scales [23,24] to express the phases
f1,2 (t ) =
√ 3 (A1 (t˜)eiω1 t ± A2 (t˜)eiω2 t + c.c ) + 3/2 f1(,12) (t ) + . . . 2 (2.2)
where f1 is taken with the positive sign and f2 with the negative sign, slow time t˜ = t and the normal mode frequencies are given by (2.1). The drive frequency ωp is related with ω2 as
ω p = 2 ω 2 + , ω 1 = ω 2 + 2 1 . Substituting the multiple time scale expression to the system of f1 , f2 generates equations containing secular terms. By requiring the secular terms to vanish and expressing the complex amplitudes using real amplitudes and phases as
A1 (t˜) = a1 (t˜)ei(x1 (t˜)+( /2−2 1 )t˜) ,
A2 (t˜) = a2 (t˜)ei(x2 (t˜)+ t˜/2) , (2.3)
we obtain the following coupled equations
1 da1 3π 3 β π hˆ = a1 a22 sin(2(x2 − x1 )) + sin(2x1 )a1 − a1 γˆ dt 2ω1 2ω1 2 3π 3 β 2 dx1
= − + 2 1 − (a1 + 2a22 + a22 cos(2(x2 − x1 ))) dt 2 2ω1 π hˆ + cos(2x1 ) 2ω1 1 da2 3π 3 β π hˆ = a2 a21 sin(2(x1 − x2 )) + sin(2x2 )a2 − a2 γˆ dt 2ω1 2ω2 2
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
dx2
3π 3 β 2 =− − (a2 + 2a21 + a21 cos(2(x1 − x2 ))) dt 2 2ω2 π hˆ + cos(2x2 ), 2ω2
dφ = dt (2.4)
where for simplicity we have omitted the tilde in the time variable t˜. In our analysis herein we will again treat the drive hˆ and the damping γˆ in (2.4) as perturbations.
y2 ∂ H0 (x, y, I )
x2 = −δ I − − (3 − δ ) − (1 − δ ) ∂I 4 2 2
which has the Hamiltonian function
δ
x4 H0 (x, y, I, φ ) = − I2 − I + 3− 2 4 4
+ x2 y2 1 −
2.1. Unperturbed equations: setting damping and drive to zero
−
For γˆ = hˆ = 0, we obtain the unperturbed system
da1 dt dx1 dt da2 dt dx2 dt
2
− 2 1 −
2
−
ω1 2 a − a21 (2 + cos(2(x1 − x2 ))), ω2 2
under the scaling
a1 −→ a1
2ω2 , 3π 3 β
a2 −→ a2
(2.5)
2ω1 . 3π 3 β
Introducing the new variables (see [25])
B=
1 2 a , 2 1
1 2
θ = x1 − x2 , I = (a21 + a22 ), φ = x2 , δ =
ω1 , t → t/2 ω2
the unperturbed system above becomes
dB dt dθ dt dI dt dφ dt
= −2B(I − B ) sin(2θ )
= 1 − I (2 − δ + cos(2θ )) + B 4 −
δ2 + 1 + 2 cos(2θ ) δ
=0 = −δ I −
4
− B(2 − δ + cos(2θ ))
(2.6)
which can be written in a Hamiltonian form as
B˙ = −∂θ H˜ 0 ,
θ˙ = ∂B H˜ 0 , I˙ = −∂φ H˜ 0 , φ˙ = ∂I H˜ 0
with the Hamiltonian
(2.7)
2
2
dI ∂ H0 (x, y, I ) =− =0 dt ∂φ
is a compact normally hyperbolic invariant manifold with boundary that admits two homoclinic manifolds W0± . Solutions on are typically periodic, as one can readily verify by substituting x = y = 0 into the φ -components of (2.9). An exception occurs for the action value
I0 = −
4δ
,
which is labeled as C, a circle of fixed points, on the manifold . The fixed points in this circle represent a 1-parameter family of nonlinear normal modes. The 2D invariant manifold , has a 3D stable and unstable manifolds Ws (), Wu (), which coincide to form a 3D homoclinic manifold = W s () ∩ W u (). Trajectories on the homoclinic manifold are homoclinic orbits that connect the origin of the (x, y) plane to itself. Moreover, the constant value of the Hamiltonian along such an orbit is equal to its value at the origin, namely H0 (0, 0, I ) = H0 (x, y, I ) = H˜0 (B, θ , I ) which yields an equation for the homoclinic orbits in terms of the action-angle variables, i.e.
Bh (t, I ) =
2δ a2 q cosh(2at ) + p
θ (t, I ) = tan (2.8)
δ2 + 1 2 δ2 + 1 −x y 2− 2δ 2δ
dy ∂ H0 (x, y, I ) δ + 1 δ2 + 1 = = x3 1 − + xy2 2 − dt ∂x 2δ 2δ − x ( I ( 3 − δ ) − 1 )
(2.10)
1 = (x, y, I, φ ) : x = y = 0, 1 < δ < 3, I > 3−δ
h
At this point we have to emphasize that I and H˜0 are constants of the motion in the unperturbed system. Note that (B, θ , I, φ ) ∈ R+ × S × R+ × S, where S is the unit circle, and R+ are the non-negative reals. √ √ Now, introducing the Cartesian variables x = 2Bcosθ , y = 2Bsinθ , the system (2.7) becomes
dx ∂ H0 (x, y, I ) =− = −y3 1 − dt ∂y + y ( I ( 1 − δ ) − 1 )
x2 y2 ( I ( 3 − δ ) − 1 ) − ( I ( 1 − δ ) − 1 ) 2 2
δ2 + 1 H˜ 0 (B, θ , I ) = − I − I + B 2 − 2 4 δ − B(I (2 − δ ) − 1 ) − B(I − B )cos(2θ ) δ
δ 2 + 1 y4 δ2 + 1 + 1− 2δ 4 2δ
The point x = y = 0 is a fixed point of our unperturbed system, as given by Eqs. (2.9a) and (2.9b), and it is a saddle point for 1 < δ < 3 and I > 1 /3 − δ . This result comes from the fact that this fixed point is saddle for values of the positive constant of motion I, that satisfy the inequality [I (1 − δ ) − 1 ][I (3 − δ ) − 1 ] < 0. The set
ω2 2 a − a22 (2 + cos(2(x2 − x1 ))) ω1 1
= a21 a2 sin(2(x1 − x2 )) =−
δ2 + 1 4δ
(2.9d)
2.2. Analytical expressions for homoclinic orbits
= −a1 a22 sin(2(x1 − x2 )) =−
135
(2.9a)
−1
I ( δ − 3 ) + 1 tanh(at ) I ( 1 − δ ) − 1
−a(δ 2 − 1 ) χ1h (t, I ) =
tan−1 − +
p2 − q2
δI −
(2.11)
4
p−q tanh(at ) p+q
t + x1 ( 0 )
φ h (t, I ) = χ1h (t, I ) − θ h (t, I )
(2.12)
(2.13) (2.14)
where
q ≡ I ( δ 2 − 1 ) + 2 δ 1 > 0 p = − 1 (δ 2 − 4δ + 1 ) − I (δ 3 − 6δ 2 + 7δ − 2 )
2
a2 = 21 − 2I 1 (δ − 2 ) − I2 (δ − 3 )(δ − 1 ). (2.9b) (2.9c)
In Fig. 2(a,b), we sketch the scenario of the dynamics explained here. In Fig. 3 we plot the formation of the homoclinic orbit to . The homoclinic manifolds that we obtained in terms of the original amplitude equations (2.4) with hˆ = γˆ = 0 are shown in Fig. 4.
136
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
We assume that every resonant curve C is connected to some other such state for γˆ = 0, hˆ = 0, i.e. there are “single-jump” fast transitions between resonances [17,25]. More precisely, we assume that the stable and unstable directions associated with C give rise to a manifold W0 (C ) of orbit homoclinic to C, see Fig. 2. We shall consider homoclinic manifolds with two components W0+ (C ) and W0− (C ). The orbits in both W0+ (C ), W0− (C ) connect the same pair of points with the phase shift vector
x˜ = lim (xh (t ) − xh (−t )) t→∞
= (x, y, I, φ ) = (0, 0, 0, lim (φ h (t ) − φ h (−t ))) t→∞
being the same for any solution xh (t) in W0+ (C ) ∪ W0− (C ), where φ = x1 − θ . The orbits given by Eqs. (2.11)–(2.14) are homoclinic to (see Fig. 3). In Fig. 2(b) we sketch some of them. At resonance, for I = I0 , the orbits are heteroclinic, connecting fixed points that are φ apart, where
φ = x 1 − θ Fig. 2. (a) The homoclinic manifold of C. (b) A sketch of the assumption that the orbits in both W0+ (C ) and W0− (C ) connect the same pairs of points, i.e. the phase shift vector. (c) A 3− chain of unperturbed orbits. See Definition 2.1 for the meaning of xj (t), j = 1, 2, . . . . The pictures are reprinted from [17].
and
2a ( δ 2 − 1 ) −1 x 1 =
tanh
p2 − q2
θ = 2 tan
−1
p−q p+q
I0 (δ − 3 ) + 1 I0 (1 − δ ) − 1
Therefore
p − q 2a ( δ 2 − 1 ) I0 (δ − 3 ) + 1 −1 φ =
tanh − 2 tan−1 p+q I0 (1 − δ ) − 1 p2 − q2
B
10
5
0 -4
28.9 28.95 -2
29 29.05
0 2
29.1 29.15
4
I
Fig. 3. Orbits homoclinic to . For I = I0 (the orbit in the middle), dφ /dt = 0 on , and the orbit is heteroclinic, connecting fixed points on that are φ apart. For IࣤI0 , dφ /dtࣤ0 on . The parameters are δ = 1.49704, = −110, 1 = 27.6085.
Definition 2.1. Let us consider for our unperturbed system (2.9) a set of solutions {x j (t )}Nj=1 each homoclinic to the torus (circle) C with the property
I=50
lim x1 (t ) = b0 ,
t→−∞
10
a2
lim x j−1 (t ) = lim x j (t ),
t→∞
t→−∞
j = 2, . . . , N
and let {Pj }Nj=1 be a sequence of +1’s and −1’s. The solutions xj (t) are said to form an N− chain with base point b0 and jump sequence {Pj }Nj=1 if ∀t ∈ R we have
9 8 7
x (t ) ∈ j
I=20 0 a1 sin(x1-x 2)
5
0 a cos(x -x ) 1
1
5
I=30
6 10
Remark 1. In resonance problems the fixed points in C typically represent same solution with different initial conditions in a moving frame. In such a case the phase shift (φ ) on an orbit homoclinic to the underlying solution will be independent of the initial phase. One can built “chains” out of the solutions xh (t ) = (x, y, I, θ ) to connect fixed points that are Nφ apart for some integer N (see Fig. 2(b,c)). The members of such a chain can be taken from both components of the homoclinic manifold, and they connect N + 1 different resonant states altogether. We shall seek actual solutions staying close to these chains, i.e., N− chains of unperturbed orbits will serve as “shadowing sets” for real orbits exhibiting N excursions from the set of resonances. We first define an N− chain more precisely.
-5
-10
-5
W0+ (C ) W0− (C )
i f Pj = +1 i f Pj = −1
We denote an N − chain with base point b0 by XN (b0 ). 3. Homoclinic intersections in the damped, driven system
2
Fig. 4. Results of a numerical integration of Equations (2.4) with hˆ = γˆ = 0, ω1 = 1.58114, ω2 = 1.22474, 1 = 35 and = −150.
In this section we want to study how the drive and the damping as perturbations affect the dynamics. Writing the drive
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
hˆ = 2εω1 h¯ /π and the damping γˆ = 2ε γ¯ , ε 1, the perturbed system (2.4) in terms of the action-angle variables takes the form
dB = −2B(I − B ) sin(2θ ) + ε Bh sin(2(φ + θ )) − ε Bγ dt dθ δ2 + 1 = 1 − I (2 − δ + cos(2θ )) + B 4 − + 2 cos(2θ ) dt δ εh + (cos(2(φ + θ ))) − δ cos(2φ ) 2 dI = ε h(B sin(2(φ + θ )) + (I − B )δ sin(2φ )) − εγ I dt dφ
h = −δ I − − B(2 − δ + cos(2θ )) + ε δ cos(2φ ), (3.1) dt 4 2 where we have dropped the bar in h¯ and γ¯ for simplicity. This system can also be written as
dB ∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) =− + ε kB = − dt ∂θ ∂θ ∂ H˜ (B, θ , I, φ ) 1 B +ε − +g
∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) ∂ H˜ 1 (B, θ , I, φ ) + ε kθ = +ε ∂B ∂B ∂B ∂ H˜ (B, θ , I, φ ) ε − 1 + gI ∂φ dφ ∂ H˜ 0 (B, θ , I ) ∂ H˜ 0 (B, θ , I ) ∂ H˜ 1 (B, θ , I, φ ) = + ε kφ = +ε (3.2) dt ∂I ∂I ∂I where
h (B cos(2(φ + θ )) + δ (I − B ) cos(2φ )) 2
and the dissipative perturbations are given by gB = −Bγ and gI = −γ I. In Cartesian variables, the perturbed system takes the form
dx ∂ H0 (x, y, I ) ∂ H1 (x, y, I, φ ) =− +ε − + gx dt ∂y ∂y ∂ H (x, y, I, φ ) dy ∂ H0 (x, y, I ) 1 = +ε + gy dt ∂x ∂x ∂ H (x, y, I, φ ) dI 1 =ε − + gI dt ∂φ
dφ = dt
H1 (x, y, I, φ ) =
h x2 + y2 2
2
+δ I − gx = −
γx 2
,
gy = −
γy 2
x x2 + y2
x +y 2
2
2 ,
cos φ −
cos(2φ )
gI = −γ ,
dI = ε hIδ sin(2φ ) − εγ I, dt
(3.4)
dφ
ε hδ = −δ I − + cos(2φ ). dt 4 2
(3.5)
According [25] and [17] in order to study the slow dynamics, which is induced by the perturbation on Mε near resonance, we √ set a slow variable I = I0 + ε η into Eqs. (3.4) and (3.5) along √ with a slow time scale τ = εt to get
√ dη = (hδ sin(2φ ) − γ0 )(I0 + ε η ) dτ √ dφ εhδ = −δη + cos(2φ ) dτ 2
(3.6) (3.7)
(3.8)
dφ = −δη dτ
(3.9)
which can be obtained from the Hamiltonian
H (η , φ ) = −
y x2 + y2
sin φ
gφ = 0
For 0 < ε 1 the perturbed system (3.3):
p1c (ξ ) = 0,
of Cr -class and ⊂ Mε , while for ε = 0, we have M0 = . • Mε has codimension-one local stable and unstable manifolds s (M ) and W u (M ). Wloc ε ε loc Solutions of (3.3) admit two different time scales near : fast in x, y-variables and slow motion in (I, φ ).
2
+ I0 γ φ +
1 I0 hδ cos(2φ ) 2
(3.10)
1 −1 γ sin 2 hδ 1 2
π − sin−1
(3.11)
γ
(3.12)
hδ
Our goal is to prove the existence of generalized Shilnikov orbits homoclinic or heteroclinic to the slow sinks [17]. We first compute the energy-function N H associated with (3.10) [17]
N H = H ( φ0 + N φ ) − H ( φ0 ) −
N k=1
+∞
−∞
→ DH0 , − g dt
(3.13)
− → where g = (−H1y + gx , H1x + gy , −H1φ + gI , H1I + gφ ) is the perturbed vector field. In order to find the energy-function N H (3.13) we calculate the following
dH˜ 1 ∂ H0 x ∂ H0 y + g + g + dt ∂x ∂y dH˜ 1 dφ dθ =− −γI −γB dt dt dt
→ DH0 , − g=−
• near C there exists a unique, locally invariant manifold
Mε = {(x, y, I, φ ) : x = xε (y, I, φ )}
δη2
The phase portrait of Eqs. (3.6) and (3.7) in the new variables (η, φ ) for different values of ε is given in Fig. 5. The slow local Hamiltonian H admits two saddle-type fixed points p1s , p2s as well as two centers p1c , p2c given by
p1c (ξ ) = 0,
(3.3)
dη = (hδ sin(2φ ) − γ0 )I0 dτ
and
∂ H (x, y, I, φ ) ∂ H0 (x, y, I ) 1 +ε + gφ ∂I ∂I
where
By setting B = 0 in the two latests Equations of (3.1) we obtain the equations that describe the dynamics on Mε near the resonance I = I0 , i.e.
ε we obtain
dθ = dt dI = dt
3.1. Dynamics near resonance
From the leading terms in Eqs. (3.6) and (3.7) independent of
∂θ
H˜ 1 (B, θ , I, φ ) =
137
∂ H0 I g ∂I
where x1 = θ + φ , B(±∞ ) = 0 and dθ /dt = −I (2 − δ − cos(2θ )) −
1 and dx1 /dt = −Iδ − /4 − B(1 − δ /2δ ). Therefore, we obtain
−
N l=1
=−
+∞
−∞ N l=1
− → DH0 , g dt +∞
−∞
dH˜ 1 −
dt
−γI
dφ dθ −γB dt dt
dt
138
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
−
2ξ N hδ
− 2 φ − θ −
L˜ I0
(3.16)
If the dissipation parameter satisfies the condition
ξ<
( N φ / 2 ) sin(Nφ ) 1 + 2 sin sin(φ /2 ) 2N hδ
− 2 φ − θ −
L˜ I0
,
then the energy function admits transverse zeros, where
L˜ =
−2I0 (δ 2 − 1 ) + 2δ 1 tan−1 √ ( δ − 1 ) δ 2 − 6δ + 1
δ−1 δ 2 − 6δ + 1
√
I0 (δ − 3 ) + 1 I0 (1 − δ ) − 1
3.2. A homoclinic connection to the sink pε We are finally in a position to show the existence of an orbit homoclinic to the sink pε . To achieve this, we first show that there exists a homoclinic orbit that approaches pε asymptotically backward in time, and approaches the perturbed annulus Mε asymptotically forward in time. We then estimate the conditions under which the perturbed counterpart of the point, which is reached on forward in time in the unperturbed system, lies within the basin of attraction of the sink pε on Mε . This gives us an estimate for the possibility of obtaining a Shilnikov orbit that connects the sink back to itself. This first step is done by finding the conditions for which the energy function has simple zeros. Theorem 1. System (3.3) depends on a vector κ of system parameters. If (i) N H ( p1c (κ ), κ ) = 0 (ii) D pc N H ( p1c (κ ), κ ) = 0, Dξ N H ( p1c (κ ), κ ) = 0 (iii) The point (0, φc (ζ ) + Nφ ) lies in the domain of attraction of the sink. Then there exist a codimension-one set M+ such that: (a) Our system admits an N−pulse Shilnikov-type orbit x , homoclinic to the point p . The takeoff point for this orbit is p , and the jump sequence of the orbits is given by Fig. 5. (a) Numerical phase portraits of Equations (3.6) and (3.7) with ε = 0. (b) Numerical phase portrait of Equations (3.6) and (3.7) with ε = 1. The parameters are δ = 1.55, = −150, h = 1, ζ = 0.45, γ = hζ δ .
=
hδ I0 2
sin
N φ 2
sin 2φ
2 sin(2φ ) sin(Nφ ) + γ I0 Nφ (3.14)
where I0 = − /4δ is the resonance, N is the number of homoclinic jumping, φ = x1 − θ and
2I0 (δ 2 − 1 ) + 2δ 1 L˜= − tan−1 √ ( δ − 1 ) δ 2 − 6δ + 1
δ−1 √ 2 δ − 6δ + 1
I0 (δ − 3 ) + 1 . I0 (1 − δ − 1 )
Next, consider again and calculate the term H (φ0 + Nφ ) − H (φ0 ) on the slow manifold, i.e.
Remark 2. According to the Theorem 1 we find solutions of N H ( p1c , κ ) = 0 where p1c center point and κ the parameters of system. Some elementary algebra turns this condition (N H ( p1c , κ ) = 0) into
1 − ζ 2 (1 − cos(2N φ ) ) = −ζ sin(2N φ )
N H
H ( φ0 + N φ ) − H ( φ 0 ) =
hI0 δ (cos(2(φ0 + Nφ )) 2 − cos(2φ )) + γ I0 Nφ .
From (3.14) and (3.15), we obtain
N H (φ ) = cos(2(φ0 + Nφ )) − cos(2φ0 ) sin + sin
N φ 2
φ 2
2 sin(2φ0 ) sin(Nφ )
(3.15)
j = 1, ..., N − 1.
(b) There also exists another codimension-one set M− that yields similar homoclinic orbits with jump sequence {−Pj }Nj=1 . (c) Both surfaces M+ and M− have open neighborhoods in the (κ , ) parameter space for which our system admits Smale horseshoes [17].
+ I0 N θ γ + N L˜γ
Pj+1 = σ sign( j H ( p1c , κ )Pj ,
P1 = +1,
sin +2
N φ 2
sin 2φ
1 sin sin−1 (ζ ) sin(Nφ ) 2
− 2 − 2 φ − θ −
1 L˜ Nζ I0
(3.17)
where ζ = ξδ .The Eq. (3.17) gives the perturbation parameter value
ζ=
1 − cos(2Nφ )
( N φ / 2 ) (1 − cos(2Nφ ))2 + −2BˆN + 2 sin sin(Nφ ) − sin(2Nφ ) sin(φ /2 )
2
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
139
(3.18) where
Bˆ = −2φ − θ −
1 L˜. I0
3
Of course, this result is meaningful only for nonzero values
k ∈ Z.
2
(3.19)
The second transversality condition in Theorem 1 gives
D ξ [ H ( N
p1c ,
κ )] = 0
B
φ
kπ = , N
1 0
(3.20)
20
1.5
to be satisfied whenever (3.18) and (3.19) hold. Carrying out the differentiation, we obtain that (3.20) is violated if
1 − ζ 2 −sin(2N φ ) + 2
N φ
sin
1
−1
sin sin 2 sin 2φ 2
(ζ ) sin(Nφ )
1 − 2 − 2φ − θ − L˜ N I0 = ζ (1 − cos(2Nφ ))
(3.21)
However, under condition (3.19), Eqs. (3.17) and (3.21) cannot be satisfied simultaneously. Therefore, the second condition of (ii) of Theorem 1 is satisfied. The first condition in (ii) is satisfied by inequality (2.6). It remains to verify condition (iii) of Theorem 1 to ensure that the landing point of N-pulse orbit taking off from a slow sink lies in the domain of attraction of one of the sinks. According to [17], the landing point lies near the angle φc (ζ ) + Nφ , where φ c (ζ ) is the φ coordinate of the center p1c . For small enough ε , the landing point will fall in the basin of a slow sink if the point
( 0 , φc ( ζ ) + N φ ) falls in D01 . This can be verified by comparing the energy of the approximate landing point to the “energy” H (η, φ ) of the boundary of D01 ∪ D02 . Note, however, that this energy is not well-defined, since H (η, φ ) is only a local Hamiltonian. To obtain a uniquely defined energy, we consider the point that lies in the interval [ π2 , 32π ] and is kπ apart from the approximation landing point. The φ -coordinate of this point can be computed as
φmN (ζ ) = φs (ζ ) + φc (ζ ) + Nφ − φs (ζ ) modπ ,
(3.22)
where φ c (ζ ), φ s (ζ ) are φ -components of center and saddle points. N (ζ ) > φ (ζ ), then redefine φ N (ζ ) by subtracting π from it. If φm s m Now, if N H ( 0 , φm (ζ )) > H (0, φs (ζ ))
cos(2φs (ζ )) − cos(2φ (ζ )) > 2ζ (φs (ζ ) − φ (ζ )) N m
(3.24)
:= ξδ , ξ =
ζ γ /h then (0, φmN (ζ )) ∈ D01 , which implies that (0, φs (ζ ) + Nφ ) ∈ D01 ∪ D02 . Now from the above we want to write an approximate condition which guarantees that this orbit approaches pε as t → −∞. To succeed in this, we find a condition for which the unperturbed heteroclinic orbit, which asymptotes to pc as t → +∞ returns back to the point on the circle of fixed points that is inside the homoclinic separatrix loop connecting the saddle qc to itself. This kind of an orbit is shown in Fig. 6. This condition is presented in terms of the difference φ between the asymptotic values of the angular variable φ as
φs < φc + φ < φ m
0.5 0 30 -0.5
Fig. 6. The heteroclinic orbit given by Equations (2.11) and (2.14) with I = I0 , superimposed with the phase portrait of the unperturbed scaled system on ε near resonance, given by Equations (3.8) and (3.9). The parameters are δ = 1.55, = −150, 1 = 35, ζ = 0.45, h = 1 and γ = hζ δ . This value of ζ sets φc = 0.233383 and we can see from the figure that φs < φc + φ < φm .
Eqs. (3.18) and (3.25) define the conditions which are responsible for the existence of orbits homoclinic to the sink pε . What we want now is to relate these results to the actual physical parameters of the coupled resonators. As we have already defined δ , the value of λ is then given by λ = (−δ 2 − 2π βδ 2 + 1 + 2π β )
/(δ 2 + 1 ). The scaled frequency
1 is 1 = (ω1 − ω2 )/ = ( 1 + 2π β − λ − 1 + 2π β + λ)/2 , so after fixing it is also determined by δ . The ratio ζ = γ /(hδ ), between the damping coefficient and the drive amplitude, has to be positive in order for the damping coefficient γ to be positive and has the standard physical meaning of energy dissipation. The parameter ζ is positive if the inequality,
−2BˆN − sin(2Nφ ) + 2
sin(Nφ /2 ) sin(Nφ ) > 0 sin(φ /2 )
(3.26)
is satisfied. The φ values of the fixed points are shown in Fig. 7(a) along with φc + φ and φ m for a particular value of . The parameter values for which these φ values satisfy the condition (3.25) are displayed in Fig. 7(b), which outlines the values of the coupling and driving frequency, for which orbits homoclinic to the sink pε exist. 4. Conclusion
(3.23)
or N m
25
J
1
(3.25)
where φ m is the maximal value of φ on the homoclinic orbit, connecting the saddle qc to itself.
In this work we studied chaotic dynamics of a pair of parametrically-driven coupled SQUIDs arranged in series and we found under which conditions it can exist by using high dimensional Melnikov theory. This was achieved by applying a method of Kovacic and Wiggins on transformed amplitude equations that were derived from the equations of motion, which model an actual experimental realization of coupled SQUIDs in metamaterials. We considered the amplitude of the drive and the damping to be small and obtained explicit expressions for orbits homoclinic to a twodimensional invariant annulus in the unperturbed equations. At resonance, we were able to calculate the energy function analytically, and provide a primary condition for having homoclinic orbits in the full, perturbed equations. By further studying the effects of perturbations on the invariant annulus near resonance, we found a secondary condition for the existence of orbits homoclinic to a fixed point of a saddle-focus type. We used a numerical scheme to verify our theoretical predictions. Such Shilnikov homoclinic orbits give rise to a particular type of horseshoe chaos, which can be expected in the dynamics of the full system for parameter values
140
M. Agaoglou et al. / Chaos, Solitons and Fractals 99 (2017) 133–140
Fig. 7. (a) The values of φs , φc , φc + φ and φ m as functions of δ , for = −150. For δ > 1.524 the condition is satisfied. (b) Parameter values for which the condition is satisfied are indicated in gray. In both figures = 0.01.
in the vicinity of those presented here. Our work provides parameter regions where homoclinic chaos can be observed that may be exploited for, e.g., chaos communication. As future work we can consider the possibility to have two SQUIDs with one above the other in which case the coupling is stronger and has opposite sign. Acknowledgments M.A. thanked Professor Ron Lifshitz for fruitful discussions. The work of M.A. has been co-financed from resources of the operational program “Education and Lifelong Learning” of the European Social Fund and the National Strategic Reference Framework (NSRF) 2007-2013 within the framework of the Action State Scholarships Foundation’s (IKY) mobility grants programme for the short term training in recognized scientific/research centers abroad for candidate doctoral or postdoctoral researchers in Greek universities or research centers. The work of V.R. has been co-financed by the European Union (European Social Fund ESF) and Greek national funds through the Operational Program “Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF) Research Funding Program D.534 MIS: 379337: THALES Investing in knowledge society through the European Social Fund. V.R and M.A acknowledge support from FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-606096). H.S. is supported by the Engineering and Physical Sciences Research Council (Grants No. EP/M024237/1). References [1] Josephson BD. Possible new effects in superconductive tunnelling. Phys Lett 1962;1:251–3. [2] Clarke J, Braginski AI, editors. The SQUID handbook. Fundamentals and technology of SQUIDS and SQUID systems, vol. I. Weiheim: Wiley-VCH; 2004. [3] Likharev KK. Dynamics of Josephson junctions and circuits. Gordon and Breach, Philadelphia; 1986. [4] Fesser K, Bishop AR, Kumar P. Chaos in rf SQUIDs. Appl Phys Lett 1983;43:123–4. [5] Bocko MF, Herr AM, Feldman MJ, M J. Prospects for quantum coherent computation using superconducting. IEEE Trans Appl Supercond 1997;7:3638–41.
[6] Yamashita T, Takahashi S, Maekawa S. Superconducting pi qubit with three josephson junctions. Appl Phys Lett 2006;88:132501–1-132501-3. [7] Lazarides N, Tsironis GP, G P. rf superconducting quantum interference device metamaterials. Appl Phys Lett 2007;16:163501–163501-3. [8] Lazarides N, Tsironis GP, Eleftheriou N. Dissipative discrete breathers in rf SQUID metamaterials. Nonlinear Phenom. Complex Syst. 2008;11:250–8. [9] Yen TJ, Padilla WJ, Fang N, Vier DC, Smith DR, Pendry JB, et al. Terahertz magnetic response from artificial materials. Science 2004;303:1494–6. [10] Caputo JG, Gabitov I, Maimistov AI. Electrodynamics of a split-ring josephson resonator in a microwave line. Phys Rev B 2012;85:205446. [11] Linden S, Enkrich C, Dolling G, Klein MW, Zhou JF, Koschny T, et al. Photonic metamaterials: magnetism at optical frequencies. IEEE J Selec Top Quant Electron 2006;12:1097–105. [12] Chen HT, Yang H, Singh R, OHara JF, Azad AK, Stuart A, et al. Tuning the resonance in hightemperature superconducting terahertz metamaterials. Phys Rev Lett 2010;105:247402. [13] Overhauser AW. Temperature dependence of the josephson current. Phys Rev Lett 1999;83:180. [14] Han S, Wolf EL. Temperature dependence of the critical josephson current in superconductor insulatornormal-metal proximity junctions near Tc . Phys Rev B 1987;35:4669. [15] Hikino S-i, Mori M, Takahashi S, Maekawa S. Temperature dependence of josephson current in a superconductor/ ferromagnet/superconductor junction. Physica C 2007;460-462:1323–4. [16] Veldhorst M, Snelder M, Hoek M, Gang T, Guduru VK, Wang XL, et al. Josephson supercurrent through a topological insulator surface state. Nat. Mater. 2012;11:417–21. [17] Haller G. Chaos near resonance. Springer, New York; 1999. [18] Menon G, Haller G, Rothos V. Silnikov manifolds in coupled nonlinear schrödinger equations. Phys Lett A 1999;263:175–85. [19] Kenig E, Tsarin YA, Lifshitz R. Homoclinic orbits and chaos in a pair of parametrically driven coupled nonlinear resonators. Phys Rev E 2011;84:016212. [20] Zhang W, Huang YT, Yao MH. Multi-pulse homoclinic orbits and chaotic dynamics of a parametrically excited nonlinear nano-oscillator with coupled cubic nonlinearities. Sci China Phys Mech Astron. 2014;57:1098–110. [21] Dana SK, Sengupta DC, Edoh KD. Chaotic dynamics in josephson junction. IEEE Trans on Circ and Syst I 2001;48:990–6. [22] Myneni K, Barr TA, Reed BR, B R. High-precision ranging using a chaotic laser pulse train. Appl Phys Lett 2001;78:1496–8. [23] Hand LN, Finch JD. Analytical mechanics. Cambridge Univ. Press, Cambridge; 1998. [24] Strogatz SH. Nonlinear dynamics and chaos. Addison-Wesley, Reading MA; 1994. [25] Kovacic G, Wiggins S, S. Orbits homoclinic to resonances, with an application to chaos in a model of the forced and damped sine-gordon equation. Physica D 1992;57:185–225.