# Continuous and discontinuous phase transitions and partial synchronization in stochastic three-state oscillators

## Abstract

We investigate both continuous (second-order) and discontinuous (first-order) transitions to macroscopic synchronization within a single class of discrete, stochastic (globally) phase-coupled oscillators. We provide analytical and numerical evidence that the continuity of the transition depends on the coupling coefficients and, in some nonuniform populations, on the degree of quenched disorder. Hence, in a relatively simple setting this class of models exhibits the qualitative behaviors characteristic of a variety of considerably more complicated models. In addition, we study the microscopic basis of synchronization above threshold and detail the counterintuitive subtleties relating measurements of time averaged frequencies and mean field oscillations. Most notably, we observe a state of suprathreshold partial synchronization in which time-averaged frequency measurements from individual oscillators do not correspond to the frequency of macroscopic oscillations observed in the population.

## 1Introduction

A great number of physical systems consist of individual entities with periodic, or nearly periodic, dynamics. Ranging from collections of chemical consitutents to groups of social entities – for example, applauding individuals whose clapping is repetitive – these systems serve as a battleground of sorts for the competition between the dynamics of individual constituents and the large scale cooperation favored in many cases by the nature of their mutual interactions. Owing to the ubiquity and certainly, in part, to the dramatic nature of the emergent synchronized behavior in such naturally oscillating settings, the subject has been intensely studied in the physics literature for several decades, with the Kuramoto oscillator and its kin serving as prototypical models on which many studies are based [1].

Recently, a simple class of models of macroscopic synchronization have provided a number of additional insights into the large-scale phenomena ocurring in noisy discrete coupled oscillators, including detailed characterizations of both the universal critical behavior of the continous phase transition [5] as well as the effects of spatial disorder in such populations [7]. While retaining many qualitative characteristics of more complex models, the discrete oscillators remain sufficiently simple to provide results unattainable in most of the paradigmatic settings.

In this paper, we generalize our class of stochastic, discrete oscillator models and detail its use in a variety of new contexts. By generalizing the form of the inter-oscillator coupling, we show that our class of mean field models encompasses oscillators which can undergo either supercritical or subcritical Hopf bifurcations, depending on the microscopic specifics of the coupling. In addition, we study dichotomously disordered populations of oscillators and show that the bifurcation can be either supercritical or subcritical depending on the degree of disorder in the population. Such behaviors are reminiscent of a number of significantly more complex oscillator models [8], including Daido’s generalized Kuramoto oscillators [15], where either disoder or microscopic coupling specifics can alter the nature (continuous or discontinuous) of the transition. However, our model provides a far simpler setting for observing both continous and discontinous transitions to synchronization.

In addition, we study the microscopic underpinnings of synchronization above threshold. In particular, we look at time-averaged frequency and its relationship to phase synchronization above threshold (which turns out to be interestingly counterintuitive). Again, we do this for a specific model from our general class (for both single and dichotomously disordered populations), but we expect the results to hold for our entire class of models undergoing a Hopf bifurcation. This is somewhat similar to the partial synchronization seen in other models [18], but again, our model is simpler and, perhaps, more transparent.

In Section 2 we present our model and highlight its essential parameters. Here we note that for systems of oscillators with identical transition rates between states the control parameter for the phase transition is the coupling strength among oscillators; when the array includes oscillators with different transition rates, the degree of disorder is also a control parameter. In Section 3 we show that depending on the values of microscopic parameters, this model can exhibit both subcritical (or first-order) and supercritical (or second-order) phase transitions, as a function of the coupling strength and also as a function of the degree of disorder. Section 4 deals with the microscopic underpinnings of the synchronization phenomenon and the connection between phase synchronization and frequency entrainment in our system. Section 5 presents a summary and further discussion of our results.

## 2The Model

In this section we present our model in some detail, repeating some of our early presentations [5] because of important (albeit simple) generalizations of the model. We begin by considering a stochastic three-state model governed by transition rates (see Figure 1), where each state may be interpreted as a discrete phase [5]. Because the transitions among states are unidirectional and do not conform to deterministic rate laws, the model retains a qualitative link with a noisy phase oscillator.

The linear evolution equation of a single oscillator is , where the components of the column vector ( denotes the transpose) are the probabilities of being in states at time , where and

The system reaches a steady state for . The oscillator’s periodicity, as contained in the timescale of the cycle , is determined by ; that is, the time evolution of our simple model qualitatively resembles that of the discretized phase of a generic noisy oscillator with the intrinsic eigenfrequency set by the value of .

To study interacting arrays of these oscillators, we couple individual units by allowing the transition rates of each unit to depend on the states of the units to which it is connected. Specifically, for identical units we choose the transition rate of a unit from state to state as

where and when , is the coupling parameter, is the transition rate parameter, is the number of oscillators to which unit is coupled, and is the number of units among the that are in state . We introduce the real constants , , and to encompass in a general way our two previous coupling functions [5]. Each unit may thus transition to the state ahead or remain in its current state, and the propensity for such a change depends on the states of its nearest neighbors. In our earlier works we considered the globally coupled system, , and also nearest neighbor coupling in square, cubic, or hypercubic arrays, ( dimensionality). Here we focus on globally coupled arrays.

For a population of identical units in the mean field (globally coupled) version of this model we can replace with the probability , thereby arriving at a nonlinear equation for the mean field probability, , with

Normalization allows us to eliminate and obtain a closed set of equations for and . We can then linearize about the fixed point , yielding a Jacobian with a set of complex conjugate eigenvalues which determine the stability of this asynchronous state. Specifically, we find that

where is a nonzero constant for all finite , , and and we introduce the abbreviation . The eigenvalues cross the imaginary axis at , yielding

with

For and (that is, ), represents a Hopf bifurcation point, indicating the emergence of macroscopic oscillations indicative of synchronization. Furthermore, we require that to ensure the bifurcation happens at a positive value of . We note that in previous studies we have used , yielding and [7], and , yielding and [5]. In addition, we stress that while a range of models may prove useful for exploring the phase transition behavior near threshold (see, for example, [5]), only models with provide physically appealing characteristics far above threshold (see, for example, [7]). Specifically, only for does the frequency of a perfectly synchronized set of oscillators maintain a nonzero finite value (). Below we explore in more detail the nature of the Hopf bifurcation associated with the class of models described by the permitted values .

In addition to the single population case, we consider globally coupled arrays of oscillators that can have one of different transition rate parameters, , . As detailed in [7], the probability vector is now -dimensional, , and the added subscript on the components of keeps track of the transition rate parameter. Explicitly, the component is the probability that a unit with transition rate parameter is in state . The evolution of the probability vector is given by the set of coupled nonlinear differential equations , with

Here

and

The function is the fraction of units which have a transition rate parameter .

Because it closely appeals to physical intuition [7] for oscillators far above threshold, we limit ourselves to the case for dichotomously disordered oscillators. We further limit our focus here to uniform distributions , but note that relaxing this constraint has been shown to preserve the qualitative features of the model [7]. For uniform distributions and , probability normalization again allows us to reduce this to a system of coupled ordinary differential equations. We can then linearize about the disordered state and arrive at a Jacobian parameterized by a collection of transition rate parameters and a coupling strength .

While it has been shown that the qualitative essence of the model remains similar for and even for completely disordered populations [7], we focus here only on the simple dichotomously disordered case, . As shown in [7], the four eigenvalues of the corresponding Jacobian are given by

where

and

Aside from an overall factor , Eqs. (Equation 8) depend only on the relative width variable , and therefore the critical coupling , that is, the value of at which Re , depends only on . As Re does not vanish for any , corresponds to a Hopf bifurcation, and our model exhibits macroscopic oscillations indicative of large-scale cooperation. We note that increases with increasing , indicating that a stronger coupling is necessary to overcome increasingly different values of and .

In what follows, we make use of the synchrony order paramter to characterize the emergence of phase synchrony [17]. This parameter is defined as

Here is a discrete phase, taken to be for state at site . The brackets represent an average over time in the steady state and an average over all independent trials. Therefore, serves as a measure of phase synchronization.

## 3Continuous and Discontinuous Transitions to Synchrony

In the mean field limit, the order of the phase transition to synchrony is closely tied to the nature of the Hopf bifurcation. Specifically, a subcritical Hopf bifurcation corresponds to a discontinuous (sometimes called first-order) phase transition, while a supercritical Hopf bifurcation indicates a continuous (second-order) transition. As such, we place special emphasis in this manuscript on the sign of , the first Lyapunov coefficient, which provides information on the nature of the Hopf bifurcation and, by extension, on the order of the phase transition.

In general, can be calculated using the projection technique given in [16], which relies on a multivariate Taylor expansion of the vector field describing the dynamics in question about an equilibrium point. For a general -dimensional dynamical system with an equilibrium point undergoing a Hopf bifurcation at parameter value , is given by [16]

where is the typical complex scalar product, is the identity matrix, and and are right and left eigenvectors of the Jacobian given by

Furthermore, is chosen so that , and and are multilinear, -dimensional vector functions corresponding to the lowest order nonlinear coefficients in the Taylor expansion of the vector field. That is,

with indicating the equilibrium point of the vector field around which we expand and the bifurcation parameter, , evaluated at the bifurcation point.

### 3.1Continuous and Discontinuous Transitions in a Single Population of Identical Oscillators

For the case of a single population of oscillators described by Eqs. (Equation 2) and (Equation 3), can be analytically calculated using the technique outlined above. Specifically, we set (without loss of generality) and consider the equilibrium point at and find and to be

independent of , , and . Then, calculating the multivariable functions and with Eq. (Equation 13) and using as defined in Eq. (Equation 4) along with Eqs. (Equation 11) and (Equation 14), we find after simplification that

As a result, the nature of the Hopf bifurcation depends on the choices , , and . Specifically, if we assume , we have

A similar result holds for , but we shall here restrict ourselves to the intuitively reasonable models positing and ; that is, the oscillators one state ahead of the one in question can only increase (or not affect) the transition rate and those behind can only decrease (or not affect) the transition rate. To verify these predictions, we show numerical solutions to the mean field equations in Figure 2; the top panel represents an example in the subcritical regime () while the bottom panel shows an example in the supercritical regime (). A clear distinction can be made in the neighborhood of the critical point. We also note that the continuous transition is characterized by the classical mean field exponent .

We further observe that the choice leads to , indicating a supercritical Hopf bifurcation and rendering the model applicable to studies of continuous phase transitions [5]. With universality in mind, we stress that any choice of parameters yielding a supercritical bifurcation should show similar critical behavior. On the other hand, the choice , while physically appealing above threshold, falls at a singular point separating the subcritical and supercritical cases (). The flexibility inherent in the choice of coefficients , , and speaks to the richness of our generic three-state oscillator and highlights its utility in studying synchronization in both supercritical (see [5]) and subcritical regimes. We proceed to study the model in the presence of dichotomous disorder and show that, for a given choice , the level of disorder can also alter the nature of the Hopf bifurcation and hence the order of the phase transition. We select the physically appealing choice and, while the behavior near the critical point may depend in some sense on this choice, we stress that our overarching goal remains unhindered. That is, we are able to provide an example indicating that disorder alone can affect the nature of the transition. The ubiquity of this phenomenon across the entire range of models remains an open question for future work, though we note that similar results are observed for all parameter choices mentioned in this paper.

### 3.2Continuous and Discontinuous Transitions in a Dichotomously Disordered Population

Interestingly, the dichotomously disordered system corresponding to Eqs. (Equation 5), (Equation 6), and (Equation 7) with and can undergo either a subcritical or supercritical bifurcation depending on the value of characterizing the individual transition rates. The transition to synchrony occurs at a single value of the coupling dependent on the relative width parameter [7]. As such, and are not truly independent parameters, and we can in principle eliminate and consider to be the bifurcation parameter of interest. Then, using the machinery of Eqs. (Equation 11), (Equation 12), and (Equation 13), it is a straightforward but tedious exercise to numerically evaluate the first Lyapunov coefficient corresponding to the Hopf bifurcation occurring at . As shown in Figure 3, the sign of varies depending on the relative width parameter (which in turn determines the critical coupling ). Hence, the phase transition to synchrony can appear continuous or discontinuous depending on the relative difference between the transition rate parameters in the two populations.

To verify these predictions, we solve the mean field equations numerically in both the subcritical () and supercritical () regimes. In the former case, we consider the case , . Figure 4 clearly indicates that the transition to synchrony is marked by a discontinuous change in the order parameter as eclipses . In addition, a small region of marked hysteresis appears just below threshold. Remarkably, this indicates that a stable disordered solution coexists with a stable, synchronized solution (the stable limit cycle) just before threshold.

By contrast, the case corresponds to a supercritical Hopf bifurcation reminiscent of a continuous phase transition. As shown in Figure 5, the transition is characterized by a continously increasing order parameter; no hysteresis is evident. We note also that the order parameter displays a power law increase near the onset of the bifurcation marked by the mean field critical exponent . This is expected both from the Hopf bifurcation theorem, which prescribes the dependence of the limit cycle radius (closely related to , the order parameter) near the onset of synchrony, and also because of the analogy with phase transitions in an infinite-dimensional space (see [5]).

Interestingly, these results indicate that the degree of spatial disorder may fundamentally alter the nature of the phase transition to synchrony. In both the subcritical and supercritical cases, synchronization is marked by the destabilization of the non-synchronous state at a single value of , giving rise to emergent oscillations in a macroscopic variable [for example, ]. Hence, both cases retain the qualitative features of synchronization in disordered populations discussed in our previous work [7]; however, the details of the onset of such cooperation distinguish the two cases.

## 4Microscopic Underpinnings of Synchronization

Having detailed two distinct mechanisms by which synchronization might arise, we now explore in detail the microscopic subtleties underlying synchronization above threshold. As detailed in [5] and mentioned above, synchronization occurs in the mean field limit via the destabilization of a nonsynchronous fixed point. Specifically, a single pair of complex conjugate eigenvalues corresponding to the linearized fixed point cross the imaginary axis at , giving rise to stable oscillations in the macroscopic variables characterizing the system [in our case, the components of ]. While the onset of this behavior is dependent on the choice and also the magnitude of disorder within the system (Sec. Section 3), the qualitative features of the synchronized state remain identical above threshold in both the subcritical and supercritical case. Hence, we limit our attention to several illustrative cases, but note that our results hold also for the supercritical case (and in fact then entire range of ). Specifically, in what follows, we take and consider a single population as well as a dichotomously disordered population with .

In particular, the threshold is marked by the onset of coherent temporal oscillations in the components of . We characterize the microscopic underpinnings of these oscillations by considering , the time-averaged frequency of oscillator in the steady state. We perform simulations on globally coupled lattices of units of a single population with and also of a dichotomously disordered population with and . As shown in Figs. Figure 6 and Figure 7, the distribution of frequencies clusters around the values prescribed by (or and for populations one and two, respectively) far below threshold (top panels). Specifically, for a deterministic oscillator with transition rate , is given by ; when (or and ), this gives the central peak of the histogram for the relevant population. We compare these histograms with the power spectra , where is the Fourier transform of . As threshold is eclipsed (middle panels), a peak arises in the power spectrum of the macroscopic variables , though the frequency of this peak does not correspond with the individual ’s of oscillators constituting the population. In the dichotomous case, this peak only roughly corresponds with the time averaged frequencies from population two and completely exceeds even the maximum characterizing population one. As is further increased, the descrepency between the time-averaged frequency histograms and the macroscopic oscillation frequencies decreases. In addition, in the disordered case, the histograms for the two populations become increasingly narrow and closer to one another (bottom panel). We note that as becomes tremendously large, the histograms become extremely narrow and begin to overlap at a frequency determined by the frequency of the macroscopic oscillations, as expected (indicative of perfect synchronization). Nonetheless, the behavior for finite, intermediate is rather counterintuitive and points to a rich microscopic dynamics underlying the cooperative behavior.

To further explore these trends, we consider the stochastic variable , the waiting time in a single state for an individual oscillator. represents the time the oscillator spends in a single state before transitioning to the subsequent state . For computational efficiency, we record for a representative subpopulation of units ( units from each population, one and two, in the disordered case). Figures Figure 8 and Figure 9 show histograms of the variable taken over this representative subpopulation once steady state was reached. Clearly, all relevant subpopulations consist of oscillators whose steps most often correspond to the frequency of the macroscopic oscillation (shown by the solid vertical line). That is, the peak of the histograms occur at a value comensurate with the frequency peak in the power spectrum of the components of . However, Figure 8 shows that the distribution of is bimodal, with a significant peak occuring at which downward biases the time-averaged frequencies of individual units. We note that as coupling increases significantly above threshold, the distribution becomes unimodal with a peak at corresponding to the frequency of macroscopic oscillation. In the disordered case, only population one, characterized by significantly lower time-averaged ’s, shows a bimodal distribution with a significant peak at . In fact, these long waiting times, while not the dominant macroscopic behavior, pervade the microscopic dynamics in such a way that the time-averaged frequencies become downward biased and no longer accurately represent the macroscopic dynamics. Interestingly, population two has become sufficiently synchronized that the second peak is effectively nonexistent, and thus the frequencies overlap more closely with the macroscopic “mean field” frequency. The right insets of Figs. Figure 8 and Figure 9 show histograms for single units chosen from the populations (or subpopulations). Again, the unit chosen from the single population case shows a bimodal distribution with significant “anamolous” peaks near . In the disordered case, the unit from population one shows a bimodal waiting time distribution characterized by occasional waiting times in the neighborhood of in addition to those corresponding to the macroscopic oscillations. Finally, in Figure 10 we show the time evolution of the subpopulations along with the macroscopic variable for each population. At any given time, the majority of oscillators in each population is synchronized, leading to the smooth oscillations of the macroscopic variable. However, isolated single units are prone to long waiting times, particularly in the less synchronized population (population one, left panel, in this example). These anamolously long waiting times, which serve to bias the time averaged frequencies of each individual unit, nevertheless do not substantially disrupt the macroscopic oscillations, largely because the occurance of coincident long waits is fairly uncommon. That is, the long waiting times do not appear in any significantly correlated way among individual constituents of the population.

## 5Discussion

We have shown that a class of simple, discrete models of stochastic phase coupled oscillators can undergo either a subcritical or supercritical bifurcation to macroscopic synchrony, depending on the chosen form of the microscopic coupling. As such, the different instances of the model can be used to study either continuous phase transitions [5] or discontinuous transitions exhibiting hysteresis, a characterstic seen in detailed theoretical models of, e.g., coupled Josephson junctions [9] but only observed in significantly more complex coupled oscillator models [15]. We stress that universality suggests that all models in this class exhibiting continous phase transitions should show similar behavior near the critical point, and this served as the basis of our eariler studies [5]. Nevertheless, it is remarkable that minor modifications in microscopic coupling can alter the nature of the bifurcation in such a fundamental way.

In addition, we have shown that in dichotomously disordered populations, both subcritical and supercritical Hopf bifurcations can occur, and the distinction is completely determined by the relative width characterizing the transition rate disorder between the two populations. While the qualitative features of the transitions within each class (subcritical and supercritical) appear identical, the distinction between classes points to fundamentally different mechanisms underlying the initial emergence of phase synchronization. In particular, it is striking that the level of disorder within a population, as measured by , can significantly alter the behavior near the critical point (though we stress that behavior even moderately above threshold is qualitatively indistinguishable).

Finally, we have studied the microscopic basis of phase synchronization above threshold. It is initially counterintuitive that phase synchronization, defined in terms of the Hopf bifurcation and temporal oscillations in the macroscopic variable (and measured in the order parameter ), is not contingent upon the existence of overlapping distributions of . That is, our results regarding the discrete oscillator model highlight the complexity of microscopic dynamics underlying macroscopic cooperation and point to a potentially misleading subtlety. Whereas phase synchronization is often considered a stronger condition than frequency entrainment – defined using an order parameter built upon the notion that a fraction of units display identical time-averaged frequencies in the oscillator population – we here report subtle microscopic features which distinguish the two without establishing a clear hierarchy. For example, Hong et al. [17] show that for disordered populations of Kuramoto oscillators, the lower critical dimension for frequency entrainment is lower than that for phase synchronization in locally coupled oscillators, indicating the relative ease with which frequency entrainment is achieved. They note that the two transitions coincide in the case of globally coupled units. Contrast that with our dichotomously disordered population, for which phase synchronization occurs without any overlap in the frequency distributions: that is, no oscillator from population one has the same frequency as any oscillator from population two. While a direct comparison is not plausible owing to the specific differences between models and order parameters, we stress that any order parameter related to time-averaged measurements of frequencies would, for our model, be misleading and provide potentially counterintuitive results. The emergence of a nonzero , which measures phase synchronization, corresponds with the loss of stability of the asynchronous fixed point (the Hopf bifurcation). This does not guarantee similar distributions of time-averaged frequencies in the two populations; in fact, we can readily see that synchronization occurs while the frequency distributions are entirely distinct. Furthermore, the frequency of the macroscopic oscillations of the mean field does not always coincide with the time-averaged frequencies of the oscillators constituting the population (or any subpopulation). Only when coupling is sufficiently large to substantially reduce the anamolously long waiting times which bias will the frequency distributions begin to overlap one another and coincide with the frequency of the mean field oscillations. Because these long waiting times appear more readily in the population with the smaller , the time-averaged frequencies of the two populations are disproportionately affected, meaning that the populations will appear to behave quite differently in terms of average frequency. This in fact underlies the stark differences in the degree of synchronization between two populations as measured by the order parameter , and provides an intuitive description capable of explaining this discrepancy. Our previous results show that completely disordered populations show qualitative similarities with the dichotomously disordered case [7]; hence, we are led to cautiously speculate that wholly disordered populations are also characterized by waiting times distributed with long tails, and hence time-averaged frequencies become downwardly biased, meaning that the order parameter for frequency entrainment, in the typical sense, will not accurately reflect the macroscopic cooperation. Further studies along these lines are currently in progress.

Finally, the results of this work raise the following question: how dependent is the above phenomenon on the choice of a discrete phase model? Would similarly counterintuitve results arise in continous phase oscillators? In fact, a recent study by Rosenblum and Pikovsky [18] suggests that a similar (though not identical) state of partial synchronization arises in continous oscillators coupled in a highly nonlinear fashion. Specifically, they find that in globally coupled oscillators, phases exist in which certain subpopulations are characterized by time-averaged frequencies which are not commensurate with the oscillations of the mean field, that is, they are not locked with the macroscopic oscillations induced in the population. While once again the differences between the models make direct comparison difficult, it is nevertheless clear that measurements of time-averaged frequencies provide potentially counterintuitive results, even in globally coupled arrays. In the case of our stochastic discrete oscillators, the behvavior is quite transparent once viewed in terms of , though it is not clear whether a similar mechanism underlies the phenomenon in the continous phase model. Uncovering the relationship between the superthreshold phase in our model and that in the continous oscillator model of [18] remains an open question, but even the superficial similarities between the results motivate continued efforts along these lines.

## Acknowledgments

This work was partially supported by the National Science Foundation under Grant No. PHY-0354937.

### References

- S. H. Strogatz,
*Nonlinear Dynamics and Chaos*(Westview Press, 1994). - A. T. Winfree, J. Theor. Biol.
**16**, 15 (1967). - Y. Kuramoto,
*Chemical Oscillations, Waves, and Turbulence*(Springer, Berlin, 1984). - S. H. Strogatz, Physica D
**143**, 1 (2000). - K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg, Phys. Rev. Lett.
**96**, 145701 (2006). - K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg, Phys. Rev. E
**74**, 031113 (2006). - K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg, Phys. Rev. E, 041107 (2007).
- F. Giannuzzi, D. Marinazzo, G. Nardulli, M. Pellicoro, and S. Stramaglia. Phys. Rev. E
**75**, 051104 (2007). - G. Filatrella, N. F. Pedersen, and K. Wiesenfeld. Phys. Rev. E
**75**, 017201 (2007). - H-A. Tanaka, A. J. Lichtenberg, and S. Oishi, Phys. Rev. Lett.
**78**, 2104 (1997). - J. A. Acebrón, L. L. Bonilla, and R. Spigler, Phys. Rev. E
**62**, 3437 (2000). - M. Y. Choi, H. J. Kim, D. Kim, and H. Hong, Phys. Rev. E
**61**, 371 (2000). - D. Pazó, Phys. Rev. E
**72**, 046211 (2005). - H. Hong, M. Y. Choi, B-G. Yoon, K. Park, and K-S Soh, J. Phys. A: Math. Gen.
**32**, L9 (1999). - H. Daido. Physica D
**91**, 24-66 (1996); Phys. Rev. Lett**73**, 760 (1994). - Yu. A. Kuznetsov,
*Elements of Applied Bifurcation Theory*, 2nd ed. (Springer, New York, 1998). - H. Hong, H. Park, and M. Choi, Phys. Rev. E
**71**, 054204 (2004); H. Hong, H. Park, and M. Choi, Phys. Rev. E**72**, 036217 (2005). - M. Rosenblum and A. Pikovski. Phys. Rev. Lett
**98**, 064101 (2007).