# Influence of Bifurcation Structures Revealed by Refinement of a Nonlinear Conductance in JosephsonJunction Element

Authors: Yuu Miino; Tetsushi Ueta
Journal: Complexity
 (2018)

---

## Abstract

We conduct a bifurcation analysis of a single-junction superconducting quantum interferometer with an external flux. We approximate the current-voltage characteristics of the conductance in the equivalent circuit of the JJ by using two types of functions: a linear function and a piecewise linear (PWL) function. We describe a method to compute the local stability of the solution orbit and to solve the bifurcation problem. As a result, we reveal the bifurcation structure of the systems in a two-dimensional parameter plane. By making a comparison between the linear and PWL cases, we find a difference in the shapes of their bifurcation sets in the two-dimensional parameter plane even though there are no differences in the one-dimensional bifurcation diagrams or the trajectories. As for the influence of piecewise linearization, we discovered that grazing bifurcations terminate the calculation of the local bifurcations, because they drastically change the stability of the periodic orbit.


---

## Body

## 1. Introduction

Josephson junctions (JJs) are devices composed of two superconductors coupled with a weak link. JJs have an extraordinary current-voltage characteristic, and circuits incorporating them show a plentiful variety of nonlinear phenomena. For example, Salam et al. [1] discovered chaotic responses in a JJ circuit and a forced JJ circuit. Hadley et al. [2] found phase locking of JJ series arrays, while Cirillo and Pedersen [3] studied bifurcation phenomena and chaos in the response of JJs. Solving the bifurcation problem is important for comprehending the properties of the system, but most of the previous studies failed to solve it or did so imprecisely, e.g., Dana et al. [4] suggested a simulation of JJ circuits defined by the piecewise linear conductance but did not solve its bifurcation problem.On the other hand, hybrid dynamical systems (HDSs) have been studied intensively by many researchers [5–7]. An HDS combines a continuous-time dynamical system and a discrete-time dynamical system, e.g., a circuit including a switch like a DC/DC converter [8], a neuron model including a firing scheme like the Izhikevich neuron model [9], or a conflict system like a bell [10]. HDSs embody a rich variety of bifurcation phenomena reflecting their interrupting properties and the discontinuity of their derivatives, and many bifurcation analysis tools and applications have been developed for them. For example, Kousaka et al. [11] studied the periodic orbits in an autonomous HDS and proposed the method to compute the local stability of their orbits and local bifurcation sets. Our previous study [12] suggested a scheme to apply Kousaka’s method to the nonautonomous HDS. Piiroinen et al. [13] discovered chaotic behavior and grazing bifurcations occurring in an HDS. Ito et al. [14] suggested a method to control chaos in HDS by perturbing the threshold value of the system.In this study, we solve the bifurcation problem of a single-junction superconducting quantum interferometer with an external flux. We assume two types of conductance in the equivalent circuit of the JJ: a linear conductance and a PWL conductance. By defining a PWL function that models the actual response of the JJ, we expect we can determine its properties. We use the HDS approach [12] to analyze the PWL system. In what follows, we define the system and its mathematical characteristics (Sections 2.1-2.2). We then describe the single-junction superconducting quantum interferometer with a vibrating external flux [15, 16]; we derive its circuit equation and normalized equation. We define the PWL function from the current-voltage characteristics of the JJ. We also consider a linear conductance, because there are some studies that use this approximation [17]. Next, we explain the local stability and local bifurcation phenomena of the periodic orbit (Section 2.3). We introduce the Poincaré map and variational equation and present criteria under which local bifurcations arise. We describe the criteria for a grazing bifurcation as well (Section 2.4). As the main result (Section 3), we reveal the bifurcation structure of the system in a two-dimensional parameter space. We discuss the bifurcations observed in the system by using one-dimensional bifurcation diagrams and trajectories. We make a comparison between the case of a linear function and the case of a PWL function, identifying the similarities and differences between them. Finally, we conclude this study (Section 4).


## 2. Materials and Methods

### 2.1. Single-Junction Superconducting Quantum Interferometer with the External Flux

Let us consider a single-junction superconducting quantum interferometer (SQUID) with an external flux [15, 16], as shown in Figure 1(a). This is also called an RF-SQUID core [18]. Taking Bext as the sum of DC and AC components [16], we get the equivalent circuit shown in Figure 1(b).Figure 1
(a) Single-junction superconducting quantum interferometer with external fluxBext. (b) An equivalent circuit of (a).
(a) (b)The circuit equations for the circuit in Figure1 are obtained by using Kirchhoff’s laws and the characteristics of each element:(1)iC+iG+iJ+iL=jt,iC=Cdvdt,iG=Gv,iJ=ICsin⁡ϕ,ħ2e·dϕdt=v,LdiLdt=v,where IC is the junction critical current, G(v) represents the voltage-current characteristic of the conductance (Figure 2(a) corresponds to the characteristic), ϕ is the phase difference of the superconducting order parameter across the junction, and ħ is the reduced Planck constant. Combining these equations yields the circuit equation of the circuit in Figure 1:(2)Cħ2e·d2ϕdt2+Għ2e·dϕdt+ICsin⁡ϕ+1L·ħ2eϕ=jt. In dimensionless form, assuming jt=I0+Icos⁡ωt, (2) becomes as follows:(3)dxdτ=y,dydτ=-kfy-cx-sin⁡x+B0+Bcos⁡ντ,where(4)x=ϕ,y=1ω0·dϕdt,τ=ω0t,ω0=2eICħC,k=1IC,fy=Gv,c=1LCω02,B0=I0IC,B=IIC,ν=ωω0.We can write the state variables x and y in vector notation: x∈R2.Figure 2
(a) Modeled conductance characteristic in the Josephson junction equivalent circuit. (b) Approximations of (a) by a piecewise linear function.
(a) (b)f is a function of y that corresponds to the conductance G(v) in the JJ equivalent circuit. The response of f is precisely modeled from experiments, as shown in Figure 2(a). In this study, we approximate f by a linear function and a piecewise linear (PWL) function. The linear approximation is an ideal case, for which we can use the legacy [19] method to analyze this smooth system. Conversely, the PWL approximation, as shown in Figure 2(b), is a realistic case; we should use the HDS approach [12] to analyze this nonsmooth system. For the linear case, we can use the following simple linear function of y:(5)f1y=y.For the PWL case, let us firstly define the state spaces as follows:(6)M1=x=x,y∈R2∣y≤-θ2,M2=x=x,y∈R2∣-θ2<y<-θ1,M3=x=x,y∈R2∣-θ1≤y≤θ1,M4=x=x,y∈R2∣θ1<y<θ2,M5=x=x,y∈R2∣θ2≤y.The PWL function in Figure 2(b) is as follows:(7)f2y=yifx∈M1∪M5,θ2y+θ1θ2-θ1ifx∈M2,0ifx∈M3,θ2y-θ1θ2-θ1ifx∈M4,where θi for i=1,2 corresponds to thresholds dividing f(y) into five segments. In this study, we will avoid dealing with the hysteresis phenomenon arising in the JJ.Anorbitx(t)∈R2 is the solution of (3) together with an initial condition x0=(x0,y0)∈R2. We define x(t) to be a function of time t and x0:(8)xt=φt,x0,x0=φ0,x0=x0.We obtain the orbits by using the Runge-Kutta method. In particular, we used a step of 10-2 for the integration in all simulations. We use the bisection method with a precision of 10-10 to calculate the exact time when the orbit collides with the boundary between Mj and Mk, namely, the time when y=±θi. An orbit x(t) isperiodic if(9)xt+L=xt,∀t∈R,∃L∈R.The period of this orbit is the minimum value of L.Since one of the state variables in (3) corresponds to an angle, the motion of this system evolves ona cylindrical surface. Let S1 be the set of points on a unit circle in R2:(10)S1=cos⁡x,sin⁡x∈R2∣x∈R,where x gives the angle of a point in S1 from the reference direction; accordingly, we define the cylindrical surface, (11)M=S1×R.Given a map(12)M⟶R2\0;x,y⟼r,θ=y,x,where y is scaled and biased appropriately to fit the polar coordinate system, and the cylindrical surface is homeomorphic to R2\{0}. We will thus observe the trajectory in the punctured plane R2\{0}.


### 2.2. Symmetry

Equation (3) is invariant to the following linear transformation:(13)xytB0⟼2nπ-x-yt+π2ncπ-B0,where n is an arbitrary integer. This means that we can transform a solution with B0=B0∗ and a solution with B0=2ncπ-B0∗ into each other by using the linear transformation (13). If the solutions are periodic, their stabilities are completely the same. Therefore, the bifurcation sets of the system (3) are symmetric with respect to the operation 2ncπ-B0⟼B0; that is, the bifurcation sets are reflectively symmetric about the line B0=ncπ.


### 2.3. Local Stability and Bifurcations

When discussing the asymptotic stability of the periodic orbit, we generally usethe Poincaré map [20] of the orbit, which is a discrete map as follows:(14)T:R2→R2;x0⟼Tx0=φ2π,x0.In the field of discrete-time dynamical systems, the point p∈R2 is calleda fixed point of the map T if(15)Tp=p.Similarly, the point p is calledanl-periodic point of the map T if(16)Tlp=p.These periodic points of T exactly correspond to the periodic orbit of the system (3).The variational equation from the fixed pointp with respect to the Poincaré map T is as follows:(17)un+1=∂T∂x0x0=pun,n=0,1,…where un is the difference between p and p′=Tn(p+u0). The multipliers μ1 and μ2 are obtained by solving the characteristic equation:(18)det∂T∂x0x0=p-μI=0,where I is a 2×2 identity matrix. When ∀i,μi≠0 and |μi|≠1, p is calleda hyperbolic fixed point. Assuming that μ1<μ2, let us classify the stability of a hyperbolic fixed point on the basis of its characteristic multipliers: (1)
A fixed point iscompletely stable if its multipliers satisfy |μ1|<1 and |μ2|<1; accordingly, we label the point D0.(2)
A fixed point isdirectly unstable if its multipliers satisfy 0<μ1<1<μ2; we label the point D1.(3)
A fixed point isinversely unstable if its multipliers satisfy μ1<-1<μ2<0; we label the point I1.(4)
A fixed point iscompletely unstable if its multipliers satisfy |μ1|>1 and |μ2|>1; we label the point D2.The index at the bottom left of the symbol indicates the number of unstable dimensions for the fixed point. Generally,D0 is calleda node, D1 and I1 are calledsaddles, and D2 is calleda source. Similarly to the fixed point, we classify the stability of the l-periodic point by using symbols, i.e., D0l, etc.We numerically obtain the pointp satisfying (15) or (16) by applying Newton’s method. We also numerically obtain the coefficient matrix of (17) by using the method mentioned in [12], which obtains the Jacobian matrix of the hybrid dynamical systems.By perturbing some of the parameters, we find the parameter sets where the stability ofp changes. This change of the stability is calleda local bifurcation of p [21], and each of the parameter sets is calleda local bifurcation set. Local bifurcations occur when one of the multipliers of p exceeds unity as a result of changing a parameter. In the system (3), two local bifurcations are possible:the tangent bifurcation (μi=1) andthe period-doubling bifurcation (μi=-1). We numerically obtain the local bifurcation sets by using the method in [12, 19].


### 2.4. Grazing Bifurcation

Let us consider the case that a periodic orbit passes near the borderx=θi. By perturbing a parameter, the periodic orbit might approach the border and finally graze it. After that, the periodic orbit cannot keep its topological structure anymore and completely disappears. This disappearance, or inversely appearance, is calleda grazing bifurcation of the periodic orbit [13], and a parameter set where a grazing bifurcation phenomenon arises is calleda grazing bifurcation set. After a grazing bifurcation arises, there might remain a periodic orbit whose shape is similar to the orbit that disappeared. In other words, such a remaining orbit has an exactly different stability from the disappeared orbit.The condition for a grazing bifurcation to occur is that a periodic orbit and the border of the system are tangent to each other at the timeτ: (19)q±ixτ=0,dq±idtt=τ=0,where q±i(x)=x±θi=0 corresponds the border of the system. We compute grazing bifurcation sets by simultaneously solving (15) and (19) with Newton’s method.


## 2.1. Single-Junction Superconducting Quantum Interferometer with the External Flux

Let us consider a single-junction superconducting quantum interferometer (SQUID) with an external flux [15, 16], as shown in Figure 1(a). This is also called an RF-SQUID core [18]. Taking Bext as the sum of DC and AC components [16], we get the equivalent circuit shown in Figure 1(b).Figure 1
(a) Single-junction superconducting quantum interferometer with external fluxBext. (b) An equivalent circuit of (a).
(a) (b)The circuit equations for the circuit in Figure1 are obtained by using Kirchhoff’s laws and the characteristics of each element:(1)iC+iG+iJ+iL=jt,iC=Cdvdt,iG=Gv,iJ=ICsin⁡ϕ,ħ2e·dϕdt=v,LdiLdt=v,where IC is the junction critical current, G(v) represents the voltage-current characteristic of the conductance (Figure 2(a) corresponds to the characteristic), ϕ is the phase difference of the superconducting order parameter across the junction, and ħ is the reduced Planck constant. Combining these equations yields the circuit equation of the circuit in Figure 1:(2)Cħ2e·d2ϕdt2+Għ2e·dϕdt+ICsin⁡ϕ+1L·ħ2eϕ=jt. In dimensionless form, assuming jt=I0+Icos⁡ωt, (2) becomes as follows:(3)dxdτ=y,dydτ=-kfy-cx-sin⁡x+B0+Bcos⁡ντ,where(4)x=ϕ,y=1ω0·dϕdt,τ=ω0t,ω0=2eICħC,k=1IC,fy=Gv,c=1LCω02,B0=I0IC,B=IIC,ν=ωω0.We can write the state variables x and y in vector notation: x∈R2.Figure 2
(a) Modeled conductance characteristic in the Josephson junction equivalent circuit. (b) Approximations of (a) by a piecewise linear function.
(a) (b)f is a function of y that corresponds to the conductance G(v) in the JJ equivalent circuit. The response of f is precisely modeled from experiments, as shown in Figure 2(a). In this study, we approximate f by a linear function and a piecewise linear (PWL) function. The linear approximation is an ideal case, for which we can use the legacy [19] method to analyze this smooth system. Conversely, the PWL approximation, as shown in Figure 2(b), is a realistic case; we should use the HDS approach [12] to analyze this nonsmooth system. For the linear case, we can use the following simple linear function of y:(5)f1y=y.For the PWL case, let us firstly define the state spaces as follows:(6)M1=x=x,y∈R2∣y≤-θ2,M2=x=x,y∈R2∣-θ2<y<-θ1,M3=x=x,y∈R2∣-θ1≤y≤θ1,M4=x=x,y∈R2∣θ1<y<θ2,M5=x=x,y∈R2∣θ2≤y.The PWL function in Figure 2(b) is as follows:(7)f2y=yifx∈M1∪M5,θ2y+θ1θ2-θ1ifx∈M2,0ifx∈M3,θ2y-θ1θ2-θ1ifx∈M4,where θi for i=1,2 corresponds to thresholds dividing f(y) into five segments. In this study, we will avoid dealing with the hysteresis phenomenon arising in the JJ.Anorbitx(t)∈R2 is the solution of (3) together with an initial condition x0=(x0,y0)∈R2. We define x(t) to be a function of time t and x0:(8)xt=φt,x0,x0=φ0,x0=x0.We obtain the orbits by using the Runge-Kutta method. In particular, we used a step of 10-2 for the integration in all simulations. We use the bisection method with a precision of 10-10 to calculate the exact time when the orbit collides with the boundary between Mj and Mk, namely, the time when y=±θi. An orbit x(t) isperiodic if(9)xt+L=xt,∀t∈R,∃L∈R.The period of this orbit is the minimum value of L.Since one of the state variables in (3) corresponds to an angle, the motion of this system evolves ona cylindrical surface. Let S1 be the set of points on a unit circle in R2:(10)S1=cos⁡x,sin⁡x∈R2∣x∈R,where x gives the angle of a point in S1 from the reference direction; accordingly, we define the cylindrical surface, (11)M=S1×R.Given a map(12)M⟶R2\0;x,y⟼r,θ=y,x,where y is scaled and biased appropriately to fit the polar coordinate system, and the cylindrical surface is homeomorphic to R2\{0}. We will thus observe the trajectory in the punctured plane R2\{0}.


## 2.2. Symmetry

Equation (3) is invariant to the following linear transformation:(13)xytB0⟼2nπ-x-yt+π2ncπ-B0,where n is an arbitrary integer. This means that we can transform a solution with B0=B0∗ and a solution with B0=2ncπ-B0∗ into each other by using the linear transformation (13). If the solutions are periodic, their stabilities are completely the same. Therefore, the bifurcation sets of the system (3) are symmetric with respect to the operation 2ncπ-B0⟼B0; that is, the bifurcation sets are reflectively symmetric about the line B0=ncπ.


## 2.3. Local Stability and Bifurcations

When discussing the asymptotic stability of the periodic orbit, we generally usethe Poincaré map [20] of the orbit, which is a discrete map as follows:(14)T:R2→R2;x0⟼Tx0=φ2π,x0.In the field of discrete-time dynamical systems, the point p∈R2 is calleda fixed point of the map T if(15)Tp=p.Similarly, the point p is calledanl-periodic point of the map T if(16)Tlp=p.These periodic points of T exactly correspond to the periodic orbit of the system (3).The variational equation from the fixed pointp with respect to the Poincaré map T is as follows:(17)un+1=∂T∂x0x0=pun,n=0,1,…where un is the difference between p and p′=Tn(p+u0). The multipliers μ1 and μ2 are obtained by solving the characteristic equation:(18)det∂T∂x0x0=p-μI=0,where I is a 2×2 identity matrix. When ∀i,μi≠0 and |μi|≠1, p is calleda hyperbolic fixed point. Assuming that μ1<μ2, let us classify the stability of a hyperbolic fixed point on the basis of its characteristic multipliers: (1)
A fixed point iscompletely stable if its multipliers satisfy |μ1|<1 and |μ2|<1; accordingly, we label the point D0.(2)
A fixed point isdirectly unstable if its multipliers satisfy 0<μ1<1<μ2; we label the point D1.(3)
A fixed point isinversely unstable if its multipliers satisfy μ1<-1<μ2<0; we label the point I1.(4)
A fixed point iscompletely unstable if its multipliers satisfy |μ1|>1 and |μ2|>1; we label the point D2.The index at the bottom left of the symbol indicates the number of unstable dimensions for the fixed point. Generally,D0 is calleda node, D1 and I1 are calledsaddles, and D2 is calleda source. Similarly to the fixed point, we classify the stability of the l-periodic point by using symbols, i.e., D0l, etc.We numerically obtain the pointp satisfying (15) or (16) by applying Newton’s method. We also numerically obtain the coefficient matrix of (17) by using the method mentioned in [12], which obtains the Jacobian matrix of the hybrid dynamical systems.By perturbing some of the parameters, we find the parameter sets where the stability ofp changes. This change of the stability is calleda local bifurcation of p [21], and each of the parameter sets is calleda local bifurcation set. Local bifurcations occur when one of the multipliers of p exceeds unity as a result of changing a parameter. In the system (3), two local bifurcations are possible:the tangent bifurcation (μi=1) andthe period-doubling bifurcation (μi=-1). We numerically obtain the local bifurcation sets by using the method in [12, 19].


## 2.4. Grazing Bifurcation

Let us consider the case that a periodic orbit passes near the borderx=θi. By perturbing a parameter, the periodic orbit might approach the border and finally graze it. After that, the periodic orbit cannot keep its topological structure anymore and completely disappears. This disappearance, or inversely appearance, is calleda grazing bifurcation of the periodic orbit [13], and a parameter set where a grazing bifurcation phenomenon arises is calleda grazing bifurcation set. After a grazing bifurcation arises, there might remain a periodic orbit whose shape is similar to the orbit that disappeared. In other words, such a remaining orbit has an exactly different stability from the disappeared orbit.The condition for a grazing bifurcation to occur is that a periodic orbit and the border of the system are tangent to each other at the timeτ: (19)q±ixτ=0,dq±idtt=τ=0,where q±i(x)=x±θi=0 corresponds the border of the system. We compute grazing bifurcation sets by simultaneously solving (15) and (19) with Newton’s method.


## 3. Results and Discussion

Let us choose the parameter values as follows:(20)k=0.2,c=0.4,ν=1.0,θ1=1.0,θ2=1.01.We will observe at most three equilibrium points if we choose c=0.4, as mentioned in [22]. (Parameter α in [22] corresponds to c.) We further assume that the parameters B0 and B are controllable. Now let us observe the bifurcation phenomena arising from the perturbations of these parameters. The bifurcation diagrams presented below contain notation used for expressing the bifurcation sets: Gl, which is plotted as a magenta solid curve, represents the tangent bifurcation set of an l-periodic point; Il, which is plotted as a blue solid curve, represents the period-doubling bifurcation set of an l-periodic point; Z1, which is plotted as a black dashed curve, represents the grazing bifurcation set of an l-periodic orbit. We plot the solution orbits in the punctured plane R2\{0}. Here, we assume that Mi′ is a subset of R2\{0} and corresponds to Mi for i=1…5. To obtain a 1-dimensional bifurcation diagram, we compute images of a point under 10,000 iterations of the Poincaré maps and take the last 500 points of the data to plot. We usethe maximum Lyapunov exponent that indicates whether the trajectory is chaotic [1, 3, 13, 19, 20] or not. We can recognize that the trajectory seems to be chaotic if the exponent is positive. We approximately compute the exponent with 10,000 iterations of the Poincaré maps.


### 3.1. Linear JJ Case

In the case of a linear JJ withf(y)=f1(y), the bifurcation structure of the system (3) is as shown in Figure 3. This figure is reflectively symmetric about the line B0=cnπ=0.4nπ because of the symmetry with respect to the operation (13) discussed in Section 2.2. Consequently, near B=0.4π, there arise many bifurcations of the solutions which are symmetric to each other with respect to the linear operation (13). Let us examine the bifurcations arising with B0=0.5π; see the 1-dimensional bifurcation diagram in Figure 4 and the trajectories in Figure 5. Increasing B from B=0.2 as shown in Figure 4(a), the 1-periodic orbit (Figure 5(a)) becomes unstable near B=0.6 and a 2-periodic orbit (Figure 5(b)) is generated at the same time; e.g., a period-doubling bifurcation I1 arises. This 2-periodic orbit also becomes unstable through another period-doubling bifurcation I2, and a 4-periodic orbit appears (Figure 5(c)). Repeating these processes, there finally arises a chaotic attractor having the maximum Lyapunov exponent of 0.32 near B=0.64 (Figure 5(d)). Upon increasing B further, there arises a 3-periodic orbit near B=1.4 (Figure 5(e)). This solution has a large amplitude compared with the solutions with small B. It becomes unstable by undergoing a period-doubling bifurcation near B=1.45 and a 6-periodic orbit appears at the same time, as shown in Figure 5(f). At B=1.9, a chaotic attractor having the maximum Lyapunov exponent of 0.14 emerges (Figure 5(g)). From the 1-dimensional bifurcation diagram, this chaotic attractor visits three different regions one after another. However, at B=2.0, the chaotic attractor starts to wander in one large region. We consider that this difference comes from the global bifurcation. After that, the chaotic attractor disappears as B increases and a 3-periodic orbit appears again near B=2.5. This orbit disappears because of the tangent bifurcation, and then the solution converges to a 1-periodic orbit, as shown in Figure 5(h). On the other hand, by decreasing B from B=3.2, we observe a little different response from the increasing case (Figure 4(b)). Near B=2.5, in the region (A’), there is a 1-periodic orbit instead of the 3-periodic orbit observed in region (A) in the case of increasing B. This means that two periodic orbits coexist at B=2.5. Upon decreasing B from here, the 1-periodic orbit undergoes a tangent bifurcation G1 and suddenly disappears. After the disappearance, the solution does not converge to the 3-periodic orbit but behaves chaotically with the maximum Lyapunov exponent of 0.43 (Figure 5(i)).Figure 3
Bifurcation diagram of system (3) in B-B0 plane with f(y)=f1(y) and parameters (20), B∈(0.2,3.2), and B0∈(0,0.8π).Figure 4
1-dimensional bifurcation diagram of system (3) with f(y)=f1(y), B0=0.5π, and parameters (20). y∗ corresponds to the y-coordinate of the Poincaré map.
(a)
By increasingB from B=0.2. (b)
By decreasingB from B=3.2.Figure 5
Trajectories of system (3) with f(y)=f1(y). The parameters and initial conditions are in (20), B0=0.5π, and (a) B=0.50, x0=(3.966,0.681), (b) B=0.60, x0=(3.316,0.432), (c) B=0.62, x0=(3.168,0.346), (d) B=0.64, x0=(3.597,0.704), (e) B=1.40, x0=(-0.018,-0.207, (f) B=1.45, x0=(0.114,2.863), (g) B=1.90, x0=(-1.311,3.594), (h) B=2.50, x0=(-0.679,1.895), and (i) B=2.35, x0=(-0.122,1.778).
(a) (b) (c) (d) (e) (f) (g) (h) (i)By the way, there is a discontinuity in the curve ofG3 near B=1.7 and B0=0.5π. We consider that this is because a cusp structure might be here; however, we have not confirmed this hypothesis yet.


### 3.2. PWL JJ Case

In the case of a PWL JJ withf(y)=f2(y), the bifurcation structure of the system (3) is as shown in Figure 6. The same as in the linear JJ case, the bifurcation diagram is symmetric about B0=0.4π. Comparing Figure 6 with Figure 3, it is clear that the precise shapes of some of the bifurcation sets differ from each other. For example, the shaded region enclosed by I2 in Figure 3 looks bent; however, it has a rectangular form in Figure 6. On the other hand, the 1-dimensional bifurcation diagram in Figure 7 reveals the strong similarity of the solutions in these two cases. Comparing Figure 7 with Figure 4, no remarkable differences appear in the responses. Consequently, we have found that the linear conductance gives almost the same bifurcation structure as the PWL conductance, which contains more realistic characteristics.Figure 6
Bifurcation diagram of system (3) in B-B0 plane with f(y)=f2(y) and parameters (20), B∈(0.2,3.2), and B0∈(0,0.8π).Figure 7
1-dimensional bifurcation diagram of system (3) with f(y)=f2(y), B0=0.5π, and parameters (20). y∗ corresponds to the y-coordinate of the Poincaré map.
(a)
By increasingB from B=0.2. (b)
By decreasingB from B=3.2.Another phenomenon in Figure6 is that the local bifurcation sets break at the intersection with grazing bifurcation sets, which are curves in the figure. This is because the solution becomes singular when it is tangent to the border, and as a result, the calculation to obtain the local bifurcation sets stops.Focusing on these grazing bifurcations, let us discuss the change in stability of periodic orbits. Figure8(a) shows an enlargement of Figure 6. Since there are too many bifurcation sets in this figure, we will divide them into two groups based on the period of the periodic orbits, as shown in Figures 8(b) and 8(c). The bifurcation sets in Figure 8(b) relate to 1-periodic and 2-periodic orbits; the bifurcation sets in Figure 8(c) relate to 3-periodic orbits. From Figure 8(b), we know that a 2-periodic orbit exists at (B0,B)=(0.5π,1.2). We can also see this orbit in both 1-dimensional bifurcation diagrams in Figure 7. This solution undergoes a period-doubling bifurcation I2 as B increases, and it soon becomes chaotic, at B=1.25 in Figure 7. In other words, a chaotic set exists in the blue shaded region. Notice that this chaotic set might or might not be attractive. On the other hand, from Figure 8(c), a 3-periodic orbit exists in the shaded region enclosed by G3. For example, we can confirm a stable 3-periodic orbit with the parameter P1, as shown in the top figure of Figure 9. This solution undergoes another period-doubling bifurcation I3 as B0 decreases and becomes unstable. After this bifurcation, i.e., with the parameter P2 in the blue shaded region enclosed by I3, a stable 6-periodic orbit appears (Figure 6). As B increases from P1, the 3-periodic orbit undergoes a grazing bifurcation Z3. When the grazing bifurcation of a periodic orbit arises, the orbit suddenly disappears and another periodic orbit appears at the same time, which is shown in Figure 6. Each of them has exactly different stability although they are very similar in shape. We see this by observing the change in the characteristic multiplier of the 3-periodic orbit. The solution orbits with B=1.61455 and B=1.61456 have different characteristic multipliers from each other, e.g., (μ1,μ2)=(-0.883093049868,-0.037654415333) versus (μ1,μ2)=(-0.670464915282,-0.037623174381), because the grazing bifurcation Z3 arises between them. From the above, we have found that there are some nonlinear phenomena including the effect of the nonlinearity of the PWL function near B0=0.4π, including a grazing bifurcation, termination of the local bifurcation sets, and discontinuous jump of the stability of the periodic orbit. These phenomena do not appear in the case of the linear conductance; in other words, they are phenomena that are typical to the case of the PWL (realistic) conductance.Figure 8
Enlargement of Figure6. (a) All bifurcation sets that we calculated; (b) bifurcation sets of 1-periodic and 2-periodic orbits; (c) bifurcation sets of 3-periodic orbits.
(a) (b) (c)Figure 9
Trajectories of system (3) with f(y)=f2(y). The parameters and initial conditions are in (20) and P1:(B0,B)=(1.45,1.60), x0=(-0.522,3.126), P2:(B0,B)=(1.42,1.60), x0=(-0.606,3.067), and P3:(B0,B)=(1.45,1.63), x0=(-0.596,3.146).


## 3.1. Linear JJ Case

In the case of a linear JJ withf(y)=f1(y), the bifurcation structure of the system (3) is as shown in Figure 3. This figure is reflectively symmetric about the line B0=cnπ=0.4nπ because of the symmetry with respect to the operation (13) discussed in Section 2.2. Consequently, near B=0.4π, there arise many bifurcations of the solutions which are symmetric to each other with respect to the linear operation (13). Let us examine the bifurcations arising with B0=0.5π; see the 1-dimensional bifurcation diagram in Figure 4 and the trajectories in Figure 5. Increasing B from B=0.2 as shown in Figure 4(a), the 1-periodic orbit (Figure 5(a)) becomes unstable near B=0.6 and a 2-periodic orbit (Figure 5(b)) is generated at the same time; e.g., a period-doubling bifurcation I1 arises. This 2-periodic orbit also becomes unstable through another period-doubling bifurcation I2, and a 4-periodic orbit appears (Figure 5(c)). Repeating these processes, there finally arises a chaotic attractor having the maximum Lyapunov exponent of 0.32 near B=0.64 (Figure 5(d)). Upon increasing B further, there arises a 3-periodic orbit near B=1.4 (Figure 5(e)). This solution has a large amplitude compared with the solutions with small B. It becomes unstable by undergoing a period-doubling bifurcation near B=1.45 and a 6-periodic orbit appears at the same time, as shown in Figure 5(f). At B=1.9, a chaotic attractor having the maximum Lyapunov exponent of 0.14 emerges (Figure 5(g)). From the 1-dimensional bifurcation diagram, this chaotic attractor visits three different regions one after another. However, at B=2.0, the chaotic attractor starts to wander in one large region. We consider that this difference comes from the global bifurcation. After that, the chaotic attractor disappears as B increases and a 3-periodic orbit appears again near B=2.5. This orbit disappears because of the tangent bifurcation, and then the solution converges to a 1-periodic orbit, as shown in Figure 5(h). On the other hand, by decreasing B from B=3.2, we observe a little different response from the increasing case (Figure 4(b)). Near B=2.5, in the region (A’), there is a 1-periodic orbit instead of the 3-periodic orbit observed in region (A) in the case of increasing B. This means that two periodic orbits coexist at B=2.5. Upon decreasing B from here, the 1-periodic orbit undergoes a tangent bifurcation G1 and suddenly disappears. After the disappearance, the solution does not converge to the 3-periodic orbit but behaves chaotically with the maximum Lyapunov exponent of 0.43 (Figure 5(i)).Figure 3
Bifurcation diagram of system (3) in B-B0 plane with f(y)=f1(y) and parameters (20), B∈(0.2,3.2), and B0∈(0,0.8π).Figure 4
1-dimensional bifurcation diagram of system (3) with f(y)=f1(y), B0=0.5π, and parameters (20). y∗ corresponds to the y-coordinate of the Poincaré map.
(a)
By increasingB from B=0.2. (b)
By decreasingB from B=3.2.Figure 5
Trajectories of system (3) with f(y)=f1(y). The parameters and initial conditions are in (20), B0=0.5π, and (a) B=0.50, x0=(3.966,0.681), (b) B=0.60, x0=(3.316,0.432), (c) B=0.62, x0=(3.168,0.346), (d) B=0.64, x0=(3.597,0.704), (e) B=1.40, x0=(-0.018,-0.207, (f) B=1.45, x0=(0.114,2.863), (g) B=1.90, x0=(-1.311,3.594), (h) B=2.50, x0=(-0.679,1.895), and (i) B=2.35, x0=(-0.122,1.778).
(a) (b) (c) (d) (e) (f) (g) (h) (i)By the way, there is a discontinuity in the curve ofG3 near B=1.7 and B0=0.5π. We consider that this is because a cusp structure might be here; however, we have not confirmed this hypothesis yet.


## 3.2. PWL JJ Case

In the case of a PWL JJ withf(y)=f2(y), the bifurcation structure of the system (3) is as shown in Figure 6. The same as in the linear JJ case, the bifurcation diagram is symmetric about B0=0.4π. Comparing Figure 6 with Figure 3, it is clear that the precise shapes of some of the bifurcation sets differ from each other. For example, the shaded region enclosed by I2 in Figure 3 looks bent; however, it has a rectangular form in Figure 6. On the other hand, the 1-dimensional bifurcation diagram in Figure 7 reveals the strong similarity of the solutions in these two cases. Comparing Figure 7 with Figure 4, no remarkable differences appear in the responses. Consequently, we have found that the linear conductance gives almost the same bifurcation structure as the PWL conductance, which contains more realistic characteristics.Figure 6
Bifurcation diagram of system (3) in B-B0 plane with f(y)=f2(y) and parameters (20), B∈(0.2,3.2), and B0∈(0,0.8π).Figure 7
1-dimensional bifurcation diagram of system (3) with f(y)=f2(y), B0=0.5π, and parameters (20). y∗ corresponds to the y-coordinate of the Poincaré map.
(a)
By increasingB from B=0.2. (b)
By decreasingB from B=3.2.Another phenomenon in Figure6 is that the local bifurcation sets break at the intersection with grazing bifurcation sets, which are curves in the figure. This is because the solution becomes singular when it is tangent to the border, and as a result, the calculation to obtain the local bifurcation sets stops.Focusing on these grazing bifurcations, let us discuss the change in stability of periodic orbits. Figure8(a) shows an enlargement of Figure 6. Since there are too many bifurcation sets in this figure, we will divide them into two groups based on the period of the periodic orbits, as shown in Figures 8(b) and 8(c). The bifurcation sets in Figure 8(b) relate to 1-periodic and 2-periodic orbits; the bifurcation sets in Figure 8(c) relate to 3-periodic orbits. From Figure 8(b), we know that a 2-periodic orbit exists at (B0,B)=(0.5π,1.2). We can also see this orbit in both 1-dimensional bifurcation diagrams in Figure 7. This solution undergoes a period-doubling bifurcation I2 as B increases, and it soon becomes chaotic, at B=1.25 in Figure 7. In other words, a chaotic set exists in the blue shaded region. Notice that this chaotic set might or might not be attractive. On the other hand, from Figure 8(c), a 3-periodic orbit exists in the shaded region enclosed by G3. For example, we can confirm a stable 3-periodic orbit with the parameter P1, as shown in the top figure of Figure 9. This solution undergoes another period-doubling bifurcation I3 as B0 decreases and becomes unstable. After this bifurcation, i.e., with the parameter P2 in the blue shaded region enclosed by I3, a stable 6-periodic orbit appears (Figure 6). As B increases from P1, the 3-periodic orbit undergoes a grazing bifurcation Z3. When the grazing bifurcation of a periodic orbit arises, the orbit suddenly disappears and another periodic orbit appears at the same time, which is shown in Figure 6. Each of them has exactly different stability although they are very similar in shape. We see this by observing the change in the characteristic multiplier of the 3-periodic orbit. The solution orbits with B=1.61455 and B=1.61456 have different characteristic multipliers from each other, e.g., (μ1,μ2)=(-0.883093049868,-0.037654415333) versus (μ1,μ2)=(-0.670464915282,-0.037623174381), because the grazing bifurcation Z3 arises between them. From the above, we have found that there are some nonlinear phenomena including the effect of the nonlinearity of the PWL function near B0=0.4π, including a grazing bifurcation, termination of the local bifurcation sets, and discontinuous jump of the stability of the periodic orbit. These phenomena do not appear in the case of the linear conductance; in other words, they are phenomena that are typical to the case of the PWL (realistic) conductance.Figure 8
Enlargement of Figure6. (a) All bifurcation sets that we calculated; (b) bifurcation sets of 1-periodic and 2-periodic orbits; (c) bifurcation sets of 3-periodic orbits.
(a) (b) (c)Figure 9
Trajectories of system (3) with f(y)=f2(y). The parameters and initial conditions are in (20) and P1:(B0,B)=(1.45,1.60), x0=(-0.522,3.126), P2:(B0,B)=(1.42,1.60), x0=(-0.606,3.067), and P3:(B0,B)=(1.45,1.63), x0=(-0.596,3.146).


## 4. Conclusions

We conducted a bifurcation analysis of a single-junction superconducting quantum interferometer with an external flux. We used two types of conductance: a linear conductance and a PWL conductance, to approximate the current-voltage characteristic of the JJ. We described a method to compute the local stability of the solution orbit and to solve the bifurcation problem of the system. For the case of the PWL function, we applied the hybrid dynamical system approach to the system.The main results of the analysis of system (3) are listed below:(1)
We exactly revealed the bifurcation structure of the system in two-dimensional parameter space. From the one-dimensional bifurcation diagram, we plainly explained what kind of bifurcation arises.(2)
The comparison of the linear function and PWL function cases indicated a difference in the shapes of the bifurcation sets in two-dimensional parameter space, despite that there were no differences in the one-dimensional bifurcation diagrams or trajectories.(3)
We discovered that grazing bifurcations terminate the calculation of the local bifurcations because they drastically change the stability of the periodic orbit.As future work, we would like to explore cases with other values ofθ1 and θ2, because we believe that they should have a large influence on the bifurcation structure. In addition, we should consider the global bifurcation problem.



---

*Source: 8931525-2018-12-02.xml*