Semiclassical dynamics of a superconducting circuit: chaotic dynamics and fractal attractors

We study here the semiclassical dynamics of a superconducting circuit constituted by two Josephson junctions in series, in the presence of a voltage bias. We derive the equations of motion for the circuit through a Hamiltonian description of the problem, considering the voltage sources as semi-holonomic constraints. We find that the dynamics of the system corresponds to that of a planar rotor with an oscillating pivot. We show that the system exhibits a rich dynamical behaviour with chaotic properties and we present a topological classification of the cyclic solutions, providing insight into the fractal nature of the dynamical attractors.


Introduction
Superconducting circuits are arguably one of the most successful quantum computing platforms [1].At the heart of these circuits is the Josephson junction (JJ), which, in a lumped-elements description, acts as a nonlinear circuit element [2,3].In a quantum setting, this nonlinearity allows isolating a two-level system within the full spectrum, and, consequently, constructing a wide collection of superconducting qubits [4][5][6], depending on the specific topology of the circuit.Beyond the field of quantum computation, the physics of JJs has also been discussed and applied in other contexts, such as the study of quantum phase transitions in condensed matter physics [7,8], sensing and metrology applications [9], where striking consequences of the tunnelling between superconductors can already be identified at the semiclassical level.One example is given by the Aharonov-Bohm effect [10] in SQUIDS.These devices, consisting of a superconducting loop interrupted by two dielectric layers and modelled as two JJs in parallel, exploit the relation between the superconducting phase and the magnetic flux to act as extremely sensitive magnetometers.
In this article, we focus on this semiclassical aspect: in a circuit somewhat analogous to the SQUID, we describe the dynamics of the superconducting phase of a device constituted, unlike the SQUID, by two JJs in series in the presence of an external voltage bias V g (see figure 1a).This setup defines a superconducting island coupled to the rest of the circuit through two (superconducting) tunnel junctions.As we will discuss in greater detail below, we note that, even if this circuit has the same topology of a superconducting single-electron transistor (SSET) at zero gate voltage [11], our analysis is valid in a regime which is complementary to the one required to observe charge quantization in SSETs.
The dynamics of the circuit is discussed in terms of Hamiltonian dynamics, allowing us to recognize how the equation of motion for the phase on the island between the JJs (see figure 1) can be mapped onto the dynamics of an equivalent mechanical system: a parametrically driven planar rotor (a planar rotor whose pivot is periodically driven).While the correspondence between the dynamics of a single JJ and a driven physical pendulum has been extensively studied in the literature [12][13][14][15][16][17], the relation between the series of two JJs and a mechanical system is, to our knowledge, novel.
The interest in the parametrically driven planar rotor lies in the fact that it can be considered as a zero-gravity parametrically driven pendulum (PDP).The latter is a system with a rich dynamical behaviour, exhibiting a wide collection of attractors and fixed points, such as the non-inverted (θ = 0) and inverted (θ = π) positions.In particular, the inverted position can be stabilized for certain values of the system parameters [18][19][20][21] -a configuration often referred to as the Kapitza pendulum [22].The general stability properties of the PDP have been studied in several works [20][21][22][23][24][25][26][27][28][29], considering gravity, driving amplitude and, in some instances, friction as free parameters.While the analytical investigation has mainly focused on the stability of the inverted position in terms of solutions of the Mathieu equation [26], numerical work has addressed various aspects of the PDP dynamics.These include the appearance of nontrivial limit-cycles (multiple nodding solutions [23,28]), the emergence of fractal basins of attraction [25], and the investigation of "rotating solutions" (i.e.unbounded rotations) [27].
In our work, we first discuss the Hamiltonian formalism used to derive the equations of motion of the electrical degrees of freedom of the circuit.We start from the Lagrangian associated with the lumped-elements description of the device, in which the external voltage source is considered as a semi-holonomic constraint on the local voltages.This allows us to identify the dynamics of the circuit of figure 1a with that of the parametrically driven planar rotor in figure 1b.Having in mind our circuit implementation, we then discuss several new aspects of the PDP.First, we distinguish between unstable, 0-stable, π-stable and limit-cycle asymptotic attractors for the superconducting phase of the circuit φ ∆ (t).By focusing on the natural conditions for the JJ setting ( φ∆ (0) = 0), we provide a description of the fractal nature of the stability diagram obtained by varying the dimensionless driving amplitude and the initial superconducting phase: we compute the Hausdorff dimension of the attractors.Moreover, we classify the limit-cycle solutions (n-cycles) in terms of the number of intersections with the zero-velocity axis in the phase space, providing then a topological classification of these trajectories.We find that the system displays a wide collection of n-cycles characterized by a chaotic distribution in the stability diagram, identifying cyclic solutions having up to n = 18.Furthermore, we analyze the case φ∆ (0) ̸ = 0 for two fixed values of the dimensionless driving amplitude, obtaining basins of attraction with fractal geometry, whose borders have non-integer Hausdorff dimension.

Semiclassical model and equations of motion
The system considered is a superconducting circuit constituted by two JJs in series, in the presence of a voltage bias V g (figure 1a).In a lumped-elements circuit description, the two Josephson junctions are characterized by capacitances C J1 and C J2 and Josephson energies E J1 and E J2 , respectively.While we initially focus on a purely reactive circuit model, we later introduce dissipation by considering a shunt resistance for each of the two JJs, using the RCSJ model [30].This resistive shunt can, for instance, describe a quasiparticle contribution to the tunnelling current.
Figure 1: (a) Equivalent circuit of the superconducting device.The Josephson junction consists of capacitances C J1 , C J2 and non-linear inductances E J1 , E J2 in parallel.We introduce shunting resistors R 1 ,R 2 to take into account dissipation.The system is voltage biased (V g ).(b) Schematic representation of a parametrically driven pendulum, consisting of a rod of length l with a pivot P oscillating in the vertical direction and a mass m attached to the free end.
The circuit we consider here has the same topology of a SSET [11].The key distinction between our setting and a SSET is that, for the latter, the focus is on the so-called Coulomb blockade regime [11], corresponding to charge quantization on the island.Such a regime is reached for a characteristic impedance larger than the resistance quantum (Z ≫ R K = h/e 2 ).Here, we focus on the opposite situation (Z ≪ R K ), in analogy to the strategy employed in the design of phase qubits [1,5,31].This can be achieved by allowing the tunnelling junctions to be large, even macroscopic, increasing their capacitance and thus lowering the charging energy E C ∝ 1/C or, alternatively, adding a large shunting capacitance: incidentally, we note that these conditions suggest that the experimental realization of this device should not pose any particular challenge.
A consequence of this choice of parameters is that the number of Cooper pairs on the island is no longer a good quantum number -contrary to the SET case-and, instead, the phase dynamics becomes the dominating effect, allowing us to treat the superconducting phase across the junctions semiclassically.
The dynamics of the electrical degrees of freedom is obtained through a Hamiltonian mechanics description of our circuit, starting from a lumped-elements model.In this approach, node fluxes and charges (time integrals of voltages and currents, respectively) are the conjugate variables that play the same role as position and momentum in a conventional classical mechanics problem.Our approach deviates from the standard method of describing the lumped-element model of a superconducting circuit [3,32] in that we directly include the voltage constraint imposed by the voltage source through the undetermined Lagrange multiplier method.To our knowledge, this approach has not been considered for the analysis of (superconducting) circuits, where, typically, the voltage source is replaced by capacitors which, in the limit of infinite capacitance, induce the correct node voltage provided by the voltage sources [3,32].We note that our analysis allows for a straightforward quantization procedure which, however, is not the focus of our work.
The general approach consists in identifying the nodes in the circuit under consideration (i ∈ {g, I, o}, in our case) and introducing flux variables ϕ i at each of them.The flux ϕ i is defined as ) where V i is the voltage at each node, implying that V i (t) = φi (t).The voltage source imposes the constraint φg − φo − V g = 0 on the circuit.Without loss of generality, we choose here the node o to be grounded, i.e.V o = 0 and hence φo = 0.The Lagrangian of the superconducting circuit can therefore be written as where is the inductive energy associated to the two JJs.The rightmost term in the Lagrangian consists of the undetermined Lagrange (charge) multiplier Q and the semi-holonomic (function of the flux derivatives) constraint φg − φo − V g = 0 on the voltages.Even if the constraint is non-holonomic, its linearity in the derivatives ⃗ φ ensures that the Lagrange multiplier method can be applied (see Appendix A).
The node charges Q i = ∂L/∂ φi , representing the momenta conjugated to the fluxes ϕ i , allow us to write the Hamiltonian where the charges are now written in terms of the number of Cooper pairs and the Lagrange multiplier as N = Q/(2e); furthermore, the phases φ i = 2πϕ i /Φ 0 are introduced (flux-phase relation).From equation ( 2), we can write Hamilton's equations for the node variables as Using the voltage constraint together with equation (3d), one can solve for the Lagrange multiplier realizing that it corresponds to the number of Cooper pairs provided by the external voltage source to the two JJ capacitances C J1 and C J2 in series, in absence of tunnelling across the junctions.
In the following, we will assume that Differentiating equation (3c) with respect to time and substituting equations (3a) and ( 4) in it, we obtain where φ ∆ = φ I − φ g /2, φ g (t) = 2πV g /Φ 0 t (we impose φ g (0) = 0) and Ω J = √ 8E C E J /ℏ is the Josephson frequency.Equation ( 5) can be written in a dimensionless form as where ε = 8Ω 2 J /ω 2 g , τ = ω g t/2 and ω g = 2πV g /Φ 0 .See Appendix A and Appendix B for further details on the derivation of the system's Hamiltonian and the equations of motion.Introducing dissipation, equation ( 6) becomes where the (dimensionless) dissipation term is given by κ = 4/ [ω g R(C J1 + C J2 )] and R describes the resistance used to model the losses within the RCSJ model: see Appendix D for a detailed derivation.
In our analysis, we will focus on the dynamics generated by equation ( 7) for different values of the parameter ε and the initial conditions [φ ∆ (0), φ∆ (0)].In terms of our device, this choice of the dynamical parameters appears natural: the parameter ε is controlled by the voltage V g ; the initial condition on the phase φ ∆ (0) corresponds to introducing a time-independent magnetic flux Φ B /Φ 0 = 2φ ∆ (0) for t < 0 through the main loop of the circuit.The initial condition of the phase derivative φ∆ (0) = 0 indicates a symmetric voltage drop across the JJs, i.e. φI (0) = V g /2.An initial condition φ∆ (0) ̸ = 0 corresponds therefore to an initially asymmetric voltage drop across the junction V 1,2 = V g /2 ∓ φ∆ (0), where V 1 (V 2 ) is the voltage drop across the first (second) JJ.
It is possible to recognize that equation ( 7) is a specific case of the equation controlling the PDP dynamics.The latter consists of a rigid planar rotor of length l whose pivot is driven harmonically along the vertical direction as Following [25], the equation of motion for the angle can be written as where β = b/l, α = g/(lω 2 ) and τ = ωt, and γ is a coefficient modelling dissipation.Setting β = ε and α = [g/(lω 2 )] = 0, we can identify the equation controlling the dynamics of our system with the PDP equations of motion for g = 0.
Interestingly, in the limit V g = 0 (ε → ∞), equation ( 5) reduces to the case of a physical (i.e.not parametrically driven) pendulum in the presence of gravity, whose value is determined by the Josephson frequency and the initial superconducting phase.From the PDP perspective, the V g = 0 condition in the JJ circuit translates into a constant acceleration of the pivot and therefore, into an apparent force applied to the oscillating mass, caused by the non-inertiality of the reference frame.The V g = 0 limit is not investigated in this work any further.

Numerical results
In this section, we perform a numerical analysis of equation (7) with two different settings of the parameters.In one case, we vary the dimensionless driving amplitude ε and the initial phase φ ∆ (0), fixing φ∆ (0) = 0.In the other case, we pick up two values of ε and we vary both the initial conditions [φ ∆ (0), φ∆ (t) ̸ = 0].In both the scenarios, we fix κ = 10 −2 .We recall that the parameter ε controls the external voltage, whereas φ∆ (0) tunes the initial voltage of the island.
We numerically determine the asymptotic behaviour of the solutions φ ∆ (t), identifying four general types of attractors as a function of the dimensionless driving amplitude ε and initial conditions [φ ∆ (0), φ∆ (0)]: 0-stable, π-stable, unstable, and limitcycle solutions.The appearance of four different types of attractors was previously discussed in the literature [23,25,28], where the rotor dynamics is analyzed as a function of the initial conditions and the system parameters: driving amplitude, gravity and, in some cases, dissipation.In our analysis, the solutions 0-stable and π-stable are such that φ ∆ = 0 and φ ∆ = π are stable fixed points (see figures 2a and 2c).The latter is reminiscent of the stabilization of the inverted position in the Kapitza pendulum  in the presence of gravity.Unstable solutions correspond to asymptotically rotating solutions [27], which are not characterized here any further (see figure 2b).The limit-cycle solutions are periodic orbits around 0 and π that exhibit a varying degree of complexity, depending on the driving amplitude and the initial conditions.In our work, we characterize these trajectories in terms of their number of intersections with the φ∆ = 0 axis in the phase space (see figures 2d and 3): if a periodic trajectory crosses the φ∆ = 0 axis 2n times in a period, it is classified as an n-cycle (see Appendix D).It is straightforward to realize that the number of intersections n coincides with the number of loops, which is characterized by the turning number t.In addition, for a closed trajectory, it is possible to define the winding number w, counting the number  of rotations around a given point ⃗ O. Choosing in our case ⃗ O = (π, 0) or (0, 0), we find both n-cycles for which w = t (figure 3d-3e-3f) and cycles for which w < t (figure 3a-3b-3c): trajectories execute loops on one side of the phase space before passing to the The diagram is divided into four regions by the main attractors, each labelled with a different colour: black stands for unstable, dark blue for 0-stable, blue for π-stable and the remaining colours represent the different n-cycles, which combine to form the limit-cycles region.The most likely are 1-cycle (light blue), 3-cycle (green) and 5-cycle (light green), divided into three different subregions.However, the limit-cycle regions exhibit stripes and domains exhibiting higher-order n-cycles.(b) Inset of (a), where a region of 2-cycles (cyan) appears along line A, standing for ϵ = 0.509253.(c) Inset of (a), exhibiting subregions of 6-cycles (yellow), 9-cycles (orange) and a narrow stripe of 15-cycles (red) for ε = 0.373 .Line B is for ε = 0.361176755.other side.We note that the same n-cycles can occur both with t = w and t > w.
In the literature, limit-cycle solutions have been classified in different ways (see [23] and [33]), but, to our knowledge, their topological properties have not been considered.In [23], these trajectories are described in terms of their multiple nodding behaviour observed in terms of the driving amplitude, dissipation and gravity, where a nodding corresponds to an instance where the pendulum has reached its local extremum away from the vertical position and starts to move back towards it: n A refers to the number of noddings counted in one-half of a limit-cycle.Comparing the limit-cycles analyzed in [23] with our topological description, we find that in [23] only multiple nodding trajectories with w = 1 are considered and that n = 2n A −w.In [33], where zero-gravity case is considered, the limit-cycle solutions are described in terms of the subharmonic resonance, a phenomenon occurring when the driving frequency of the pivot is an integer multiple of the frequency of the cyclic trajectory of the pendulum.However, a clear connection between a specific subharmonic and a trajectory with certain number of nods is not drawn.
We have extensively analyzed the case in which ε and φ ∆ (0) are varied and we fix φ∆ (0) = 0, i.e. the JJs share the same initial voltage drop V g /2.The results are summarized in the stability diagram of figure 4, where the distribution of the attractors and n-cycles in the parameters space is reported.The main attractors -unstable, 0-stable, π-stable and limit-cycles (composed of all the n-cycles together)-form four different regions in the stability diagram, showing a fractal geometry and high sensitivity to initial conditions (see figure 4): the same behaviour is also encountered for the ncycles subregions.From figure 4a it is possible to see that 1-cycle, 3-cycle and 5-cycle solutions identify three disjoint subsets, still exhibiting intricate geometries.A closer examination of the 1-cycle (figure 4c) and 3-cycle (figure 4b) regions displays a richer scenario where higher-order n-cycles appear.For instance, in the 1-cycle region, a 2cycle stripe appears along ε ≃ 0.509253 (line A in figure 4c), along with (less prominent) 5-cycle and 3-cycle stripes for ε ≃ 0.5 and ε ≃ 0.53 (see figure 4c).In the 3-cycle region (see figure 4b), respectively around ε ≃ 0.385 and ε ≃ 0.377, 9-cycle and 6-cycle stripes appear, in addition to a (faint) 15-cycle stripe for ε ≃ 0.373 and a narrow line of 18-cycles below the 9-cycle stripe (for a detailed description of the attractors and the n-cycles analysis, see Appendix D).The attractors in the stability diagram exhibit an intricate distribution in terms of the parameters and high sensitivity to initial conditions.Our numerical investigation provides further insight into these properties of the stability diagram: we characterize the fractal dimension of each region in the stability diagram as a function of the two parameters ε and φ ∆ (0).We do this by estimating the Hausdorff dimension of their borders in figure 4a with the box-counting method [34]: the unstable solutions have a Hausdorff dimension δ H ≃ 1.520, 0-stable δ H ≃ 1.327, π-stable δ H ≃ 1.307 and limit-cycles (comprised of all the n-cycles) δ H ≃ 1.201.The list of different n-cycles observed in our simulations is reported in table 1, along with the fractal dimensions of the regions they form in the stability diagram in figure 4 (for a detailed discussion about the estimation of the Hausdorff dimensions, see Appendix D).   4b).The colour-coding is the same as figure 4.
The case φ∆ (0) ̸ = 0, i.e. asymmetric initial voltage drop V g /2 ∓ φ∆ (0) across the two JJs, has been analyzed for ε ≃ 0.509253 (line A in figure 4c) and ε ≃ 0.361176755 (line B in figure 4b) varying both the initial conditions [φ ∆ (0), φ∆ (0)].The attractors appearing from these parameter settings are of the same type as the case φ∆ (0) = 0: the results of the numerical analysis are shown in figure 5, where we report the basins of attraction related to the different attractors.In addition, we note that the number n of the intersection of the n-cycles seem not to change when we move away from the condition φ∆ (0) = 0, at least for the lines A and B: comparing figures 4c and 5a, line A exhibits 1-and 2-cycles both for φ∆ (0) = 0 and φ∆ (0) ̸ = 0 and the same happens for line B, which exhibits only 1-cycles throughout its dynamics (see figures 4b and 5b).
In [25], the authors show the basins of attraction of the rotor dynamics obtained as a function of the initial conditions and for fixed values of the driving amplitude and the gravity.Even though we use a different parameters space, our basins of attraction in figure 5 exhibit a fractal structure similar to the one reported in [25]: here, we emphasize this property by computing the Hausdorff dimension of the borders of the attractors, in the same way as the case φ∆ (0) = 0.For the case depicted in figure 5a the unstable solutions have dimension δ H = 1.233, δ H = 1.036 for 0-stable, δ H = 1.250 for 1-cycles and δ H = 1.123 for 2-cycles, whereas for figure 5b the Hausdorff dimensions are δ H = 1.434 for unstable solutions, δ H = 1.414 for 0-stable, δ H = 1.310 for π-stable and δ H = 1.246 for 3-cycle.

Conclusions
Through a classical mechanics approach that takes into account voltage sources as semi-holonomic constraints, we have shown that the semiclassical dynamics of a circuit comprising two Josephson junctions in series can be described by the same dynamical equations of the planar (zero-gravity) parametrically driven pendulum.In our numerical analysis, we identify four different types of attractors -0-and π-stable, unstable and limit-cycle-accessible by tuning the dynamical parameters of the system.Moreover, we classify the stable limit-cycle solutions as n-cycles, where n equals the turning number of the limit-cycle.In addition, we describe the distribution of the attractors in the stability diagram and basins of attraction, which reveals the chaotic behaviour of the system dynamics by characterizing the different attractors according to their Hausdorff dimension.

Appendix B. Deriving Hamilton's equations and the dynamical equation for superconductive phase
The canonical conjugate variables in our formalism are the flux ϕ i (t) and the charge Q i (t) at each node i = g, I.However, in this case, we write the Hamiltonian in terms of the number of Cooper pairs n i = Q i /(2e) and the phase φ i = 2πϕ i /Φ 0 .For this reason, taking f and g as two generic functions of the canonical variables, the usual Poisson brackets become ℏ The voltage constraint ω g = φg − φo = φg gives that the time evolution of the superconducting phase (related to the voltage) on the lead g given by equation (B.3d) is a constant ω g = 2πV g /Φ 0 and we can solve for the Lagrange multiplier The physical interpretation of the multiplier N can be given by considering a slightly different setting of the circuit, in which the tunnelling across junctions is neglected and the external voltage source is replaced by the capacitor Thus, the Lagrange multiplier N corresponds to the number of Cooper pairs on the capacitance C ∥ , which, considering the original circuit with the voltage source, can be regarded as the charge stored by the external voltage V g on the junctions capacitances C J1 and C J2 in series.Furthermore, we note that the solution of the Lagrange multiplier allows us to write the Hamiltonian in equation (A.8) in a more familiar form

Appendix C. Derivation of the dissipative equation of motion
Let us derive the dissipative equation of motion for the superconductive phase using the resistively and capacitively shunted junction model, where a Josephson junction is placed in parallel, as the name suggests, with a capacitor and a resistor to approximate loss mechanisms in the junction.To this end, we utilize the two Josephson relations where I C is the critical current, δ the superconducting phase difference across the junction, and V the voltage across the junction.For both JJs in our circuit, the current through them is Let us write the phase differences across the junctions using the explicit node phases so that δ 1 = φ g − φ I and δ 2 = φ I − φ o .The two JJs are in series, so the electric current through them must be equal We note that node 3 is grounded, i.e. φ o = φo = φo = 0, and due to the voltage source in the circuit, the phase of node 2 satisfies Assuming identical junctions with R 1 = R 2 = R and recalling that the critical current relates to the Josephson energy so that I C = 2π ϕ 0 E J , we can solve equation (C.4) to obtain φI = 1 (C.7) Due to the discreteness of the points in the stability diagram, Hausdorff dimension is estimated via linear fits (see figure D1), exploiting the formula where N b (r) is the number of boxes having radius r needed to cover the border of the regions.
In this scenario, limit-cycle solutions emerge.This kind of solution, after a transient time due to the dissipation, is trapped in stable and regular trajectories around φ ∆ = π or 0. In our case, a limit-cycle solution is classified as an n-cycle, where the number n represents the number of the points of the trajectory having zero velocity, i.e. the number of times the trajectory intersects the zero velocity axis, on one side of the phase space.The n-cycles collected in the stability diagram are labelled with different colours (see figure 4), while in figure 3 the phase space portraits of six different n-cycles are shown.For instance, we analyze the trajectory in figure 2d, having one intersection with zero velocity axis on the left and one on the right side of the point φ ∆ = π: it represents a 1-cycle.The system dynamics exhibits complex cyclic trajectories, such as 5-cycles (see figures 3b and 3e) and 9-cycles (see figures 3c and 3f).Furthermore, we find that n-cycles exhibit different topological properties and they can be distinguished through the turning number t and the winding number w.The turning number represents the number of loops in a cyclic trajectory, which in our case is equal to the number of intersections n.The winding number counts the number of oscillations around a specific phase space point, which in our case can be (0, 0) or (π, 0).We note that n-cycle can occur both with t = w and t > w (see figure 3) for different values of the parameters.
We estimate numerically the number of intersections n, working on the phase space trajectories of the n-cycles.First, we estimate their period through the analysis of the component in the frequency space, i.e. through the peaks of the Discrete Fourier Transform of φ ∆ (t) − φ∆ , where φ∆ (t) is the mean value of the phase: the smallest positive frequency ω * is inverted for the determination of the period T * = 2π ω * : n is then determined by counting the times in which the velocity φ∆ (t) crosses the φ = 0 axis in one period T * .The results of the classification are the n-cycle regions in figures 4a-4c, where each n-cycle is labelled with a different colour in the stability diagram.In particular, it is possible to identify three main subregions with 1-cycle, 3-cycle, and 5-cycle.However, narrow stripes of limit-cycle solutions with higher numbers of nods occur for some specific values of the parameter ε, as reported in figures 4b and 4c.For instance, in a neighbourhood of ε = 0.385 a large band of 9-cycles appears.
For the case φ∆ (0) ̸ = 0 we have obtained two basins of attraction (see figure 5), where a basin of attraction consists in the collection of all the initial conditions leading to a specific asymptotic attractor, obtained for constant values of the parameters.The estimation of the Hausdorff dimensions of the basins of attraction follows the same procedure of the case φ∆ (0) = 0 previously discussed.

Figure 4 :
Figure 4: (a) Stability diagram for a parameters space comprised of 896 × 2048 points of φ ∆ ∈ [ π 2 , π] in the horizontal axis, ε ∈ [0.3, 0.56] in the vertical axis and φ∆ (0) = 0.The diagram is divided into four regions by the main attractors, each labelled with a different colour: black stands for unstable, dark blue for 0-stable, blue for π-stable and the remaining colours represent the different n-cycles, which combine to form the limit-cycles region.The most likely are 1-cycle (light blue), 3-cycle (green) and 5-cycle (light green), divided into three different subregions.However, the limit-cycle regions exhibit stripes and domains exhibiting higher-order n-cycles.(b) Inset of (a), where a region of 2-cycles (cyan) appears along line A, standing for ϵ = 0.509253.(c) Inset of (a), exhibiting subregions of 6-cycles (yellow), 9-cycles (orange) and a narrow stripe of 15-cycles (red) for ε = 0.373 .Line B is for ε = 0.361176755.

B. 1 ) 1 ℏ 8 ℏ
from which we derive Hamilton's equations for the variables at nodes i = g, Iṅi = {n i , H} = − 1 J1 sin(φ I − φ g ) + E J2 sin(φ I − φ o )] ,(B.3a) ṅg = E J1 sin(φ I − φ g ), (B.3b) φI = E C2 (n I + n g − N ), (B.3c) φg = 8 with a voltage drop V g : the last term in the expression of N represents then the amount of Cooper pairs stored on C ∥ .We can write the number of Cooper pairs on g and I in terms of the numbers of Cooper pairs n 1 and n 2 on the two JJ capacitances as n g = n 1 + C ∥ Vg 2e and n I = n 2 − n 1 .Substituting these expressions in equation B.4, we obtain C.2) j ∈ {g, I}.Using the second Josephson relation, equation (C.1b), we can reformulate the equation for the current through the junction as

Figure D1 :
Figure D1: Estimation of the fractal dimensions δ H via linear best fits (blue lines) of the four different regions in figure 4a.The different linear fits correspond to (a) 0−stable solutions (d H ≃ 1.327), (b) unstable solutions (d H ≃ 1.520), (c) π−stable solutions (d H ≃ 1.307),(d) limit-cycle solutions, comprised of all the different n-cycles (d H ≃ 1.201).

Table 1 :
List of the n-cycles encountered in the stability diagram of figure4with their absolute, relative frequency and fractal dimension.