Anomalous Phase Dynamics of Driven Graphene Josephson Junctions

Josephson junctions with weak-links of exotic materials allow the elucidation of the Josephson effect in previously unexplored regimes. Further, such devices offer a direct probe of novel material properties, for example in the search for Majorana fermions. In this work, we report on DC and AC Josephson effect of high-mobility, hexagonal boron nitride (h-BN) encapsulated graphene Josephson junctions. On the application of RF radiation, we measure phase-locked Shapiro steps. An unexpected bistability between $\pm 1$ steps is observed with switching times on the order of seconds. A critical scaling of a bistable state is measured directly from the switching time, allowing for direct comparison to numerical simulations. We show such intermittent chaotic behavior is a consequence of the nonlinear dynamics of the junction and has a sensitive dependence on the current-phase relation. This work draws connections between nonlinear phenomena in dynamical systems and their implications for ongoing condensed matter experiments exploring topology and exotic physics.

The ground state wavefunction of a superconductor is endowed with an emergent phase on the macroscopic scale [1]. Josephson pointed out that two superconductors separated by a non-superconducting weak link produces a supercurrent that depends on the phase difference between the superconductors [2]. Today, eponymously known as the Josephson effect, it remains a hallmark of the superconducting state and has led to a large number of scientific and technological advances [3][4][5]. As one of the few phase-sensitive devices available in condensed matter physics, it has found use in the investigation of novel superconducting materials [6] and more recently in the identification of topological superconductors and Majorana bound states [7][8][9][10][11][12]. In particular, the Josephson effect in the presence of RF radiation (the AC Josephson effect) has been instrumental in advancing our knowledge of underlying physics in these systems. The incorporation of high-quality, tunable two-dimensional materials into a Josephson junction (JJ) now allow for these nonlinear devices to be probed in previous unexplored regimes. As the JJ becomes more prevalent in use both as a detector of novel materials properties and exotic circuit elements, a thorough understanding of the role of nonlinear dynamics in these devices is essential.
The understanding of JJs driven by externally applied currents and radio frequency (RF) radiation rests on simple circuit models that capture the superconducting phase dynamics. The relationship between the supercurrent (I) and the phase difference across the JJ (φ) is encoded in a current-phase relationship (CPR), which subsumes the microscopic physics of the weak-link [13]. The dynamics of a JJ then follow the evolution of a phase particle in an effective potential dependent on the physical parameters of the junction, such as normal state resistance, capacitance and externally applied currents [14,15]. At the heart of such models is a dynamical system described by a set ordinary differential equations made nonlinear by the CPR [1,16].
Nonlinear dynamical systems are interesting in their own right because of the rich set of phenomena they display. The phase dynamics of a driven JJ is similar to that of a driven, damped pendulum -a system known to be chaotic. Nonlinear effects in JJs have been studied extensively in theory [17][18][19][20][21][22][23]. As a result, a zoo of nonlinear effects including chaos, intermittency and strange attractors can occur [24], and these effects should be manifest in JJs. Hence, it is critical to detail these phenomena to have a thorough understanding of the physics of JJs and to ensure these effects do not masquerade as any exotic effects that could erroneously be ascribed to the novelty of the weak link material.
In this paper, we draw a bridge between experiments on graphene JJs and the unexpected consequences of the nonlinear phase dynamics. Graphene encapsulated in hexagonal boron nitride (h-BN) offers a platform to fabricate junctions with ballistic carrier transport. Moreover, gate electrodes and magnetic fields can be used to tune the junction parameters. We focus on measurements of the junction resistance as a function of a DC current bias and RF power, referred to as Shapiro diagrams. Measurements the Shapiro diagram show strong deviations from the expected behavior. In particular, portions of these diagram disappear in certain regimes, a result of chaotic behavior possessing timescales many orders of magnitude slower than have ever been observed before in JJs. This slow time scale allows for this nonlinear effect to be tracked in real time. Through comparison to simulations, the time scale is shown to depend sensitively on the CPR of the junction, becoming large for skewed CPRs like that predicted for weak links materials possessing a Dirac (linear) energy dispersion. These results should inform future experiments using RF radiation to probe exotic materials in JJs [7][8][9] and devices that exploit JJs in RF environments [25,26].

DEVICE CHARACTERISTICS
h-BN encapsulated graphene is fabricated into JJs, proximitized with molybdenum-rhenium alloy (MoRe) superconducting leads [ Fig. 1(a)]. MoRe is a type-II superconductor with an upper critical field of 8 T. Previously, such junctions have been used to study supercurrent in the quantum Hall regime and have been shown to have ballistic carrier transport [27]. The junction length along the direction of current flow is ∼500 nm and the width is ∼2.7 µm. In addition to MoRe contacts, we also have Cr/Au top gates and a backgate to electrically tune electron density. The top gates are kept grounded for the entirety of this experiment.
The device was cooled down in a cryostat with a base temperature of 50 mK. The magnetic field B was applied out of plane and perpendicular to the junction. Four terminals at the chip level measure the JJ in a twoterminal geometry (which removes the series resistance of the wires down the cryostat). The differential resistance R = dV dI is measured with lock-in amplifier at a frequency of 19 Hz and excitation of 1 nA as a function of back gate voltage (V BG ) and the applied DC current bias (I dc ). Fig. 1(b) shows the V BG dependence of the critical current. The charge neutrality point (CNP) is found to be at V BG = −0.75 V, corresponding to the lowest critical current and the largest normal state resistance. Fig. 1(c) shows the Fraunhofer pattern R(B) at the CNP. The regularity of the pattern implies the supercurrent is distributed uniformly along the contacts at the CNP.

MEASUREMENT RESULTS
In order to probe the phase dynamics of the JJ, R is measured as a function of I dc and RF power. When junctions are periodically driven, phase locked solutions produce a a fixed voltage across the JJ. These are known as Shapiro steps [28], and show up as zero differential resistance in a lock-in measurement. The measurement of the Shapiro diagram yields information about the CPR and has been used previously to study weak links of topological materials [7,8,[10][11][12][29][30][31].  Fig. 2(a) shows a Shapiro diagram at RF frequency (f ) of 6.350 GHz and V BG = −0.75 V. We observe regions of zero differential resistance corresponding to Shapiro steps separated by resistive transitions between steps. This Shapiro diagram is distinct in two principal ways from previous measurements on graphene [32] and other materials [7,8,[10][11][12][29][30][31]. The first is the resistive transitions are elongated node-like regions, not single points that are expected for the Bessel function step width as a function of RF power [4]. These elongated nodes are a generic feature of underdamped JJs (explained below) and can be used to infer underdamped behavior even when samples lack hysteresis in measurements of a I − V curve. Second, the first node separating the ±1 steps at I DC =0 disappears (indicated by the white box) for certain RF powers. Switching behavior on a timescale slower than the integration time for a lock-in measurement can lead to such a feature. Since the integration time used was around a second, this implies a switching timescale on the order of seconds. Fig. 2(b) shows a closeup of the broken node region. In addition to the broken node, the details of resistive transition is dependent on the junction parameters. Fig. 2(c) shows the changes when I DC is reduced by B, where additional structure in R is observed (white arrows). For a more extensive visualization of the changes to the broken node region in response to small changes in f , V BG and B, see the Appendix A.

Shapiro Diagram and Broken Node
To investigate the broken node region further, we switch to direct DC voltage measurements as a function of time. Figs. 2(d,e) show a histogram of the voltage values measured as function of I dc and RF power. For a given I dc or RF power, we measure the DC voltage V (t) at regular time intervals. These values are then binned. The counts are reported as the color in Figs. 2(d,e). When the count distribution is bi-modal, a bistability between the two values can be inferred. We clearly observe a bistability between the two ±1 Shapiro steps at the broken node as a function of I dc  However, in the transition region between −1 and +1, we see a bistable region. V (t) randomly switches between ±1 steps. τ is defined as the average time V (t) spends on a particular step. τ is extraordinarily long -on the order of seconds -much longer than any junction timescales. For comparison, a typical timescale associated with JJ dynamics is the characteristic frequency of the junction, defined as f c = 2eIcR , where I c is the critical current and R is the normal state resistance. f c for this junction is ∼ 50 GHz. Fig. 3(b) shows the dependence of τ on I dc . τ changes by almost three orders of magnitude in the bistable region. A similar scaling of τ now as a function of RF power is shown in Fig. 3(c). Note the log scale on the τ axis in Figs. 3 Such scaling behavior is unexpected. One would expect that the dynamics smoothly change between ±1 steps. The scaling with RF power also shows a maximum in τ at around −7.4 dBm in Fig. 3(c) and quickly decays away from this maximum. These observations point out that it is crucial to understand exact nature of the phase locked solutions though a model, which we do in the following section.

RCSJ MODEL
A model for the evolution of φ(t) under external DC and RF drives can be used to understand the Shapiro diagram. A simple circuit model for JJs, called the resistively and capacitively shunted junction (RCSJ) model [ Fig. 4(a)] [14,15], has been widely used for understanding phase dynamics in JJs. In this section, we describe the parameters for our model and the modifications required when graphene serves as the weak-link. φ can be related to experiment by calculating voltage across the junction using the Josephson relation V = 2e dφ dt [2]. A change from being entirely on −1 to +1 Shapiro step is seen. Near zero current bias, the voltage stochastically switches between the two steps. A long switching timescale, of order 10 seconds, is observed. (b, c) Dependence of the switching timescale τ on I dc and RF power measured in the region of the broken node.

Model Definition
In junctions with a conducting weak link, the supercurrent flow is mediated by Andreev bound states (ABS), formed by Andreev reflection at the two weak linksuperconductor interfaces [33,34]. The position of ABSs depends on φ, the electron/hole dispersion in the weak link region and on the coupling strength between the superconductors and the weak link. For a short and wide ballistic graphene junction at the Dirac point, the CPR is [16], where ∆ is the superconducting order parameter, W is junction width perpendicular to current flow and L is the junction length. The numerical factor of 2 1.33 ensures that the maximum as a function of φ is I c . See Fig. 4(b) for a plot of this CPR. Once the CPR is established, the phase dynamics can be modeled by RCSJ model and the governing equation is: where C is the junction capacitance, R is the normal state resistance, and I dc and I ac are the externally applied DC and AC currents. The equation can be cast into dimensionless form by transforming t → t 2eIc C : where Q = 2eIcR 2 C is called the quality factor and Ω = ω C 2eIc is the normalized frequency. Depending on the value of the capacitance and consequently the Q value, the junction can be overdamped (C → 0 ∴ Q → 0) or underdamped (Q 1).
An equivalent first order autonomous form dy dt = F (y) can be obtained by defining y = (φ, dφ dt , Ωt). This implies that the phase space for this system has to be three-dimensional, and hence solutions are not limited to fixed points and limit cycles, but strange attractors and chaotic effects are possible (Poincare-Bendixson theorem [24,35]). On the other hand, in the limit C → 0, the second order term vanishes and we are left with two-dimensional phase space, where the dynamics is restricted to fixed points and limit cycles. In the C → 0 limit, we recover the overdamped or resistively shunted junction (RSJ) model, whose dynamics are wellunderstood. See the Appendix B for a discussion on Shapiro diagrams in the overdamped limit.
The flow on y defined by F (y) has a Jacobian equal to −1 Q implying a dissipative system. It can be shown that in the presence of dissipation, the dynamics of a JJ with a sinusoidal CPR is equivalent to the driven pendulum. This correspondence to nonlinear dynamics of driven-damped pendulum has been explored before in simulations and connections to JJ experiments have been drawn, as we elaborate on in the next subsection.

Previous Work on Underdamped, Driven Chaotic Systems
Early motivation to consider the nonlinear dynamics of Josephson junctions arose from the large values of electric noise observed in experiment. It was observed that JJbased parametric amplifiers resulted in broad-band voltage fluctuations [36] with a noise temperature on the order of 10 4 K, which was ascribe to the chaotic behavior of Eq. 4 [17]. It was suggested that this arises from the presence of a strange attractor in phase space for certain set of parameters characterizing the junctions.
Experimentally, Poincare maps and bifurcations present in the parameter space were first measured with phase-locked electrical loops whose dynamics are similar to JJs and the driven simple pendulum equation [21]. A similar system was used to investigate symmetry breaking in solutions of the equation as a period-doubling route to chaos. Chaotic and intermittency was reported in measurements of bifurcation cascades of phase-locked loops [37]. Intrinsic and noise-induced intermittency due to crisis has been been studied with numerical simulations. The authors reported an intermittent switching behavior between positive and negative voltage solutions with characteristic power-law type relation for the switching timescale [38]. The basin boundary for these solutions was found to be fractal.
Extensive numerical simulations were carried to investigate the chaotic behavior of JJs and the effect of noise [18,19]. In experiments with a Josephson junction, a correlation was found between the noise and the fractal dimension of the boundary between periodic solutions [39,40]. Despite the extensive numerically simulations of the driven pendulum system, experimental verification remain scarce and are restricted to qualitative comparison [40][41][42][43].

Simulation Results for the Underdamped RCSJ Model
We now describe the simulation of a Shapiro diagram under the RCSJ model. The voltage across the junction is given by V (t) = 2e dφ dt . Since the voltage is measured in the DC limit, the average voltage is given as: where T = 2π Ω is the drive period and n is large number to average out transient behavior. For a phase locked solution, φ advances by an integer number of 2π period, and dφ dt at the m th Shapiro step is given by mΩ. The DC limit is achieved when n → ∞. It is also useful to define a window average for a a fixed number of periods which we denote as V n = dφ dt n . We will use both the DC limit and fixed window averages in the following sections. The values of V n = dφ dt n are normalized in units of Ω so that m th Shapiro step corresponds to V = dφ dt n = m. It is also important to consider the intial conditions for solving the ODE given by Eq. 4. Since, it is of second order, a unique solution requires two initial numbers and we use the initial condition (φ(t = 0),φ(t = 0)) = (0, 0) unless otherwise stated. In order to understand the effects of noise, a Wiener process in time is added to the RHS of Eq. 4. This models a Gaussian white current noise source with zero mean. As a result, the RCSJ model with noise becomes a stochastic differential equation (SDE) as follows: where η(t)dt is a Weiner process and σ n controls the noise strength. It is useful to think of σ 2 n as being proportional to a energy or temperature scale, which would hold exactly if the noise sources were due to Johnson-Nyquist noise in resistors. In simulations with noise, we integrate this SDE with the Euler-Maruyama method.
To understand our experimental results, we consider junctions in the weakly underdamped limit where Q ≈ 1. Since the junction capacitance is not measured directly, we indrectly infer that the value of Q must be close to 1. Fig. 4(c) shows a Shapiro diagram measured at 6.350 GHz, V BG = +0.25 V and magnetic field of 0.7 mT. We see that results are different from those at CNP and the broken node is not present. Fig. 4(d) shows a simulated Shapiro diagram with Q = 1 and Ω = 0.6; qualitative reproduction of the measured features is observed. In particular, the contours separating the steps, meet in points (±1 steps) or show a small gap between them. Since the values for V BG and B chosen in Fig. 4(c) cause a reduction in R and I c and hence a reduction in Q, we expect the dynamics at CNP and zero magnetic field to be described by a model with Q > 1. Fig. 5(a) shows a simulated Shapiro diagram for Q = 1.5 and Ω = 0.5. As opposed to diagrams in the overdamped limit (Appendix B), there are many distinctive features. The resistive contours for the steps meet in elongated notes as opposed to single points, similar to the experimentally observed features [ Fig. 2(a)]. For a more detailed comparison between the fine features of the measurements and simulations, see the Appendix C.

Discrete symmetry of underdamped RCSJ model
We now focus on the explanation on the broken node feature. Since the node occurs on the I dc = 0 line, we consider the following equation: Eqn. 9 has a discrete symmetry among the possible solutions. For φ → −φ and t → t + π Ω , the equation is left invariant. Note that, we also need i(φ) = −i(−φ) and this holds from time-reversal symmetry [13]. This symmetry implies, Under this operation, for given values of Q, Ω and i ac , if V = m is possible, so is V = −m. Hence if m = 0, it is possible for the system to have two distinct solutions with an inequivalent observable (the DC voltage) in a measurement. Which of these two phase locked solution is achieved depends on the initial conditions (φ(t = 0),φ(t = 0)). It is convenient to think of a basin structure [24] for each of the two solutions, where each initial condition eventually settles into one of these two possible solutions. We wlll later see how an interplay between the basin structure and the noise level can lead to the switching characteristics measured in Fig. 3(b,c).

Poincare Maps and Bistability
A handy construction for elucidating the behavior of a dynamical system is a Poincare map or section. For a dynamical system, it is the intersection of a periodic orbit with a lower-dimensional subspace, transverse to the flow of the system. In our case, for the RCSJ model, we consider (φ(t), dφ dt (t)) for t = n 2π Ω where n is an integer. In Fig 5(b), the Poincare map for dφ dt with parameters Q = 1.5 and Ω = 0.5 is plotted. i dc fixed to 0 and i ac is varied in a range around where the ±1 steps meet. The Poincare map has a rich structure with periodic attractors, separated by windows where the dynamics causes dφ dt to spread uniformly over a range of values. The behavior in these windows is similar to the periodic doubling bifurcation cascade in the logisitc map [24] and we refer to these windows as chaotic windows. In general, the phase space structure of the model along i ac axis consists of periodic attractors separated by chaotic windows.
The phase dynamics φ(t) is challenging to measure directly in an experiment and the only quantity accessible is the DC voltage. To this end, we now discuss Poincare map with an average over a fixed window. We simulate dynamics for 10000 periods, average over 200 consecutive periods and plot the values taken by the window averaged for each value of i ac . We call this the DC limit of the Poincare map and it is pragmatic in explaining our observations. Fig. 5(c) shows the DC limit of the map from Fig. 5(b). The chaotic windows average out to noisy regions while the periodic attractors attain their DC values. An attractor to the +1 step is found to occur for i ac ∈ (1.23, 1.33). Fig. 5(c) shows a closeup of the DC Poincare map of this attractor. As i ac is increased in this range, we move from being on an attractor with dφ dt = 0 to dφ dt = +1 and back to 0 again, separated by noisy regions due to the chaotic windows.
The symmetry argument made in the last section, implies that the attractor in the range i ac ∈ (1.23, 1.33) also has solutions corresponding to the −1 step. Which attractor the system settles in depends on the initial conditions. However, the presence of noise removes any bias corresponding to the initial conditions and causes the system to switch between ±1 attractors. Fig. 5(e) show this effect, where a noise term with σ n = 4×10 −2 is included. It is useful to recall the basin structure, where each set of initial conditions lead to one of the attractors. The presence of noise causes a drift in initial conditions after every period and so it becomes possible for the solution to switch from one attractor to another. Fig. 5(f,g) show the histogram of experimentally measured V (t) values measured as a function of RF power in the broken node region (f, B) = 6.322 GHz, B=0 (f) and 6.350 GHz, B=0.5 mT (g). We see that the observed behavior of the DV voltage in the broken node region is qualitatively similar to the dynamics in the DC limit of Poincare maps in Fig. 5(b-e). An attractor to ±1 steps is sandwiched between attractors to the 0 step separated by noisy regions corresponding to the chaotic windows leads to the broken node feature observed in the experiment.

Scaling of switching time
Having understood the origin of bistability in simulation, we now focus on the changes in switching time τ as a function of DC and AC currents. We reported that τ undergoes changes over orders of magnitude when I dc and RF power are changed in Fig. 3(b,c). We now reproduce similar scaling in simulation. Fig. 6 shows the changes in τ as a function of i ac , noise level σ n and i dc . We use the same set of parameters and the ±1 attractor discussed in the last section.
In Fig. 6(a,b), we plot the switching time τ for the +1 and −1 steps respectively. We also compare the effect of the graphene CPR to a simple sin(φ) CPR. While fixing all other parameters, the CPR of graphene causes the switching time to be much larger. This feature also occurs with the scaling of switching time as a function of noise level and i dc [ Fig. 6(b-e)]. This can be attributed to the graphene CPR leading to an effective potential for φ which has steeper and deeper minima, making escape for a given noise level less likely. This increase in τ resulting from the CPR of graphene serves as an explanation why the bistability might have been harder to directly measure in other underdamped systems.
The presence of chaotic windows for a range of parameters complicates an analytic approach that might elucidate the origin of the observe scaling. Lacking any exact theoretical basis, we now give a heuristic argument for the scaling of the switching time. It is helpful to recall the basin structure picture, where the noise causes a drift of initial conditions, allowing for the solution to switch from the basin of one attractor to another. Formalizing this notion, suppose, the noise causes a diffusion at time t as ∆φ(t) ∼ (σ n √ t) ν , where the term σ n √ t has been chosen because it is dimensionless and ν is a diffusion exponent. ν = 1 corresponds to standard Brownian motion, where net displacement is proportional to √ t. Suppose f (i ac , i dc ) represents a scale in the basin structure; when initial conditions are changed by this scale, a switch from one attractor to another occurs. We can therefore write a condition for the switching timescale τ : Using the diffusion expression for ∆φ, the following expression is obtained: (e,f) Effect of i dc on τ for the +1 (e) and −1 (f) steps. As i dc is increased from negative to positive values, the solution prefers to stay more on +1 step than on −1 step; the corresponding changes in τ are evident. In all the plots, the effect of the graphene to sinusodinal CPR is compared. Though both CPRs show a scaling of τ with iac, the overall magnitude is larger for the graphene CPR. A window size of 10 periods was used to take the DC limit and the system was simulated for a total of 10 5 periods.
The details of the function r depends on the exact nature of the CPR, noise level σ n as well as the parameters Q and Ω. Though we don't have a way to calculate f (i ac , i dc , Eq. 12 show why such relationships might exist in the first place. Further theoretical work is necessary to elucidate the scaling relationship of τ in the model.

CONCLUSION AND OUTLOOK
In summary, we have reported on the AC Josepshon in h-BN encapsulated graphene JJs. A number of features distinct from previous measurements from unencapsulated graphene junctions as well as other systems have been observed. The resistive contours separating two steps meet in elongated nodes as opposed to points expected from RSJ model. Incorporating the capacitance and considering the full RCSJ model in the weakly underdamped limit leads to a provide an origin of these elongated nodes. Moreover, we observe a bistability at the first elongated node. The switching timescale is extraordinarily long, on the order of seconds, and is much slower than any timescale relating to the junction dynamics. A change in τ over almost three orders of magnitude is observed with changes in external applied DC currents and RF power.
Simulation of the RCSJ model reveals a rich dynamical structure where periodic attractors are interleaved with chaotic windows. In the DC limit, this leads to attractors separated by noisy regions. In particular, we show the presence of an attractor for ±1 steps whose dynamics in the presence of noise is similar to the bistability observed in the broken node region. Simulation the switching time in the bistable region as a function of i ac , i dc and noise level offers an explanation of the scaling behavior observed in the experiment. We also show how the graphene CPR effects the overall magnitude of the switching timescale. We offer a heuristic argument for the existence of such scaling from the interplay between the basin structure and the noise level.
A driving theme in contemporary condensed matter research is the realization of Majorana particles and topologically non-trivial ground states [44][45][46][47]. One of the principal ways to detect topology in these systems is the 4π Josephson effect and has been used to this effect in certain experiments [7][8][9]. Yet, chaotic behavior can mimic this effect: simulation above show that period doubling is possible for certain values of i ac . Our work points to a key feature which can indicate the potential for chaotic behavior (i.e. underdamped) in a JJ -the elongation of nodes in a Shapiro diagram. This seems to be more accurate that the standard means of inferring underdamped junctions, hysteresis, which may be absent in slightly underdamped (Q ∼ 1) JJs. For example, work on HgTe JJs estimate an overdamped junction from device parameters [48] and as a result rule out the possibility chaotic behavior in the JJs [8]. Yet observe elongated nodes seem to contradict this assessment. Further, gaps in resistive nodes have been ascribed to a 4π-periodic contribution to the CPR [10,12]. In this work, we observe this gap to arises entirely from nonlinear behavior in the junction [ Figs. 4(c,d)]. As the coupling between the superconductor and the material forming the weak link is improve and as the mobility of the weak link is increase, larger values of I C and Q are expected. Hence, a systematic delineation of phase dynamics of driven JJs is crucial for using the AC Josephson effect as a probe of topology and exotic physics.
The underdamped RCSJ model is dissipative with a strength parametrized by 1 Q . In the limit Q 1, the RCSJ model becomes a conservative Hamiltonian system. We have let φ be a classical variable, though it is possible to quantize the system. It would be interesting to study the quantum dynamics of a Hamiltonian system resulting from the RCSJ model. Further theoretical work in this direction is necessary to understand the full parameter space of RCSJ model in this limit.
In conclusion, driven JJs where nonlinear effects can lead to the presence of periodic attractors and chaos offers a new route to study nonlinear phenomena. The graphene junction presented in this work is a tunable platform well-suited for such studies. We hope it will be a useful toolbox to study novel phenomena at the intersection of nonlinear dynamics and condensed matter physics.

APPENDIX A: VARIATION IN SHAPIRO DIAGRAM OF THE BROKEN NODE
In order to understand the origin of the broken node region, we make small changes in the external parameters that affect the junction physics. Fig. 7 shows the changes in the broken node region for small changes in f , V BG and magnetic field. In Fig. 7(a), from left to right, f increases as from 6.310 GHz, 6.320 GHz, 6.330 GHz and 6.370 GHz while V BG is fixed to CNP and B=0. In Fig. 7(b), from left to right, V BG changes as −0.55 V, −0.85 V, −1.05 V and −1.15 V while f is fixed to 6.350 GHz. In Fig. 7(c), we ramp the magnetic field to 0.3 mT, fix f to 6.350GHz and change V BG as −0.45 V, −0.55 V, −0.65 V and −0.75 V.

APPENDIX B: SIMULATED SHAPIRO DIAGRAM IN THE OVERDAMPED REGIME
A junction is said to be overdamped if the capacitance in the RCSJ model is small. In the limit C = 0, we can rewrite the evolution equation for φ as: 2eR dφ dt + I(φ) = I dc + I ac sin(ωt) (13) We recast this equation into a dimensionless form by transforming t → t 2eRIc : dφ dt + i(φ) = i dc + i ac sin(Ωt) (14) where Ω = ω 2eRIc is a normalized frequency. Note that as compared to the full RCSJ model (eqn. 4), this system defines a flow in a two-dimensional phase space and its dynamical behaviour is constrained to take on fixed points and limit cycles by the Poincare-Bendixson theorem [24,35]. In fact, it is possible to solve this equation analytically, and the Shapiro step widths take a form given by Bessel functions [1]. For completeness, we reproduce Shapiro diagrams as a function of varying frequency Ω in Fig. 8. Such overdamped diagrams have been observed before in graphene junctions [32].

APPENDIX C: UNDERDAMPED SHAPIRO DIAGRAMS
In addition to the broken node feature described in the main text, the RCSJ model in the weakly underdamped limit (Q ∼ 1) can also be used understand the features in the resistive contours separating the Shapiro steps. Fig. 9a and c show a simulated Shapiro diagram with DC voltage and differential resistance respectively. Fig. 9b shows linecuts at fixed i ac values. In addition to integer steps, fractional steps are also seen in the transition region between two integer steps. The parameters used for this simulation were Q = 1.5, Ω = 0.5, and this is a closeup version of the diagram in Fig. 4a near the intersection of the ±1 step. Fig. 9d and e show measured Shapiro diagrams at f of 4.000 GHz and V BG set to CNP and −0.45 V respectively. The opening and closing up of the resistive contours in the measurement is qualitatively similar to the fine features observed for i a c ∈ (1.2, 1.6) in the simulation.

ADDITIONAL DATA ON SWITCHING TIME SCALING WITH RF POWER
The nature of the scaling relationship of switching time τ with RF power on the bistable step varies with applied frequency and the magnetic field. Fig. 10(a) and (b) shows τ vs RF power for (f, B, I dc ) = (6.338GHz, 0mT, 0nA) and (f, B, I dc ) = (6.35GHz, 0.5mT, 1nA) respectively. The back gate was set to CNP.