# Noise Effects on birhythmic Josephson Junction coupled to a Resonator

###### Abstract

We study the effect of noise on a Josephson junction that, coupled to a linear resonator, can oscillate at two frequencies. To establish the global stability of the attractors, we estimate the position of the separatrix, an essential information to establish the stability of the attractor for this multidimensional system, from the analysis of the mean first passage time. We find that the frequency locked to the resonator is most stable at low bias, and less stable at high bias, where the resonator exhibits the largest oscillations. The change in the birhythmic region is dramatic, for the effective barrier changes of an order of magnitude and the corresponding lifetime of about seven decades.

###### pacs:

05.40.-a;05.45-a; 05.45.Xt;05.40.Ca; 85.25.Cp## I Introduction

The contemporary presence of two frequencies for the same set of parameters, or birhythmicity, is encountered in some biochemical systems decroly82 (); morita89 (); haberichter01 (); sosnovtseva02 (); abou-jaoude11 (), nonlinear electronic circuits enjieu07 (); zakharova10 (); yamapi10 (); ghosh11 (); yamapi12 (); yue12 (), and extended distributed systems stich02 (); casagrande05 (). The experimental observation of birhythmic systems is, however, less frequent hounsgaard88 (); geva-zatorsky06 (); ventura07 (); gonzalez08 (). In this context the superconducting circuit consisting of Josephson Junctions (JJ) coupled to a cavity hadley88 (); filatrella92 (); ozyuzer07 (), as in Fig. 1, represents a preeminent example of birhythmic system that is also interesting for applications. The coupling among the junctions is supposed to be provided by a resonant cavity grib06 (); gross13 (); grib13 (), thus when all the junctions are entrained it is essential to have a large current in the cavity, such that the junctions can be entrained through the current in the resonator grib13 (). The state with a large current coexist with a state at lower power; the two states are clearly characterized by two different frequencies. This is the essential feature of birhythmicity, the coexistence of two attractors characterized by two different amplitudes and frequencies: depending on the initial conditions, the system can produce oscillations at two distinct periods. Being the attractors locally stable, the system would however stay at a single frequency, the one selected by the choice of the initial. Thus the system exhibits an hysteretic behavior: the displayed frequency depends upon the initial conditions. In the presence of noise the system can switch from an attractor to the other under the influence of the random term. Birhythmicity is therefore a nonlocal phenomenon that cannot be investigated by linear analysis landauer78 (). In this work we aim to determine the global stability of the two states at different frequencies and of the on Fig. 2, to ascertain the birhythmic properties induced by the circuit. From the simulated of Fig. 2 it is evident that at the same bias point, e.g. , two frequencies appear , viz. and , depending on the initial conditions. The first frequency is reached increasing the bias current from zero on the Josephson supercurrent, while the second is obtained decreasing the current from high values on the resistive McCumber branch; the selection of the frequency actually displayed is thus determined by the initial conditions. The features of the IV depend upon other factors such as the number of JJs and the features of the resonator grib06 (); grib13 (). Also, heating effects are believed to be relevant for synchronization wang10 (), as well as coupling through charge transfer through the Josephson channel ovchinnikov13 (). In this work we consider the simplest case of a single JJ coupled to a high cavity, and we neglect heating, that occurs at a much slower time scale.

A switch from an attractor to the other is of central interest for devices based on synchronization of JJ through an circuit, in particular for BSCCO stacks for THz generation tachiki11 (). Applicationwise, it is undesirable an uncontrolled switch from the state locked to the (the high power generation) to the other (the low power emission) grib13 (). Unfortunately, the analysis of large fluctuations, as large as to carry the system from an attractor to another, is not easy, for it goes beyond the linear stability kautz94 (); lin11 () given by Lyapunov exponents lin11 (); tsang91 (). In equilibrium non dissipative systems, global stability is given by the time to escape from the energy potential at a given a noise level . Arrhenius law predicts that the average escape time exponentially depends upon the ratio between the energy barrier and the noise intensity hanggi90 ():

(1) |

In nonequilibrium systems, or when the potential energy is not available, a possibility is to reverse the logic and to define a pseudopotential energy barrier proportional to the logarithm of the lifetime dykman79 (); graham85 (); graham86 (); dykman94 (); dykman01 () , viz.

(2) |

In fact, under general assumptions, it can be postulated that the escape time between the two attractors exponentially depends upon a quantity (the pseudopotential energy) and it is inversely proportional to the noise intensity dykman79 (); graham85 (); graham86 (); dykman01 (). This approach has been used to determine the energy barrier of vortex motion kautz94 (); girotti07 () in Josephson systems, and has been employed to determine the energy barrier for anharmonic oscillators with cubic dykman94 () and quintic yamapi10 (); yamapi12 () nonlinearities. The same methods has been also used to investigate Shapiro steps kautz88 (). At variance with irradiated JJ where one frequency is given by an external drive, in the present system the system self-generates the two frequencies. Moreover, chaos can occur in rf-fields kautz88 (), as well as in several JJ coupled together kolachi13 (). Instead we prefer to focus on a simpler system, where the switch only occurs because of noise, between two otherwise (locally) stable attractors.

If the energy barrier is to be determined by means of the lifetime, as per Eq.(2), it is crucial to locate the separatrix between the two basins of attraction of the stable states. To determine the basins of attraction requires the knowledge of the initial conditions that lead to one or the other of the stable solutions, and therefore demands a detailed exploration of the phase space. However, being this exploration very difficult in the four-dimensional system of Fig. 1, we propose to exploit the fact that is a Lyapunov function graham85 () to estimate the separatrix. As will be shown in Sect. IV, the method we propose is capable to determine an effective threshold for the escape time, and therefore our approach constitutes a method for the estimate of the pseudopotential when the boundary of the basin of attraction is not exactly known. With this approach, we find that the stability of the attractor is not uniform: at the bottom of the step the trapping energy is high, and decreases at the top. This behavior is somehow counterintuitive, in that the global stability is enhanced when the frequencies of the two attractors get closer.

The work is organized as follows. In the next Section we describe an underdamped JJ coupled to a resonator and subject to external bias and noise. In Sect. III we discuss the locally stable attractors characterized by two frequencies, and how a transition from an attractor to the other can occur under the influence of noise. In Sect. IV we describe the method to locate the separatrix between the two attractors, an essential information to reconstruct the activation barrier. The methodological premises permit to determine the stability properties of the JJ in the birhythmic region. Section V concludes.

## Ii Model of a Josephson junction coupled to a resonator

Figure 1 schematically describes the model used in our analysis: an underdamped JJ connected in parallel to an resonator. Both elements are supposed in the temperature controlled vessel, while the bias current is supplied by a device at room temperature.

In this configuration the noise from the bias supply dominates respect to the Johnson noise from the resistors and . Alternatively, one could add a random term for each resistor, as done for instance in Ref. lin11 (). However, the noise is but a tool. Our goal is to determine the pseudoenergy; the principle of minimum energy graham86 (); kautz88 () assures that the contributions from the minimal trajectory determines the height of the trapping potential, and therefore one does not expect substantial changes with a different noise source.

The electrical model consists of the capacitor , the resistor , and the ideal Josephson element, connected in parallel. The nonlinear relation between the current and the gauge invariant phase difference across two superconductors:

(3) |

together with the Josephson voltage relationship

(4) |

determines that a JJ is an active oscillator that converts a dc current into an ac drive for the resonator. To derive the equations governing the system, we indicate with the current flowing through the circuit and with the charge on the capacitor. The JJ and the resonator are both biased by a current generator affected by a noise current that split in the current through the JJ element and the current through the . If we indicate with the current through the JJ resistor and with the current through the junction capacitance, we obtain the current balance:

(5) |

The Kirchhoff law for the loop voltage

(6) |

completes the model, that is thus described by two second order coupled differential equations:

(7) |

Introducing the Josephson frequency , Eqs.(7) can be cast in the normalized units and :

(8) |

where

The statistical features of the noisy term are determined by:

(9) |

The noise is due to an external source, therefore it does not obey the fluctuation-dissipation theorem and it is independent of the resistance. Equations (8,II) are simulated with the Euler algorithm fox88 (). Deterministic results have also been simulated with a Runge-Kutta algorithm. The curves have been obtained slowly increasing the bias current, with a step , and using the final state at the previous current step as the initial state for the increased (or decreased) current bias. At each current step a transient of about normalized time is discarded. The averages are also calculated over the same time. The time step is, through all simulations, for the Euler algorithm, and for the Runge-Kutta deterministic simulations. The stochastic results are also averaged over as many realizations as it is necessary to have reliable results. Finally, the Box-Mueller algorithm knuth69 () is used to generate Gaussian white noise from two random numbers and which are uniformly distributed on the unit interval . Thus for each step , is randomly distributed as follows:

(10) |

The system (8), depending on the initial conditions, can exhibit oscillations at two distinct periods. If the resonator and the JJ are weakly coupled (), the two frequencies and are substantially unperturbed and correspond to the resonant frequency of the and the unperturbed () frequency of the junction, respectively lin11 (); filatrella03 (). In contrast with previous studies filatrella92 (); filatrella03 (), this is the strong coupled limit (), and therefore shukrinov12a (); shukrinov12b (), as shown in Fig. 2 by the normalized voltage in correspondence of the applied current . The shift in voltage due to the interaction has been estimated in the limit case of a non dissipative resonator shukrinov12a (); in our normalizations it reads:

(11) |

The quantitative agreement is poor, as expected for a dissipative cavity. However, the pure cavity correctly predicts the trend towards an increase of the resonant frequency.

A word about normalization. First, in these units the normalized voltage and the normalized frequencies are expressed in the same units; therefore (see Fig. 2) the resonant frequency should be comparable with the characteristic frequency of the junction, thus:

Also, the resonator capacitance and the junction capacitance are connected by the relation:

For large coupling () the capacitance of the resonator should much larger than the junction capacitance. Finally, the relation between the resistance of the resonator and the resistance of the junction reads:

thus for high ( in these simulations) and large coupling (), we get . For underdamped junctions (here ) the resistance of the resonator is much less than the resistance of the junction. However, some care should be taken: in the equivalent circuit of Fig. 1 the resistance of the JJ is in parallel, while the resonator is modeled by series lumped elements.

The resonant step locked to the cavity is shown in more detail in Fig. 3: In the range the system stays on one or the other frequency, depending on the initial conditions (that are controlled by the bias sweep). To each frequency corresponds a different attractor, as will be analyzed in the next Section.

## Iii Attractors properties of birhythmic Josephson Junctions

Figure 4 displays the projection of the phase space in the plane . The branch locked to the resonator (characterized by the frequency ) – quite naturally – exhibits much larger excursions of the charge oscillations respect to the unlocked branch (characterized by the frequency ). In between, one postulates the existence of an unstable orbit with frequency that represents the separatrix. Fig. 4(ii) confirms that one can identify the attractor by the amplitude of the oscillations. In fact in Fig. 5 it is shown the behavior of the amplitude of the voltage (Fig. 5i ) and charge (Fig. 5ii) oscillations while the bias is sweeped. The amplitudes are defined as the largest excursions of the phase derivative (proportional to the JJ voltage) and of the charge (proportional to the capacitor voltage):

(12) |

In Fig. 5 it is evident that a sudden change of the amplitudes and occurs both for low and high bias. It is exactly this sharp change that we want to exploit to retrieve the escape rate. Let us consider the dynamics under the influence of noise, Fig. 6.

The attractors are deformed, but still well separated, see Fig. 6i; we can therefore tentatively locate the separatrix at . Being the system 4-dimensional the separatrix is a volume in 4-D space, whose projection in the plane ought not to be a line, and hence the dashed segment of Fig. 6 is but a rough approximation.

We postulate that when the charge passes the threshold a switch occurs to the other attractor, as shown in Fig. 7. In general the construction of the whole pseudopotential landscape to identify the separatrix requires the solution of a variational problem dykman94 (); kautz94 (). We instead adopt a simpler procedure based on the observation that the pseudopotential is a Lyapunov function graham85 (), and therefore becomes negative beyond the separatrix, as schematically shown in Fig. 8. In fact Fig. 7 illustrates that a sudden switch occurs when the fluctuations exceed a threshold, or when , i.e. when the system passes into the descending part of the pseudopotential in Fig. 8. However, exactly because of the switch towards the attractor, it is not necessary to accurately know the value of : any value behind leads to a similar estimate of the MFPT (Mean First Passage Time) hanggi90 (), see Fig. 9. We emphasize that the mean first passage time across any point in the vicinity of the separatrix has two distinct behaviors: i) it increases exponentially when the threshold point is set before the separatrix, and ii) it increases very weakly when the threshold is beyond the maximum of the potential. The different behavior is shown in Fig. 9, and therefore from the change in the slope of the MFPT we estimate the position of the separatrix. In summary, in the descending region beyond the separatrix (the dashed part of the pseudopotential) the system quickly runs ”downhill”, and the time elapsed in the dashed part is negligible respect to the time necessary to reach, under the influence of noise, the peak of the pseudopotential. This conjecture is confirmed by the MFPT with different choices of the threshold (see Fig. 9): there is a region where the average time is almost independent of the choice of the threshold. We conclude that the knee of the MFPT can be used as an effective separatrix to estimate the pseudoenergy activation barrier.

This is practically implemented in Fig. 10 for different values of the bias current. The linear relationship between the logarithm of the escape time and the inverse of the noise intensity offers the estimate of an effective energy barrier, see Eq.(2):

(13) |

Equation (13) is the main result of this part of the paper: to characterize with an activation energy the metastable states in the birhythmic region.

## Iv Energy barrier and lifetime of the induced step

In this Section we collect the results on the analysis of the birhythmic region of the IV curve in Fig. 2.

We use Eq.(13) to retrieve the behavior of the activation energy as a function of the bias , see Fig. 11. The energy barrier for low bias is large and the attractor is bounded in a stable well. When the current is increased along the step the energy barrier decreases, and almost disappears at the top of the step, where the frequency splitting is at a maximum, see Fig. 3. In the same process, the energy dissipated by the cavity increases, for the normalized power linearly rises along the step:

(14) |

(in the dimensionless equations we are using the voltage is expressed in units). At the bottom of the step () the power is low, while the energy barrier is at the maximum. In the region the effective energy barrier decreases of about an order of magnitude ( passes from to ). A relevant feature of Fig. 11 is that the change of the pseudopotential in the birhythmic region cannot be ascribed to a difference in the frequencies. In fact the pseudopotential is at a maximum when the difference is at a minimum. Thus, the switch from the resonant step back to the IV curve of the unperturbed dynamics occurs because the activation pseudopotential vanishes. The change is dramatic if one considers that time is normalized respect to , that is typically above 100GHz. A lifetime of the order of a second therefore entails a noise level as low as to reach , or . From the behavior shown in Fig. 10 one estimates for , and for . Put it another way, at a fixed noise level the lifetime decreases of at least seven decades when the bias current passes from to . It is also noticeworthy that the Arrhenius-like behavior implied by the existence of the pseudopotential greatly simplifies the numerical problem to find the noise intensity at which the desired lifetime is reached. In fact the relation Eq.(2) allows to extrapolate very long lifetimes from Fig. 10, whereas direct simulations of such long lifetimes are prohibitive.

## V Conclusions

We have found that the global stability analysis of JJ coupled to a resonator shows a striking change in the birhythmic region: the attractor characterized by a frequency locked to the resonator is most stable for low bias current, when the power dissipated in the cavity is small. The system is, unfortunately, less stable at the top of the step, when the current in the resonator is at a maximum and the two frequencies are most separated ( a similar conclusion has been reached for another nonlinear birhythmic system yamapi10 ()). Thus the analysis of large excursions, as large to drive the system from the desired attractor to another, indicates that the global stability is weak where most power is available. This observation is, in some sense, bad news for applications, inasmuch it shows that stability and high power are contradictory requirements. However, the detailed behavior of the global stability demonstrates that the deterioration occurs at the middle of the step, where the stability is still relatively high – see Fig. 11. From a more general point of view, we have shown that the stability of a dynamic state can be analyzed in terms of the pseudopotential also when the separatrix is not known and the variational approach is difficult to apply. Instead we propose a simpler method to (approximately) determine the position of the separatrix from the change in the slope of the MFPT.

A number of cautions are in order, however. In the first place, we have analyzed a single junction coupled to a cavity, while for applications such as BSCCO stacks one should consider many junctions tachiki11 (); welp13 (). Second, we are using lumped elements for both the junctions and the cavity, whereas a distributed description tachiki11 (); welp13 (); sakai93 (); madsen04 (); madsen08 () is more appropriated. Finally, thermal effects cause self-heating and back-bending of the usually associated to emission ozyuzer07 (); gross13 (), whereas we have here only addressed the effect of random fluctuations. Nevertheless, the calculations of this work point to a conceivable danger: that global stability properties are of crucial importance to determine the region of parameters where large power devices could possibly work.

## Acknowledgements

R.Yamapi undertook this work with the support of the ICTP (International Centre for Theoretical Physics) in the framework of Training and Research in Italian Laboratories (TRIL) for AFRICA programme, Trieste, Italy. He also acknowledges the hospitality of the Dipartimento di Fisica ”E.R. Caianiello” of the Università di Salerno, Fisciano, Italy.

The authors acknowledge partial financial support from PON Ricerca e Competitività 2007-2013 under grant agreement PON NAFASSY, PONa3_00007.

## References

- (1) 0. Decroly and A. Goldbeter, Proc. Natl. Acad. Sci. USA 79, 6917 (1982).
- (2) M. Morita, K. Iwamoto, and M. Sen, Phys. Rev. A 40, 6592 (1989).
- (3) T. Haberichter, M. Marhl, and R. Heinrich, Biophysical Chemistry 90, 17 (2001).
- (4) O.V. Sosnovtseva, D. Setsinsky, A. Fausboll, and E. Mosekilde, Phys. Rev. E 66, 041901 (2002).
- (5) W. Abou-Jaoudé, M. Chaves, J.-L. Gouz, PLoS ONE 6, e17075 (2011).
- (6) H. G. Enjieu Kadji, J. B. Chabi Orou, R. Yamapi, and P. Woafo, Chaos, Solitons and Fractals 32, 862 (2007).
- (7) A. Zakharova, T. Vadivasova, V. Anishchenko, A. Koseska, and J. Kurths, Phys. Rev. E 81, 011106 (2010).
- (8) R. Yamapi, G. Filatrella, and M. A. Aziz-Alaoui, Chaos 20, 013114 (2010).
- (9) P. Ghosh, S. Sen, S. S. Riaz, and D. S. Ray, Phys. Rev. E 83, 036205 (2011).
- (10) R. Yamapi, G. Filatrella, M. A. Aziz-Alaoui, and Hilda A. Cerdeira, Chaos 22, 043114 (2012).
- (11) X.L. Yue, W. Xu , L. Wang, and B. Zhou, Probabilistic Engineering Mechanics 30, 70 (2012).
- (12) M. Stich, M. Ipsen, and A. S. Mikhailov, Phys. Rev. Lett. 86, 4406 (2001); Physica D 171, 19 (2002).
- (13) V. Casagrande and A. S. Mikhailov, Physica D 205, 154 (2005).
- (14) J. Hounsgaard, H. Hultborn, B. Jespersen, and O. Kiehn, J. Physiol. 405, 345 (1988).
- (15) N. Geva-Zatorsky, N. Rosenfeld, S. Itzkovitz, R. Milo, A. Sigal, E. Dekel, T. Yarnitzky, Y. Liron, P. Polak, G. Lahav, and U. Alon, Mol. Sys. Biol. 2, 2006.0033 (2006).
- (16) A. Ventura, D. G. Kirsch, M. E. McLaughlin, D. A. Tuveson, J. Grimm, L. Lintault, J. Newman, E. E. Reczek, R. Weissleder, and T. Jacks, Nature 445, 661 (2007).
- (17) H. González, H. Arce, and M. R. Guevara, Phys. Rev. E 78, 036217 (2008).
- (18) P. Hadley, M. R. Beasley, and K. Wiesenfeld, Phys. Rev. B 38 , 8712 (1988).
- (19) G. Filatrella, G. Rotoli, N. Grønbech-Jensen, R.D. Parmentier, and N.F. Pedersen, J. Appl. Phys. 72, 3179 (1992).
- (20) L. Ozyuzer, A. E.Koshelev, C.Kurter, N. Gopalsami, Q. Li, M. Tachiki, K. Kadowaki, T. Yamamoto, H. Minami, H. Yamaguchi, T. Tachiki, K. E. Gray, W.-K. Kwok, and U. Welp, Science, 318, 1291 (2007).
- (21) B. Gross, J. Yuan, D. Y. An, M. Y. Li, N. Kinev, X. J. Zhou, M. Ji, Y. Huang, T. Hatano, R. G. Mints, V. P. Koshelets, P. H. Wu, H. B. Wang, D. Koelle, and R. Kleiner, Phys. Rev. B 88, 014524 (2013).
- (22) A. Grib, M. Mans, J. Scherbel, M. Buenfeld, F. Schmidl, and P. Seidel, Supercond. Sci. Technol. 19 200 (2006).
- (23) A. Grib, M. Mans, M. Buenfeld, J. Scherbel, F. Schmidl, and P. Seidel, Superconductive Electronics Conference (ISEC), 1, IEEE (2013).
- (24) R. Landauer, Phys. Today 31, 23 (1978);
- (25) H. B. Wang, S. Guenon, B. Gross, J. Yuan, Z. G. Jiang, Y. Y. Zhong, M. Grunzweig, A. Iishi, P. H. Wu, T. Hatano, D. Koelle, and R. Kleiner, Phys. Rev. Lett., 105, 057002 (2010).
- (26) Yu. N. Ovchinnikov and V. Kresin, Phys. Rev. Lett. 88, 214504 (2013).
- (27) M. Tachiki, K. Ivanovic, K. Kadowaki, and T. Koyama, Phys. Rev. B 83, 014508 (2011).
- (28) R.L. Kautz, J. Appl. Phys. 76, 5538 (1994).
- (29) S. Z. Lin, X. Hu, and L. Bulaevskii, Phys. Rev. B 84, 104501 (2011).
- (30) K. Y. Tsang, S. H. Strogatz, and K. Wiesenfeld, Phys. Rev. Lett. 66, 1094 (1991)
- (31) P. Hänggi, P. Talkner and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
- (32) M.I. Dykman, and M. A. Krivoglaz, JETP 50, 30 (1979).
- (33) R. Graham and T. Tél, Phys. Rev. A 31, 1109 (1985).
- (34) R. Graham and T. Tél, Phys. Rev. A 33, 1322 (1986).
- (35) M. I. Dykman, D. G. Luchinsky, R. Mannella, P. V. E. McClintock, N. D. Stein, and N. G. Stocks, Phys. Rev. E 49, 1198 (1994).
- (36) M. I. Dykman, B. Golding, L. I. McCann, V. N. Smelyanskiy, D. G. Luchinsky, R. Mannella, P. V. E. McClintock, CHAOS 11, 587 (2001).
- (37) G. Filatrella, S. Girotti, and G. Rotoli, Phys. Rev. B 75, 054510 (2007).
- (38) R.L. Kautz, J. Appl. Phys.52, 3528 (1981); Phys. Rev. A 38, 2066 (1988).
- (39) M.R. Kolachi, Yu.M. Shukrinov, M. Handipour, A.E. Botha, and M. Suzuki, Physica C 491, 63 (2013).
- (40) R.F. Fox, I. R. Gatland, R. Roy, and G. Vemuri, Phys Rev A, 38,5938 (1988).
- (41) D. E. Knuth, The art of Computer Programming Vol. 2 (Addison-Wesley, Reading, MA, 1969)
- (42) G. Filatrella, N. F. Pedersen, C. J. Lobb, and P. Barbara, Eur. Phys. J. B 34, 3 (2003).
- (43) Yu. M. Shukrinov, I. R. Rahmonov, and K. V. Kulikov, JETP Letters, 96, 588 (2012).
- (44) Yu. M. Shukrinov, P. Seidel, E. Il’ichev, W. Nawrocki, M. Grajcar, P. A. Plecenik, I. R. Rahmonov, and K. V. Kulikov, IOP Publishing, Journal of Physics: Conference Series 393 012020 (2012).
- (45) U. Welp, K. Kadowaki, and R. Kleiner, Nature Photonics 7, 702 (2013).
- (46) S. Sakai, P. Bodin, and N.F. Pedersen, J. Appl. Phys. 73, 2411 (1993).
- (47) S. Madsen, G. Filatrella, and N. F. Pedersen, Eur. Phys. J. B 40, 209 (2004).
- (48) S. Madsen, N. Grønbech-Jensen, N. F. Pedersen, and P. L. Christiansen, Phys. Rev. B 78, 174525 (2008).