J an 2 00 6 Statistics of transition times , phase diffusion and synchronization in periodically driven bistable systems

The statistics of transitions between the metastable states of a periodically driven bistable Brownian oscillator are investigated on the basis of a twostate description by means of a master equation with time-dependent rates. The theoretical results are compared with extensive numerical simulations of the Langevin equation for a sinusoidal driving force. Very good agreement is achieved both for the counting statistics of the number of transitions per period and the residence time distribution of the process in either state. The counting statistics corroborate in a consistent way the interpretation of stochastic resonance as a synchronization phenomenon for a properly defined generalized Rice phase. PACS numbers: 02.50.Ga, 05.40.-a, 82.20.Uv, 87.16.Xa Submitted to: New J. Phys. Transition times and phase diffusion 2


Introduction
Time-dependent systems that are in contact with an environment represent an important class of non-equilibrium systems. In these systems effects may be observed that cannot occur in thermal equilibrium, as for example noise-sustained signal amplification by means of stochastic resonance [1], or the rectification of noise and the appearance of directed transport in ratchet-like systems [2]. In the case of stochastic resonance the transitions between two metastable states appear to be synchronized with an external, often periodic signal that acts as a forcing on the considered system [3,4]. This force alone, however, may be much too small to drive the system from one state to the other. Responsible for the occurrence of transitions are random forces imposed by the environment in form of thermal noise. This kind of processes can conveniently be modeled by means of Langevin equations which are the equations of motion for the relevant degree (or degrees) of freedom of the considered system and which comprise the influence of the environment in the form of dissipative and random forces [5]. Unfortunately, neither the Langevin equations nor the equivalent Fokker-Planck equations for the time evolution of the system's probability density can be solved analytically for other than a few rather special situations [6]. Therefore, most of the prior investigations of stochastic resonance are based on the numerical simulation of Langevin equations [1], or the numerical solution of Fokker-Planck equations either by means of continued fractions [6] or Galerkin-type approximation schemes [7,8]. Analytical results have been obtained in limiting cases like for weak driving forces [1], or weak noise [7].
Interestingly enough, the same qualitative features that had been known from numerical investigations also resulted from very simple discrete models [9,10] that contain only two states representing the metastable states of the continuous model. The dynamics of the discrete states is Markovian and therefore governed by a master equation. The external driving results in a time-dependent modulation of the transition rates of this master equation. In recent work [11] it was shown that such master equations correctly describe the relevant aspects of the long-time dynamics of continuously driven systems with metastable states provided that the intra-well dynamics is fast compared to both the driving and the noise-activated transitions between the metastable states.
In the present work we will pursue ideas [4] suggesting that a periodically driven bistable system can be characterized in terms of a conveniently defined phase that may be locked with the phase of the external force, yielding in this way a more precise notion of synchronization for such systems. Various possible definitions of phases for a bistable Brownian oscillator have been compared in Ref. [4] with the main result that the precise definition does not matter much in the rate-limited regime where the master equation applies. Therefore, one will use the definition for which it is simplest to determine the corresponding phase. If we think of the system as of a particle moving in a bistable potential under the influence of a stochastic and a driving force then the Rice phase provides a convenient definition under the condition that the particle has a continuous velocity [4,12]. This phase counts the number of crossings of the potential maximum separating the two wells. The necessary continuity of the velocity is guaranteed if the hypothetical particle has a mass and obeys a Newtonian equation of motion. The continuity of the velocity is lost if the particle is overdamped and described by an "Aristotelian equation of motion," i.e. by a first-order differential equation for the position driven by a Gaussian white stochastic force. Then the trajectories are known to be nowhere differentiable and to possess further level crossings in the close temporal neighborhood of each single one.
We will allow for this jittering character of Brownian trajectories by setting two different thresholds on the both sides of the potential maximum and only count alternating crossings. By definition, the generalized Rice phase increases by π upon each counted crossing of either threshold. We will choose the threshold positions in such a way that an immediate recrossing of the barrier from either threshold is highly unlikely. Put differently, up to exceptional cases, the particle will move from the threshold to the adjacent well rather than to immediately jump back over the barrier to the original well. In the two-state description of the process we then can identify the phase by counting the number of transitions between the two states. With each transition the phase grows by π. In previous work [13] one of us gave explicit expressions for various quantities characterizing the alternating point process, comprised of the entrance times into the states of a Markovian two-state process with periodic time-dependent rates. The definition of the Rice phase in that work was based on the statistics of the entrance times into one particular state. Some of the statistical details then depend on the choice of this state. To avoid this ambiguity we here base our definition on the transitions between both states.
A further quantity that has been introduced in order to characterize stochastic resonance is the distribution of residence times [14]. It is also based on the statistics of the transition times and characterizes the duration between two neighboring transitions.
The paper is organized as follows: In Section 2 we introduce the model, formulate the respective equivalent Langevin and Fokker-Planck equations and give the formal solutions of the two-state Markov process with time-dependent rates that follow from the Fokker-Planck equation. In Section 3 several statistical tools from the theory of point processes are introduced by means of which the switching behavior resulting from the master equation can be characterized. Explicit expressions for various quantitative measures in terms of the transition rates and solutions of the master equation are derived. These are the first two moments of the counting statistics of the transitions from which the growth rate of the phase characterizing frequency synchronization, its diffusion constant and the Fano factor [15] follow. Further, the probabilities for n transitions within one period of the external driving force and an explicit expression for the residence time distribution in either state are determined. In Section 4 these quantities are estimated from stochastic simulations of the Langevin equations and compared with the results of the two-state theory. The paper ends with a discussion.

Rate description of a Fokker-Planck process
We consider the archetypical model of an overdamped Brownian particle moving in a symmetric double well potential U(x) = x 4 /4 − x 2 /2 under the influence of a timeperiodic driving force F (t) = A sin(Ωt). Throughout the paper we use dimensionless units: mass-weighted positions are given in units of the distance between the local maximum of U(x) and either of the two minima. The unit of time is chosen as the inverse of the so-called barrier frequency ω 0 = −d 2 U(0)/dx 2 = 1. The particle's dynamics then can be described by the Langevin equatioṅ where ξ(t) is Gaussian white noise, with ξ(t) = 0 and ξ(t)ξ(s) = δ(t − s). In the units chosen the inverse noise strength β equals four times the barrier height of the static potential divided by the thermal energy of the environment of the Brownian particle. The equivalent time-dependent Fokker-Planck equation describes the time-evolution of the probability density ρ(x, t) for finding the process at time t at the position x where the time-dependent potential V (x, t) contains the combined effect of U(x) and of the external time-dependent force F (t), i.e., We shall restrict ourselves to forces with amplitudes A being small enough such that V (x, t) has three distinct extrema x 1 (t), x 2 (t) and x b (t) for any time t, where x 1 (t) and x 2 (t) are the positions of the left and the right minimum, respectively, and x b (t) that of the barrier top. The barrier heights as seen from the two minima are denoted by V α (t) = V (x b (t), t) − V (x α (t), t) for α = 1, 2. We further assume that the particle only rarely crosses this barrier. This will be the case if the minimal barrier height V m = min t V α (t) is still large enough such that βV m > 4.5. Under this condition the time-scales of the inter-well and the intra-well motion are widely separated and the longtime dynamics is governed by a Markovian two-state process where the states α = 1, 2 represent the metastable states of the continuous process located at the minima at x 1 (t) and x 2 (t) [11]. The transition rates r α,α ′ (t), i.e. the transition probabilities from the state α ′ to the state α per unit time, are time-dependent due to the presence of the external driving. Explicit expression for the rates are known if the driving frequency is small compared to the relaxational frequencies ω α (t), α = 1, 2 and ω b (t) that are defined by In the limit of high barriers these rates take the form of Kramers' rates for the instantaneous potential [16], i.e.
A more detailed discussion of the validity of these rates in particular with respect to the necessary time-scale separation is given in Ref. [11]. More precise rate expressions that contain finite-barrier corrections are known [16] but will not be used here.
In the semi-adiabatic regime both the time-dependent rates and the driving force are much slower than the relaxational frequencies but no relation between the driving frequency and the rates is imposed [11,17]. Then, the long-time dynamics of the continuous process x(t) can be reduced to a Markovian two-state process z(t) governed by the master equation with the instantaneous rates (5) where p α (t), α = 1, 2 denotes the probability that z(t) = α. These equations can be solved with appropriate initial conditions to yield the following expressions for the conditional probabilities p(α, t | α ′ , s) to find the particle in the metastable state α at time t provided that it was at α ′ at the earlier time s ≤ t : where If the conditions are shifted to the remote past, i.e. for s → −∞ they become effectless at finite observation times t. The asymptotic probabilities then read We note that the asymptotic probabilities are periodic functions of time with the period T = 2π/Ω of the driving force, i.e.
p as Together with the conditional probabilities in eq. (7) the asymptotic probabilities p as α (t) describe the switching behavior of the process x(t) at long times.
Finally, we introduce the conditional probabilities P α (t | s) that the process stays uninterruptedly in the state α during the time interval [s, t), provided that z(s) = α. They coincide with the waiting time distribution in the two states which can be expressed in terms of the transition rates

Statistics of transition times and measures of synchronization
The two-state process z(t) is completely determined by those times at which it switches from state 2 into state 1, and vice versa from 1 to 2. These events constitute an alternating point process . . . t n < t * n < t n+1 < t * n+1 < . . . consisting of the two point processes t n (2 → 1) and t * n (1 → 2), see Ref. [18]. This alternating point process can be characterized by a hierarchy of multi-time joint distribution functions. Of these, we will mainly consider the single-and two-time distribution functions. The single-time distribution function W (t) gives the averaged number of transitions between the two states 1 and 2 in the time interval [t, t + dt), i.e., W (t)dt = #{s | s = t n or s = t * n , s ∈ [t, t + dt)} . This density of transition times is given by the sum of the entrance time densities W α (t) specifying the number of entrances into the individual state α The entrance time densities can be expressed in terms of the transition rates r α,α ′ (t) and the single-time probabilities p α (t), see Ref. [13] W 1 (t) = r 1,2 (t)p 2 (t) and W 2 (t) = r 2,1 (t)p 1 (t).
The density W (t) determines the average number of transitions N(t, s) between the two states within the time interval (t, s) In the limiting case described by the asymptotic probability p as α (t), see eq. (9), the entrance time densities W (t) also is a periodic function of time. Then, the average N(t, s) becomes periodic with respect to a joint shift of both times t and s by a period T and, moreover, the number N(s + nT, s) is independent of the starting time s and grows proportionally to the number of periods n.
There are two two-time transition distribution functions f (2) (t, s) and Q (2) (t, s) that specify the average product of numbers of transitions within the infinitesimal interval [s, s + ds) and the number of transitions in the later interval [t, t + dt). These two functions differ by the behavior of the process between the two times s and t. For the distribution function f (2) (t, s) the process may have any number of transitions within the time interval [s + ds, t). It is given by the sum of the two-time entrance distribution functions f α,α ′ (t, s) that determine the densities of pairs of entrances into the states α and α ′ at the prescribed respective times t and s. For t > s they can be expressed by the transition rates and the conditional probability p(ᾱ, t | α ′ , s), see Ref. [13] where the bar over a state,ᾱ, denotes the alternative of this state, i.e.1 = 2 and2 = 1.
Note that the conditional probability p(ᾱ, t | α ′ , s) allows for all possible realizations of the two-state process starting with z(s) = α ′ at time s up to the time t with any number of transitions in between. For t < s the function f α,α ′ (t, s) follows from the symmetry In the second two-time distribution function Q (2) (t, s) transitions at the times s and t are taken into account only if no further transitions occur between the prescribed times. It is again given by a sum of respective two-time distribution functions for the individual transitions from 1 to 2, and vice versa, and hence reads where the two-time entrance distribution function gives the density of an entrance into the stateᾱ at time s and an entrance into state α at time t conditioned upon processes z(t ′ ) =ᾱ, s < t ′ < t that stay constant in the time between s and t. Therefore, these distribution functions depend on the waiting time distribution in the stateᾱ as given by eq. (12), see also Ref. [13].
According to the theory of point processes, see Ref. [19], the second moment of the number of transitions within the time interval [s, t) results from the two-time distribution function f (2) (t, s) as Subtracting the squared average number N(t, s) 2 one obtains the second moment of the number fluctuations (δN(s + τ, s)) 2 . It is given by an analogous expression as N 2 (t, s) where If the time difference t − s becomes of the order of the maximal inverse rate, max t r 1,2 (t) −1 , i.e. of the order of the time scale on which the process becomes asymptotically periodic, the two-time distribution function factorizes into the product W (t)W (s) and g(t, s) vanishes for large t − s. Consequently, the double integral on the right hand side of the eq. (22) grows linearly with t in the asymptotic limit (t − s) → ∞. This leads to a diffusion of the transition number fluctuations, i.e. an asymptotically linear growth that can be characterized by a diffusion constant: which is a periodic function of the initial time s with the period T of the external driving force. This time-dependence reflects the non-stationarity of the underlying process.
It provides a quantitative measure of the relative number fluctuations and it assumes the value F (s) = 1 in the case of a Poisson process. Here, it is a periodic function of s which attains minima at the same time s as D(s). As already mentioned in the introduction, the number of transitions between the two states can yet be given another interpretation as a generalized Rice phase of the random process x(t). At each time instant the process has switched to another state, the phase grows by π. The simplest definition would be the linear relation Φ(t, s) = πN(t, s) where the phase is set to zero at the initial time t = s. With this definition the phase changes stepwise. A linear interpolation would lead to a continuously varying phase, but will not be considered here. Independently from its precise definition, in the asymptotic periodic regime, both the average and the variance of the phase increase linearly in time, superimposed by a modulation with the period of the driving force.
A further coarse-graining of the considered periodically driven process can be obtained by considering the number of transitions during an interval [s, t). By p α (n; t, s) we denote the probability that the process assumes the value z(s) = α at time s and undergoes n transitions up to the time instant t > s. Keeping in mind the significance of the waiting time distributions P α (t | s), given by eq. (7), as the probability of staying uninterruptedly in the state α and of the transition rates r α,β (t) as the probability per unit time of a transition from β to α one finds the following explicit expressions for the first few values of n invoking the basic rules of probability theory: where the states α andᾱ are opposite to each other. The probabilities with values of n > 0 can be determined recursively from the following set of differential equations: ∂ ∂t p α (n + 1; t, s) = − r α n+1 ,αn (t)p α (n + 1; t, s) + r αn,α n−1 (t)p α (n; t, s) p α (n + 1; s, s) = 0 where α n = α for n even α for n odd.
The hierarchy starts at n = 0 with p α (0; t, s) as defined in eq. (28). Accordingly, the probability P (n; t, s) of n transitions within the time interval [s, t) consists in the sum of the individual probabilities p α (n, t, s) P (n; t, s) = p 1 (n; t, s) + p 2 (n; t, s).
Finally, we like to emphasize that with the information at hand also the distributions of residence times can be determined. The residence times, say of the state α, are defined as the duration of the subsequent episodes in which the process dwells in state α without interruption. For a non-stationary process these times must not be confused with the life times of this state. Rather, the distribution of the residence times R α (τ ) is the life-time distribution in the state α averaged over the starting time with the properly normalized entrance density of the considered state [20,21], i.e., Here we assumed that the process has reached its asymptotic periodic behavior and expressed the entrance time distribution according to eq. (14). The denominator guarantees for the proper normalization of the residence time distribution.

Simulations
We performed numerical simulations of the Langevin equations by means of an Euler algorithm with a step size of 10 −3 and with a random number generator described in Ref. [22]. The comparison with different random number generators [23] did not produce any sensible differences. We always started the simulations at the left minimum x = −1 at time t = 0 and determined the first switching time t 1 as the first instant of time at which the value x = 1/2 was reached. The switching time t * 1 to the left well was defined as the first time larger than t 1 at which the opposite threshold at x = −1/2 was reached. For t 2 we waited until the positive threshold at x = 1/2 was again crossed, and so on. In this way a series of switching times t 1 , t * 1 , t 1 , t * 2 , t 3 , . . . was generated. We stopped the simulations after the time at which either 10 4 , or at least 10 3 periods of the driving force and at least 500 transitions to the right hand side of the potential had occurred. The left and right thresholds were taken in such a way that any multiple counting due to the inherent irregularities of the Langevin trajectories was excluded. Because of the time-scale separation between the inter-and intra-well dynamics the precise value of the thresholds at x = ±1/2 is immaterial as long as they stay in a finite distance from the top of the barrier. Then a fast return of a trajectory to the previously occupied metastable state after a crossing of the threshold can safely be excluded.
In the simulations the amplitude of the force was always A = 0.1. For the frequency we considered the two values Ω = 10 −3 and Ω = 10 −4 and the inverse temperature β was varied from 20 to 55 in integer steps. Note again that the barrier height of the reference potential U(x) is β/4 in units of the thermal energy. In Fig. 1 up to a time t as a function of t are depicted for different inverse temperatures together with averages estimated from the full trajectory and those averages following from the two-state model, see eq. (15). The average growth behavior is determined by the average number of transitions per period in the asymptotic state, N(T + s, s) as , which is independent of s. In Fig. 2 the estimated value of this number is compared with the two-state model prediction as a function of β. The agreement between theory and simulations is excellent even for the relatively high temperature values and corresponding low transition barriers for the values of β < 25. The average number of transitions monotonically decreases with falling temperature. At the temperature where N(T + s, s) as = 2 the system is optimally synchronized with the driving force in the sense that the Rice phase increases on average by 2π per period of the driving force. This optimal temperature depends on the driving frequency and becomes lowered for slower driving. In the vicinity of the optimal temperature the decrease of N(T + s, s) as with inverse temperature is smallest. The emerging flat region resembles the locking of a driven nonlinear oscillator and becomes more pronounced for smaller driving frequencies.
The fluctuations of the number of transitions δN(t) ≡ N(t, 0) − N(t, 0) for times up to t = 500T in a single simulation together with the theoretical average following from eq. (22) are shown in Fig. 3. The average behavior is characterized by the number fluctuations per period δN 2 (T, 0). This quantity was estimated from the simulations as a function of the inverse temperature. In Fig. 4 we compare the prediction of the two-state model, see eq. (22) with numerics. These number fluctuations exhibit a local minimum very close to the optimal temperature where two transitions per period occur on average, see Figs. 2 and 4. This minimum is more pronounced for the lower driving frequency. The fluctuations of the Rice phase also assume a minimum at this optimal temperature. This means that the phase diffusion is minimal at this temperature.
For the Fano factor F (0) we obtain an absolute minimum near the optimal  temperature, see Fig. 5. For higher temperatures the Fano factor may become larger than one, whereas it approaches the Poissonian value F = 1 for low temperatures because then, transitions become very rare and independent from each other. Also here, theory and simulations agree indeed very well.
Next, we consider the probabilities P (n) = P (n; T, 0) for finding N(T, 0) = n transitions per period in the asymptotic, periodic limit. For this purpose we count the number k of transitions within each period [nT, (n + 1)T ) and determine their relative frequency occurring in a simulation. A comparison with the prediction of the twostate model determined from the numerical integration of eq. (28) and from eq. (30) is collected in Fig. 6 for different temperature values. The agreement between simulations and theory is within the expected statistical accuracy. For large temperatures there is a rather broad distribution of n-values around a most probable value n * . With decreasing temperature, the most probable value moves to smaller numbers whereby the width of the distribution shrinks. Once n * = 2 the probability P (2) further increases at the cost of the other k values with decreasing temperature until k = 0 gains the full weight in the limit of low temperatures.   The symbols are the same as in Fig. 2, the solid lines display equation (22). Note that δN 2 (T ) is proportional to the diffusion constant of the Rice phase. The vertical thin black lines indicate the optimal inverse temperatures as found from the synchronization of the averaged phase (see Fig. 2). They coincide remarkably well with the positions of the minima of δN 2 (T ) . So we find that at the optimal temperature also the fluctuations of the number of transitions and consequently of the generalized Rice phase are suppressed.  Fig. 2. The minimum positions of the Fano factor again nicely coincide with the optimal temperature values indicated by the thin black lines. With decreasing frequency these minima move to lower temperatures and become broader and deeper. For sufficiently low temperatures the transitions become very rare and almost Poissonian leading to the asymptotic low temperature limit F = 1.  Figure 6. The probabilities P (n) = P (n; T, 0) are shown for the driving the driving strength A = 0.1 and the driving frequency Ω = 10 −4 ; for two temperatures that are higher than the optimal one corresponding to β ≈ 40 in panel (a) and and for the optimal and a lower temperature in panel (b). As one expects for a good synchronization of the system with the driving force at the optimal temperature one finds two transitions per period with about 90% of the total weight. For high temperatures the distribution becomes rather broad with a maximum at some large number of transitions. On the contrary, for low temperatures the probability is largest at n = 0 and decreases with n. The agreement of the two-state theory resulting from equation (30) with the simulations is remarkably good.
The degree of synchronization of the continuous bistable dynamics with the external driving force can be characterized by the value of the probability P (2). As a function of the inverse temperature β it has a maximum close to the corresponding optimal temperature, see Fig. 7. For the longer driving period the maximum is at a lower temperature and its value is higher. Fig. 8 depicts P (2) as a function of the frequency at a fixed temperature. It has a maximum at some optimal value of the frequency.
Finally, we come to the residence time distributions which can be estimated from histograms of the simulated data. In the Fig. 9 these histograms are compared with the results of the two-state model. Theory and simulations are in good agreement and show a transition from a multi-modal distribution with peaks at odd multiples of half the period to a mono-modal distribution at zero as one would expect. At stochastic resonance taking place at the inverse temperature β ≈ 40 the distribution is also monomodal with its maximum close to half the period.

Discussion and conclusions
In this work we investigated the statistics of transitions between the two states of a periodically driven overdamped bistable Brownian oscillator. We compared results from extensive numerical simulations with theoretical predictions of a discrete Markovian two-state model with time-dependent rates. For the studied slow driving frequencies the rates can be taken as the adiabatic, or frozen, rates. Although the minimal barriers separating the states from each other may become rather small measured in units of the thermal energy (βV min ≈ 4 for the lowest inverse temperature β = 20) we did not take into account finite barrier corrections for the rates and still obtained very good for two transitions within one period is shown as a function of the inverse temperature β. Results from the simulations are depicted by red (×) and blue (+) crosses for two frequencies. They agree well with the respective theoretical predictions of the two-state theory resulting from equation (30) (solid lines). The probabilities show the typical stochastic resonance behavior with a maximum very close to the optimal temperature. As for the minimum of the Fano factor, this maximum is more pronounced at the lower driving frequency.  (30)) for two transitions within one period is shown as a function of the frequency at the inverse temperature β = 35. It also shows a resonance like behavior with a maximum at an optimal frequency.  Figure 9. Residence time distributions for either of the two states are shown for Ω = 10 −4 and different values of the temperature. Theoretical results obtained from (eq. 31) are displayed as solid lines, estimates from the simulations as symbols. In panel (a) R(τ ) is shown for higher temperatures and in panel (b) for the optimal and lower than optimal temperatures. The peak structure and its temperature dependence is in complete agreement with previous findings, see Ref. [1]. agreement for all considered quantities even at the highest chosen temperatures. We considered different quantities from the counting statistics of transitions such as the average and variance, as well as the probability of a given number of transitions in one period. These quantities are directly related to the average Rice phase, its diffusion and the probability of a given increase of this phase. We note that, apart from the average phase, these quantities depend on the position of the chosen period for which we took s = 0 as the starting time. A different, more sophisticated, way of analyzing the counting statistics would have been to position the averaging window of the length of a period such that the variance of the number of transitions is minimal. But already with the present simpler method we could detect significant signatures of synchronization in a noisy nonlinear system such as the locking of frequency and phase which both occur at basically the same parameter values.
The novel analytical expression for the residence time distribution derived for the two-state model also agrees very well with the simulation results.