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

We conduct a bifurcation analysis of a single-junction superconducting quantum interferometer with an external flux. We approximate the current-voltage characteristicsof 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.


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][6][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 singlejunction 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).

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 ext as the sum of DC and AC components [16], we get the equivalent circuit shown in Figure 1(b). The circuit equations for the circuit in Figure 1 are obtained by using Kirchhoff 's laws and the characteristics of each element: where is the junction critical current, (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: In dimensionless form, assuming ( ) = 0 + cos , (2) becomes as follows: We can write the state variables and in vector notation: is a function of that corresponds to the conductance (V) in the JJ equivalent circuit. The response of is precisely modeled from experiments, as shown in Figure 2(a). In this study, we approximate 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 : For the PWL case, let us firstly define the state spaces as follows: The PWL function in Figure 2(b) is as follows: where for = 1, 2 corresponds to thresholds dividing ( ) into five segments. In this study, we will avoid dealing with the hysteresis phenomenon arising in the JJ.
An orbit ( ) ∈ 2 is the solution of (3) together with an initial condition 0 = ( 0 , 0 ) ∈ 2 . We define ( ) to be a function of time and 0 : 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 and , namely, the time when = ± . An orbit ( ) is periodic if The period of this orbit is the minimum value of . Since one of the state variables in (3) corresponds to an angle, the motion of this system evolves on a cylindrical surface. Let 1 be the set of points on a unit circle in 2 : where gives the angle of a point in 1 from the reference direction; accordingly, we define the cylindrical surface, Given a map where is scaled and biased appropriately to fit the polar coordinate system, and the cylindrical surface is homeomorphic to 2 \{0}. We will thus observe the trajectory in the punctured plane 2 \{0}. (3) is invariant to the following linear transformation:

Symmetry. Equation
where is an arbitrary integer. This means that we can transform a solution with 0 = * 0 and a solution with 0 = 2 − * 0 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 2 − 0 → 0 ; that is, the bifurcation sets are reflectively symmetric about the line 0 = .

Local Stability and Bifurcations.
When discussing the asymptotic stability of the periodic orbit, we generally use the Poincaré map [20] of the orbit, which is a discrete map as follows: In the field of discrete-time dynamical systems, the point ∈ 2 is called a fixed point of the map if Similarly, the point is called an ℓ-periodic point of the map if These periodic points of exactly correspond to the periodic orbit of the system (3). The variational equation from the fixed point with respect to the Poincaré map is as follows: where is the difference between and = ( + 0 ). The multipliers 1 and 2 are obtained by solving the characteristic equation: where is a 2×2 identity matrix. When ∀ , ̸ = 0 and | | ̸ = 1, is called a 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 is completely stable if its multipliers satisfy | 1 | < 1 and | 2 | < 1; accordingly, we label the point 0 .
The index at the bottom left of the symbol indicates the number of unstable dimensions for the fixed point. Generally, 0 is called a node, 1 and 1 are called saddles, and 2 is called a source. Similarly to the fixed point, we classify the stability of the ℓ-periodic point by using symbols, i.e., 0 ℓ , etc.
We numerically obtain the point 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 of changes. This change of the stability is called a local bifurcation of [21], and each of the parameter sets is called a local bifurcation set. Local bifurcations occur when one of the multipliers of exceeds unity as a result of changing a parameter. In the system (3), two local bifurcations are possible: the tangent bifurcation ( = 1) and the period-doubling bifurcation ( = −1). We numerically obtain the local bifurcation sets by using the method in [12,19].

Grazing Bifurcation.
Let us consider the case that a periodic orbit passes near the border = . 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 called a grazing bifurcation of the periodic orbit [13], and a parameter set where a grazing bifurcation phenomenon arises is called a 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 : where ± ( ) = ± = 0 corresponds the border of the system. We compute grazing bifurcation sets by simultaneously solving (15) and (19) with Newton's method.
We will observe at most three equilibrium points if we choose = 0.4, as mentioned in [22]. (Parameter in [22] corresponds to .) We further assume that the parameters 0 and 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: , which is plotted as a magenta solid curve, represents the tangent bifurcation set of an -periodic point; , which is plotted as a blue solid curve, represents the period-doubling bifurcation set of anperiodic point; 1 , which is plotted as a black dashed curve, represents the grazing bifurcation set of an -periodic orbit. We plot the solution orbits in the punctured plane 2 \{0}. Here, we assume that is a subset of 2 \{0} and corresponds to for = 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 use the 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.

Linear JJ Case.
In the case of a linear JJ with ( ) = 1 ( ), the bifurcation structure of the system (3) is as shown in  Figure 4 and the trajectories in Figure 5. Increasing from = 0.2 as shown in Figure 4(a), the 1-periodic orbit ( Figure 5(a)) becomes unstable near = 0.6 and a 2-periodic orbit ( Figure 5(b)) is generated at the same time; e.g., a period-doubling bifurcation 1 arises. This 2-periodic orbit also becomes unstable through another period-doubling bifurcation 2 , 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 = 0.64 ( Figure 5(d)). Upon increasing further, there arises a 3-periodic orbit near = 1.4 ( Figure 5(e)). This solution has a large amplitude compared with the solutions with small . It becomes unstable by undergoing a period-doubling bifurcation near = 1.45 and a 6-periodic orbit appears at   the same time, as shown in Figure 5(f). At = 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 = 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 increases and a 3-periodic orbit appears again near = 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 from = 3.2, we observe a little different response from the increasing case ( Figure 4(b)). Near = 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 . This means that two periodic orbits coexist at = 2.5. Upon decreasing from here, the 1-periodic orbit undergoes a tangent bifurcation 1 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)).
By the way, there is a discontinuity in the curve of 3 near = 1.7 and 0 = 0.5 . We consider that this is because a cusp structure might be here; however, we have not confirmed this hypothesis yet.

PWL JJ Case.
In the case of a PWL JJ with ( ) = 2 ( ), 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 0 = 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 2 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. Another phenomenon in Figure 6 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. Figure 8(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 ( 0 , ) = (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 2 as increases, and it soon becomes chaotic, at = 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 3 . For example, we can confirm a stable 3-periodic orbit with the parameter 1 , as shown in the top figure of Figure 9. This solution undergoes another period-doubling bifurcation 3 as 0 decreases and becomes unstable. After this bifurcation, i.e., with the parameter 2 in the blue shaded region enclosed by 3 , a stable 6-periodic orbit appears ( Figure 6). As increases from 1 , the 3-periodic orbit undergoes a grazing bifurcation 3 . 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.    Figure 8: Enlargement of Figure 6. (a) All bifurcation sets that we calculated; (b) bifurcation sets of 1-periodic and 2-periodic orbits; (c) bifurcation sets of 3-periodic orbits. the grazing bifurcation 3 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 0 = 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.

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,

10
Complexity despite that there were no differences in the onedimensional 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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.