text:
# On the collective behavior of parametric oscillators

## I Introduction

Over the past two decades there has been an increasing interest in the nonequilibrium behavior of spatially extended systems modeled as ensembles of simple dynamical units coupled to each other. The collective evolution of such discrete coupled systems often exhibits qualitatively different behavior from that of the single units. An example of a system that has attracted a great deal of attention is a collection of a large number of coupled limit-cycle (phase) oscillators with randomly distributed natural frequencies . This system has been invoked as a simple model for coupled chemical, biological or physical oscillators. A most spectacular collective phenomenon discovered with this model is a synchronization phase transition involving mutual entrainment of the oscillators through frequency and phase locking. However, in this model the coupling is assumed to be sufficiently weak that the amplitude is not affected. As a result, the model cannot describe amplitude instabilities. A system that does exhibit a rich variety of amplitude instabilities consists of coupled parametric oscillators and is the subject of this paper .

There is a large literature on single parametric oscillators with periodically  or stochastically  modulated frequencies, perhaps also subject to thermal fluctuations . A wide variety of linear and nonlinear, deterministic and stochastic, systems that exhibit energetic instabilities can be modeled as simple parametric oscillators undergoing the amplitude instability known as “parametric resonance”. It is therefore rather surprising that there is actually little work on coupled parametric oscillators. Work previous to our own that we are aware of is limited to a global parametric modulation that acts on all the oscillators in exactly the same way. Examples include parametrically pumped electrons in a Penning trap , pattern formation under global resonant forcing , time-periodic loading of an elastic system , and nonlinear parametrically driven lattices .

In a previous paper  we introduced a very simple, analytically solvable model of coupled parametric oscillators consisting of an infinite set of globally coupled harmonic oscillators subject to time-periodic piecewise-constant modulations with randomly distributed quenched phases (“quenched” in this context means that the phase of the frequency modulation of each oscillator is set at time $`t=0`$ and then remains unchanged). We showed the appearance of a collective parametric instability: even though each individual oscillator is in its stable parameter domain, the average amplitude of the coupled system may diverge monotonically. This collective instability occurs when the amplitude of parametric modulation is sufficiently large that the instantaneous frequency of the oscillators temporarily becomes imaginary. The instability is re-entrant with respect to the strength of the coupling of the oscillators and persists in the overdamped limit. In the presence of a saturating nonlinearity , it generates a pitchfork bifurcation, corresponding to a genuine second-order nonequilibrium phase transition (implying the spontaneous breaking of spatial symmetry and ergodicity).

The purpose of this paper is to show that in addition to the instability described by monotonic growth of the mean amplitude, the globally coupled infinite system can also undergo transitions to a collective oscillatory instability with an intrinsic frequency that is not connected in a simple way to either the frequencies of the individual oscillators or to the frequency of the external modulation. A saturating nonlinearity turns this instability into a Hopf bifurcation generating a limit-cycle with the concomitant breaking of temporal symmetry and ergodicity . Furthermore, we show that even when the modulation amplitude is very small, monotonic growth can still occur but via a mechanism entirely different from that identified previously.

In trying to explain these collective behaviors, one is impressed by the similarities with a simple model introduced earlier  of two coupled parametric oscillators modulated periodically with a fixed phase difference $`\theta `$ between their modulations. The behavior of this dimer when $`\theta =\pi `$ is remarkably similar to that of the globally coupled infinite system, and the roots of the instabilities observed in the latter are already present in this simple system. At an even more primary level the seeds of these collective instabilities can be traced back to behaviors of single oscillators. In general an individual oscillator tends to synchronize with the external modulation whereas the coupling induces mutual synchronization between oscillators. These two synchronization tendencies cannot always be satisfied simultaneously. Coupling between oscillators can then be seen as leading to a sort of “selection” among the single oscillator modes, enhancing some (destabilization) and smoothing out others (stabilization).

In Section II we introduce the globally coupled model. Section III is a review of the properties of a single linear parametric oscillator. Section IV establishes the mathematical setting for the analysis of the globally coupled system as a mean field problem, and typical numerical results are presented in Section V. In Section VI we collect these results in the form of a phase diagram that characterizes the behavior of the coupled system as the modulation parameters are varied. Here we discuss the boundaries between stable and unstable behavior and also between different instability regimes. Section VII offers a comparison between the mean field system and the dimer. Conclusions and a summary are presented in Section VIII.

## II The basic linear model

Consider a set of $`N`$ unit mass linear parametric oscillators with displacements $`\{x_i\}`$, each with a periodically modulated frequency and all harmonically coupled to one another. Here we restrict our analysis to coupled linear parametric oscillators. The nonlinear case will be presented elsewhere . The equation of motion of the $`i`$-th oscillator is given by

$$\ddot{x}_i+\gamma \dot{x}_i+\omega _0^2[1+\xi _i(t)]x_i=\frac{k}{N}\underset{j=1}{\overset{N}{}}(x_ix_j),$$
(1)
with $`i=1,\mathrm{},N`$.

Analytic results are possible with a simple piecewise-constant periodic modulation (period=$`T_p`$), defined as

$$\xi _i(t)=A\text{sgn}[\mathrm{sin}\omega _p(t+\tau _i)]$$
(2)
where $`\omega _p=2\pi /T_p`$, and the initial phase $`\tau _i`$ is chosen at random for each oscillator. We are mainly interested in the mean amplitude

$$x=\frac{1}{N}\underset{i=1}{\overset{N}{}}x_i$$
(3)
as a measure of the macroscopic behavior of the system. In the thermodynamic limit $`N\mathrm{}`$, this site average (3) is equivalent to the average with respect to the random phase $`\tau _i`$ of the displacement of a single oscillator $`i`$,

$$x=\frac{1}{T_p}_0^{T_p}x_i𝑑\tau _i,$$
(4)
which is independent of $`i`$. Equation (1) can then be reduced to a single mean-field differential equation,

$$\ddot{x}+\gamma \dot{x}+\omega _0^2[1+\xi _\tau (t)]x=k(xx)$$
(5)
where we have dropped the index $`i`$ and have now indexed the modulation with the initial phase chosen at random for each oscillator. Note that the average $`x`$ must be evaluated self-consistently using Eqs. (4) and (5).

## III Review of The Single Parametric Oscillator

The further analysis of the coupled system is clarified by first briefly examining the properties of a single parametric oscillator subject to the piecewise-constant frequency modulation (2, for which the equation of motion is simply

$$\ddot{x}+\gamma \dot{x}+\omega _0[1+\xi _\tau (t)]x=0.$$
(6)
The arbitrary phase $`\tau `$ can be set to zero, but is retained explicitly in order to make clearer the connection with the coupled system. The equation of motion can be solved using Floquet theory (see the Appendix), and the result may be written as

$$\left(\begin{array}{c}x(t)\\ \dot{x}(t)+\frac{\gamma }{2}x(t)\end{array}\right)=𝖦_\tau (t)\left(\begin{array}{c}x(0)\\ \dot{x}(0)+\frac{\gamma }{2}x(0)\end{array}\right)$$
(7)
where $`𝖦_\tau (t)`$ is the time-evolution operator, whose explicit expression is given in the Appendix.

Floquet theory provides information about the state of the system at the end of each modulation period, but it does not determine the time evolution within each period. If, as in most of the literature, stability conditions were our only concern, finding the Floquet eigenvalues would be sufficient. However, we are interested not only in the conditions giving rise to instabilities but also in the types of instabilities. Since exponential decay or growth of the oscillator amplitude is expected in any case, the time evolution may be conveniently investigated by Laplace transformation,

$`\stackrel{~}{x}(s)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}e^{st}x(t)𝑑t`$ (8)
$`=`$ $`[\stackrel{~}{𝖦}_\tau (s)]_{11}x(0)+[\stackrel{~}{𝖦}_\tau (s)]_{12}\left({\displaystyle \frac{\gamma }{2}}x(0)+\dot{x}(0)\right),`$ (9)
where $`[\stackrel{~}{𝖦}_\tau (s)]_{ij}`$ is a matrix element of the Laplace transform of the time-evolution operator. In general, if $`\stackrel{~}{x}(s)`$ has poles at $`s_j=\mathrm{\Lambda }_j+i\mathrm{\Omega }_j`$, the temporal behavior of the amplitude is expressed as

$$x(t)=\underset{j}{}c_je^{\mathrm{\Lambda }_jt}e^{i\mathrm{\Omega }_jt}$$
(10)
where the $`c_j`$’s are constants determined by the inverse Laplace transform and the initial conditions,

$$c_j=2\pi i\underset{ss_j}{lim}(ss_j)\left\{\left[\stackrel{~}{𝖦}_\tau (s)\right]_{11}x(0)+\left[\stackrel{~}{𝖦}_\tau (s)\right]_{12}\left(\frac{\gamma }{2}x(0)+\dot{x}(0)\right)\right\}.$$
(11)
The Laplace transform of the time-evolution operator is derived in the Appendix as

$$\stackrel{~}{𝖦}_\tau (s)=\frac{\frac{1}{2}e^{\gamma T_p/2}\left[e^{(s+\gamma )T_p}I𝖦_\tau ^1(T_p)\right]}{\mathrm{cosh}\left[\left(s+\frac{\gamma }{2}\right)T_p\right]R}_0^{T_p}e^{st}𝖦_\tau (t)𝑑t,$$
(12)
where $`I`$ is the identity matrix and the resonance parameter $`R`$ is defined as

$$R=\mathrm{cos}\left(\frac{\omega _+T_p}{2}\right)\mathrm{cos}\left(\frac{\omega _{}T_p}{2}\right)\frac{\omega _+^2+\omega _{}^2}{2\omega _+\omega _{}}\mathrm{sin}\left(\frac{\omega _+T_p}{2}\right)\mathrm{sin}\left(\frac{\omega _{}T_p}{2}\right)$$
(13)
with

$$\omega _\pm =\omega _0\sqrt{1\pm A\left(\frac{\gamma }{2\omega _0}\right)^2}.$$
(14)
The poles of Eq. (12) are determined by the condition

$$\mathrm{cosh}\left[\left(s+\frac{\gamma }{2}\right)T_p\right]R=0,$$
(15)
and their real parts are explicitly given by

$$\mathrm{\Lambda }_j=\{\begin{array}{cc}\frac{\gamma }{2}\pm \frac{1}{T_p}\mathrm{ln}\left[R+\sqrt{R^21}\right]\hfill & \text{for}R1\hfill \\ \frac{\gamma }{2}\hfill & \text{for}|R|1\hfill \\ \frac{\gamma }{2}\pm \frac{1}{T_p}\mathrm{ln}\left[R+\sqrt{R^21}\right]\hfill & \text{for}R1,\hfill \end{array}$$
(16)
while the imaginary parts are

$`{\displaystyle \frac{\mathrm{\Omega }_j}{\omega _p}}=\{\begin{array}{cc}j\frac{1}{2}\hfill & \text{for}R1\hfill \\ j\hfill & \text{for}R1,\hfill \end{array}`$ (19)
$`{\displaystyle \frac{\mathrm{\Omega }_j^\pm }{\omega _p}}=j\pm {\displaystyle \frac{1}{2\pi }}\mathrm{arccos}(R)\text{for}|R|1.`$ (20)
Here $`\mathrm{arccos}(R)`$ lies in the range $`[0,\pi ]`$, $`j`$ is an integer. Note that the real parts of the poles are in fact independent of $`j`$, and therefore one can drop the subscript $`j`$ on $`\mathrm{\Lambda }`$ and rewrite the amplitude (10) as

$$x(t)=e^{\mathrm{\Lambda }t}\underset{j}{}c_je^{i\mathrm{\Omega }_jt}.$$
(21)
Correspondingly, the relative weights of the different modes remain the same for all times: they all die ($`\mathrm{\Lambda }<0`$), or they all explode ($`\mathrm{\Lambda }>0`$), or they all keep their initial amplitudes ($`\mathrm{\Lambda }=0`$). None becomes relatively more dominant with time.

Figure 1 illustrates the dependence of the poles on the resonance parameter. Consider first the real part $`\mathrm{\Lambda }`$. When $`|R|<1`$ it is single-valued, $`R`$-independent, and negative (unless there is no damping $`\gamma `$, in which case it vanishes and the motion is purely oscillatory). At the bifurcation points $`R=\pm 1`$ the real parts become $`R`$-dependent but remain negative until $`|R|`$ reaches the value

$$R_c\mathrm{cosh}\left(\frac{\gamma T_p}{2}\right)$$
(22)
which is greater than unity unless $`\gamma =0`$. Beyond $`|R|=R_c`$ one of the $`\mathrm{\Lambda }`$’s becomes positive, leading to exponential growth of the oscillator amplitude. The condition $`|R|=R_c`$ thus corresponds to the onset of parametric resonance or, more appropriately, of parametric instability.

Consider next the imaginary parts $`\mathrm{\Omega }_j`$, shown in Fig. 1 for $`j=2,\mathrm{}2`$ alternately as solid and dashed curves. The condition $`|R|=1`$ here also marks an interesting boundary at which the behavior of the imaginary parts changes from being $`R`$-independent when $`|R|>1`$ to $`R`$-dependent when $`|R|<1`$. On the large-$`|R|`$ side of this boundary the frequencies $`\mathrm{\Omega }_j`$ are simply proportional to the frequency $`\omega _p`$ of the modulation. When $`|R|<1`$, however, the oscillator frequencies change continuously with $`R`$ and bear no simple relation to either $`\omega _p`$ or the natural frequency $`\omega _0`$. In Fig. 2 we present the regions in parameter space $`(A,T_p/T_0)`$ where $`|R|>1`$. The darker regions correspond to $`R>1`$ and the lighter to $`R<1`$. The stability boundaries $`|R|=R_c`$ are also indicated by solid lines, the oscillator being unstable inside these boundaries. Note that the $`|R|=R_c`$ and $`|R|=1`$ frontiers almost coincide because of the very low damping ($`\gamma =0.1`$). As is well known, for small $`A`$ the instability appears in the vicinity of $`T_pjT_0/2`$ ($`j`$ integer).

Although the boundary $`|R|=R_c`$ is important in determining the transition from stable to unstable behavior for the single parametric oscillator, it does not play the same role for the coupled system, for which the boundary $`|R|=1`$ turns out to acquire further significance. Indeed, as we will see in the following section, each individual pole $`s_j=\mathrm{\Lambda }+i\mathrm{\Omega }_j`$ gives rise to a collective pole in the coupled system, and these collective poles have different $`\mathrm{\Lambda }`$’s. Therefore, some modes may become unstable ($`\mathrm{\Lambda }>0`$), and therefore dominant, even when $`|R|<1`$, while others remain stable ($`\mathrm{\Lambda }<0`$). A very striking example is that of the single oscillator mode with a zero imaginary part, i. e., the one with $`j=0`$ when $`R>1`$. If $`\mathrm{\Lambda }`$ is negative, it simply provides a monotonically decaying contribution to the oscillator displacement, and if $`\mathrm{\Lambda }`$ is positive it provides a monotonically growing contribution. However, its relative weight in the sum (21) is extremely small and therefore this mode is practically not detectable in simulations of a single oscillator. Its contribution is overwhelmed by those of the oscillatory modes. However, in the coupled system this non-oscillatory mode may become dominant because it may become unstable while the oscillatory modes are still stable. This will lead to a monotonic explosion of the mean.

## IV Collective Instabilities

The mean field equation (5) can be rearranged as

$$\ddot{x}+\gamma \dot{x}+\{\omega _0^2[1+\xi _\tau (t)]+k\}x=kx$$
(23)
which describes a single parametric oscillator of frequency

$$\omega _k=\omega _0\sqrt{1+\frac{k}{\omega _0^2}}$$
(24)
driven by an “additive force” $`kx`$. Therefore, it can be solved with the time-evolution operator (A13) of a single oscillator except that Eq. (14) must be replaced with the new shifted frequencies

$$\omega _\pm =\omega _0\sqrt{1\pm A+\frac{k}{\omega _0^2}\left(\frac{\gamma }{2\omega _0}\right)^2}.$$
(25)
The general solution of Eq. (23) can be written as

$$\left(\begin{array}{c}x(t)\\ \dot{x}(t)+\frac{\gamma }{2}x(t)\end{array}\right)=𝖦_\tau (t)\left(\begin{array}{c}x(0)\\ \dot{x}(0)+\frac{\gamma }{2}x(0)\end{array}\right)+k𝖦_\tau (t)_0^t𝖦_\tau (t^{})^1\left(\begin{array}{c}0\\ x(t^{})\end{array}\right)𝑑t^{}.$$
(26)
Taking the average of Eq. (26) with respect to the random phase, we obtain a self-consistent equation for the mean amplitude,

$$x(t)=[𝖦_\tau (t)]_{11}x(0)+[𝖦_\tau (t)]_{12}\left(\frac{\gamma }{2}x(0)+\dot{x}(0)\right)+k_0^tK(tt^{})x(t^{})𝑑t^{},$$
(27)
where the kernel is defined by

$$K(tt^{})=\left[𝖦_\tau (t)𝖦_\tau ^1(t^{})\right]_{12}.$$
(28)
The kernel (28) depends only on the time difference because of the uniform distribution of the initial phases $`\tau `$. One can therefore solve the integral equation (27) by Laplace transformation. The solution in Laplace space is given by

$$\stackrel{~}{x}(s)=\frac{[\stackrel{~}{G}_\tau (s)]_{11}x(0)+[\stackrel{~}{G}_\tau (s)]_{12}\left(\frac{\gamma }{2}x(0)+\dot{x}(0)\right)}{1k[\stackrel{~}{G}_\tau (s)]_{12}}.$$
(29)
Equation (29) has a set of poles determined by the condition

$$[\stackrel{~}{𝖦}_\tau (s)]_{12}=\frac{1}{k},$$
(30)
which differ from the ones derived for a single parametric oscillator. When $`k0`$ the right hand side of Eq. (30) diverges and the poles of the coupled system clearly approach those of the single oscillator. By explicitly performing the time integral in Eq. (12) and taking the average over the random phase $`\tau `$, one can obtain a rather complicated but explicit expression of Eq. (30) as

$`T_p`$ $`\omega _+\omega _{}\left[\mathrm{cosh}\left({\displaystyle \frac{\stackrel{~}{\gamma }T_p}{2}}\right)R\right]`$ (36)
$`\times \left(\omega _+^2+{\displaystyle \frac{\stackrel{~}{\gamma }^2}{4}}\right)^2\left(\omega _{}^2+{\displaystyle \frac{\stackrel{~}{\gamma }^2}{4}}\right)^2\left[{\displaystyle \frac{k\left(\omega _+^2+\omega _{}^2+\frac{\stackrel{~}{\gamma }^2}{2}\right)}{2\left(\omega _+^2+\frac{\stackrel{~}{\gamma }^2}{4}\right)\left(\omega _{}^2+\frac{\stackrel{~}{\gamma }^2}{4}\right)}}1\right]`$
$`k\omega _+(\omega _+^2\omega _{}^2)^2\left(\omega _{}^2{\displaystyle \frac{\stackrel{~}{\gamma }^2}{4}}\right)\mathrm{sin}\left({\displaystyle \frac{\omega _{}T_p}{2}}\right)\left[\mathrm{cosh}\left({\displaystyle \frac{\stackrel{~}{\gamma }T_p}{4}}\right)\mathrm{cos}\left({\displaystyle \frac{\omega _+T_p}{2}}\right)\right]`$
$`k\omega _{}(\omega _+^2\omega _{}^2)^2\left(\omega _+^2{\displaystyle \frac{\stackrel{~}{\gamma }^2}{4}}\right)\mathrm{sin}\left({\displaystyle \frac{\omega _+T_p}{2}}\right)\left[\mathrm{cosh}\left({\displaystyle \frac{\stackrel{~}{\gamma }T_p}{4}}\right)\mathrm{cos}\left({\displaystyle \frac{\omega _{}T_p}{2}}\right)\right]`$
$`k\omega _+\omega _{}\stackrel{~}{\gamma }(\omega _+^2\omega _{}^2)^2\mathrm{sinh}\left({\displaystyle \frac{\stackrel{~}{\gamma }T_p}{4}}\right)\left[2\mathrm{cosh}\left({\displaystyle \frac{\stackrel{~}{\gamma }T_p}{4}}\right)\mathrm{cos}\left({\displaystyle \frac{\omega _+T_p}{2}}\right)\mathrm{cos}\left({\displaystyle \frac{\omega _{}T_p}{2}}\right)\right]`$
$`=\mathrm{\hspace{0.33em}0}`$
where

$$\stackrel{~}{\gamma }\gamma +2s=\gamma +2\mathrm{\Lambda }+2i\mathrm{\Omega }.$$
(37)
Finding the complex roots of Eq. (36) is difficult even numerically. Therefore, we first investigate the collective modes graphically. Graphical inspection not only provides qualitative understanding of collective modes but also helps to identify suitable numerical algorithms.

Coupling between the oscillators causes a frequency shift of each oscillator as indicated in Eq. (25). This is a trivial effect that is accounted for by calculating all single oscillator quantities such as $`R`$ and the poles $`\mathrm{\Lambda }+i\mathrm{\Omega }`$ using the shifted frequency. The interesting questions concern the way in which the couplings between the shifted oscillators affect the dynamics, and, in particular, whether instabilities and synchronization may be caused or suppressed by the coupling. A great deal can be learned about these dynamics starting from the single oscillator roots shown in Fig. 1. We thus investigate the coupled system in three typical regimes: (1) A single oscillator with the given (shifted) parameters is unstable ($`RR_c`$); (2) The single oscillator is stable and in the regime where its frequency is determined by the modulation frequency ($`1RR_c`$); (3) The single oscillator is stable and its frequency bears no simple relation to the natural (shifted) oscillator frequency nor to the modulation frequency ($`|R|1`$). The results are presented in the next section.

## V Results

To find the poles associated with collective modes we need to solve the coupled equations

$`k\text{Re}[\stackrel{~}{𝖦}_\tau (s)]_{12}1=0`$ (38)
$`\text{Im}[\stackrel{~}{𝖦}_\tau (s)]_{12}=0`$ (39)
for $`\mathrm{\Lambda }`$ and $`\mathrm{\Omega }`$ in $`s=\mathrm{\Lambda }+i\mathrm{\Omega }`$. In the next set of figures we present results for the three cases listed at the end of the last section. Each case is presented in two or three figures.

In one set of figures we plot the left-hand sides of Eqs. (38) and (39) as contour lines in the space ($`\mathrm{\Lambda }/\omega _p`$, $`\mathrm{\Omega }/\omega _p`$) \[the solid lines for Eq. (38), and the broken ones for Eq. (39)\]. The parameters $`k=1.28`$, $`\omega _0=0.4`$, $`\gamma =0.01`$, and $`A=1.0`$ are used in all the figures. The thick lines indicate zero contours. Therefore, the solutions we seek are the intersections of these two sets of thick lines, indicated by open circles. The poles of single oscillators with shifted frequency are indicated by black solid circles. There is of course an infinite number of solutions but we only exhibit those that fall within the scale of our figures. Associated with each case is also a figure representing a number of trajectories as a function of time that characterize the particular situation. For some of the examples we also show the associated phase space trajectories. Note that in this presentation we make a careful distinction between a single oscillator and an individual oscillator. The former refers to an independent oscillator with parameters $`\omega _k`$ and $`\gamma `$, while the latter refers to one of the oscillators in a coupled system with parameters $`\omega _k`$, $`\gamma `$, and $`k`$.

Figure 3 shows the contours and poles for case (1). Here $`\omega _k=\omega _p/2`$, which corresponds to $`R=1.0062`$. Since $`R<1`$, the poles of single oscillators appear on the lines $`\mathrm{\Omega }_j/\omega _p=(j\frac{1}{2})`$ each as a pair because in this regime there are two values of $`\mathrm{\Lambda }`$ associated with each $`\mathrm{\Omega }`$. Only the poles for $`j=1`$ are shown because the others (and there are an infinite number of them) are off the scale of the figure. Since for these parameters $`|R|>R_c=1.00077`$, one of each pair of poles has a positive $`\mathrm{\Lambda }`$. The amplitude of a single oscillator would therefore diverge exponentially. However, as the numerical calculations indicate, all the collective modes have $`\mathrm{\Lambda }=\gamma /2`$ and thus the mean amplitude decays to zero despite the instability of single oscillators.

Computer simulation results for various trajectories of a system of 100,000 oscillators with these parameters are shown in Fig. 4. We see that the mean amplitude is indeed zero, and that the deviation $`\mathrm{\Delta }x=\sqrt{x^2x^2}`$ diverges. The inset shows the same two trajectories as well as the diverging trajectory of an individual oscillator. Phase space trajectories of individual oscillators in the coupled system are shown in Fig. 5. Each circle indicates a snapshot of an individual oscillator. Solid circles represent oscillators with positive modulation and open circles are ones with negative modulation at the time of the snapshot. Only 1000 oscillators out of 100,000 are shown.

The six snapshots show that with increasing time the phase volume increases (note the different scales in each snapshot), which is consistent with the divergent behavior of each oscillator and with the growth of the deviation $`\mathrm{\Delta }x`$, and provide confirmation that there is indeed no mutual synchronization nor other kinds of organized collective motion. However, the persistent separation of solid and open circles into separate quadrants indicates that individual oscillators are synchronized with the external modulation. Note that any individual oscillator moves clockwise, switching colors accordingly. However, since the phase of the modulation is chosen at random for each oscillator, at any time half of the oscillators have positive modulation and the other half negative, so that the pattern in preserved and an average over the phase gives $`x=0`$. In this case, the effect of the coupling is only the frequency shift (25). When the coupling strength is increased at fixed $`T_p/T_k`$, mutual synchronization will eventually be established. Then, the individual oscillators must lose synchronization with the modulation and consequently parametric instability is suppressed. Therefore, the coupling stabilizes the system in this parameter region.

Interesting behavior is observed when $`\omega _k=\omega _p`$, that is, the single oscillator frequency is equal to the modulation frequency. This corresponds to case (2) on our list. The contours and poles for this case are shown in Fig. 6. Since $`R=1.0001`$, the single oscillator has poles at $`\mathrm{\Omega }_j=j\omega _p`$. Only the poles for $`j=0`$ are shown. This condition is close to parametric resonance, but since $`R`$ is just below $`R_c=1.0031`$ all the single oscillator modes have negative $`\mathrm{\Lambda }`$’s and therefore single oscillator trajectories decay to the absorbing state $`x=0`$. However, one of the collective modes has a pole with $`\mathrm{\Omega }=0`$ and a positive $`\mathrm{\Lambda }`$. All the other collective modes are dominated by this $`j=0`$ unstable mode and therefore the mean amplitude diverges monotonically.

Computer simulation results for the associated trajectories are shown in Figure 7. The mean decays slowly at the beginning and then diverges monotonically. The deviation diverges as well, and does so more rapidly. The individual oscillator trajectory also diverges; in this example, although each single oscillator would be stable, the coupling causes individual oscillators in the system to become unstable. In other words, each oscillator is driven by the diverging mean on the right hand side of Eq. (23). The phase space trajectories in Fig. 8 show that after an initial transient (first three panels, where open circles are hiding most of the solid circles), individual oscillators in the coupled system oscillate with increasing amplitude about $`x`$, while the mean $`x`$ is moving away from the origin. In the long time limit, each oscillator “forgets” its initial conditions and it is driven by the mean. Therefore, while the phase of each individual oscillator is determined by the phase of the modulation, the amplitude in phase space is determined by the mean. Correspondingly, the oscillators become “amplitude-synchronized” through the mean. Until the synchronization is well established, the individual oscillators decay because the single oscillator modes have negative $`\mathrm{\Lambda }`$.

As in the previous case, the phases of all single-oscillator modes with $`\mathrm{\Omega }_j0`$ become synchronized with the external modulation. The separation of solid and open circles in the three later-time panels in Fig. 8 indicates this synchronization. However, in contrast with the previous case, there is now a zero-frequency mode which does not have a phase to be synchronized. This mode is therefore not affected by the phase of each oscillator nor by the phase of the modulation. The zero-frequency mode shifts the center of oscillation away from $`x=0`$ in either the positive or the negative direction. In the presence of coupling the oscillators tend to follow the mean, and therefore shift in the same direction, thus breaking the system symmetry. Note that the open and solid circles no longer lie entirely in separate quadrants. It is precisely the excess of open circles relative to solid ones in the positive quadrants (and vice versa in the negative quadrants) that drives the motion of the mean further to the right. This mode can now dominate the behavior of the system because, contrary to the situation of a single oscillator, each mode has a different $`\mathrm{\Lambda }`$. In this example the effect of coupling is thus not only to shift the frequency according to (25), but also the more interesting collective symmetry-breaking monotonic divergence of the mean amplitude and the mutual synchronization of the individual oscillator amplitudes.

Since the role of the coupling in this case is to suppress the single-oscillator modes with $`\mathrm{\Omega }_j0`$ and to allow the mode with $`\mathrm{\Omega }_j=0`$ to grow, this instability persists even for very large $`k`$. It should be noted that this collective instability is a different phenomenon from the monotonic growth reported previously. That instability occurs only for large $`A`$ and disappears when the coupling strength $`k`$ is increased (re-entrant transition).

Finally, case (3) occurs when $`\omega _k`$ is between the two previous cases, so that $`|R|<1`$ and single oscillators are not in parametric resonance condition. However, Fig. 9 indicates that at least one of the collective modes has a positive $`\mathrm{\Lambda }`$ as well as non-zero $`\mathrm{\Omega }`$. Therefore, the mean $`x`$ oscillates with diverging amplitude. Recall that when $`|R|<1`$ the eigenfrequencies of single oscillators vary continuously with $`R`$ but do not match the frequency of the modulation nor the natural frequency of the oscillator. Therefore, the individual synchronization to the external modulation plays no role and the phases of the oscillators are free to synchronize to one another through a synchronization to the phase of the mean. Computer simulations shown in Fig. 10 confirm the oscillatory instability, and the phase space points of individual oscillators shown in Fig. 11 testify to this phase synchronization. Although solid and open circles again form separate groups, the entire ring of open and solid circles alternates between the positive quadrants and the negative quadrants. Therefore, all oscillators are mutually synchronized.

In the next section we collect the results for the coupled system into a phase diagram indicating regions of stability and instability of different types.

## VI Phase Diagram

A convenient way to summarize observations and characterize the instabilities systematically is in the form of appropriate “phase diagrams” in which the stability boundaries are presented as a function of the system parameters. Since there are many parameters in this model the full diagram would involve a many-dimensional representation. Instead, we present the diagram in the two-dimensional space $`(A,T_p/T_k)`$ that characterizes the external modulation for a given set of oscillator parameters $`\omega _k`$, $`k`$, and $`\gamma `$.

Figure 12 shows the phase diagram for the system parameters indicated in the caption. The colored regions denote unstable regimes, each color coding a particular type of instability. The yellow region denotes parameter ranges where the individual oscillators are unstable but the mean amplitude for the coupled system is zero (“incoherent unstable oscillations”). This region is associated with divergences of the first term in Eq. (26) but not of the second term. In this case the distinction between “single oscillators” and “individual oscillators” becomes moot, since the right hand side of Eq. (23) plays no role. Oscillatory instabilities of the mean with positive $`\mathrm{\Lambda }`$ and nonzero $`\mathrm{\Omega }`$ are denoted in pink (“unstable spirals”). Blue regions (“ saddle nodes”) denote monotonic divergence of the mean with one positive $`\mathrm{\Lambda }`$ and zero $`\mathrm{\Omega }`$. Green regions (“unstable nodes”) also denote monotonic instabilities but with two positive $`\mathrm{\Lambda }`$’s and zero $`\mathrm{\Omega }`$. These latter three instabilities involve divergences of the second term in Eq. (26). The first term may diverge as well, i.e., the $`\mathrm{\Lambda }`$’s of the single oscillator and of the collective mode may be simultaneously positive. When this is the case we always observe the $`\mathrm{\Lambda }`$ of the collective mode to be larger than that of the single oscillator modes; the collective mode therefore dominates the dynamics. In establishing the terminology for various instabilities we have loosely followed the usual conventions of nonlinear dynamics.

We can associate each specific case discussed in the previous section with a location on this sort of phase diagram. Thus the trajectories in Fig. 4 are in a yellow regime of incoherent oscillatory divergence of each oscillator. Those of Fig. 7 are in a blue or green regime of monotonic divergence of the mean, and those of Fig. 10 are in the pink regime of unstable spirals.

It is helpful to follow the behavior of the oscillator system across the various collective instability boundaries by considering the signs of $`\mathrm{\Omega }`$ and $`\mathrm{\Lambda }`$ as one increases the modulation amplitude $`A`$ (thus moving upward vertically along the phase diagram) for different fixed values of the modulation period $`T_p`$. Various associated bifurcation diagrams presenting $`\mathrm{\Lambda }`$ (solid lines) and $`\mathrm{\Omega }`$ (dotted lines) as a function of $`A`$ are shown in Fig. 13. Consider first the period $`T_p=0.75T_k`$, shown in panel (a). As $`A`$ increases, $`\mathrm{\Lambda }`$ changes sign, becoming positive at $`A_c=1.14`$ while $`\mathrm{\Omega }`$ remains positive. This represents a transition from a stable spiral to an unstable spiral (pink region in Fig. 12).

Consider next the period $`T_p=2.04T_k`$, shown in panel (b). Here $`\mathrm{\Lambda }`$ becomes positive at $`A_c=3.06`$ while $`\mathrm{\Omega }`$ remains positive. This therefore again marks a transition from a stable spiral to an unstable spiral (pink region). However, with a further increase in amplitude, $`\mathrm{\Omega }`$ eventually goes to zero at $`A_c^{}=3.51`$, where $`\mathrm{\Lambda }`$ bifurcates to two positive values (an unstable node, green region) via a saddle-node bifurcation. The oscillatory instability thus switches to a monotonic instability at this point.

A different transition pattern is seen when $`T_p=3.0T_k`$, shown in panel (c). It begins with a stable spiral and switches to a stable node at $`A_c=2.91`$. With a further increase in amplitude the system undergoes a transition to a saddle node (blue region) at $`A_c^{}=3.04`$.

A more complex transition pattern is shown in panel (d), in which the character of the instability changes several times along the $`T_p=3.9T_k`$ line. As usual, at low amplitudes there is a stable spiral. At the point $`A_c=2.99`$ the system moves into an unstable spiral (pink region). The unstable spiral becomes an unstable node (very small green region in the phase diagram) via a saddle-node bifurcation at $`A_c^{}=3.30`$. A further transition to a saddle node (blue) occurs at $`A_c^{\prime \prime }=3.4`$.

We note that the phase diagram just described is quite rich and intricate. Instabilities cover even larger regions in parameter space, and do so with increasing intricacy, as the damping $`\gamma `$ decreases. A typical phase diagram for low damping is shown in Fig. 14.

Increasing coupling or damping also leads to greater stability: the stability boundaries move “up” in the phase diagram when $`\gamma `$ and/or $`k`$ are increased, indicating that a stronger modulation is needed to cause unstable behavior. Furthermore, the oscillatory instabilities eventually disappear with increasing modulation period $`T_p`$, leaving only the monotonic collective instabilities. However, it should be noted that the monotonic instabilities for sufficiently large modulation period are simply due to an inversion of the effective harmonic potential and hence not due to any collective effects. That there is an inversion can already be suspected from the fact that at least one of the shifted frequencies in Eq. (25) becomes imaginary.

An analysis of the system for large $`T_p`$ is fairly simple and instructive in elucidating the source of instabilities more explicitly. In the adiabatic limit $`T_p\mathrm{}`$ the single oscillator frequencies are frozen in time, half of them at frequency $`\stackrel{~}{\omega }_+`$ and the other half at frequency $`\stackrel{~}{\omega }_{}`$, where

$$\stackrel{~}{\omega }_\pm =\omega _0\sqrt{1+\frac{k}{2\omega _0^2}\pm A}.$$
(40)
The mean field equations of motion then are

$`\ddot{x}_++\gamma \dot{x}_++\stackrel{~}{\omega }_+^2x_+`$ $`=`$ $`{\displaystyle \frac{k}{2}}x_{}`$ (41)
$`\ddot{x}_{}+\gamma \dot{x}_{}+\stackrel{~}{\omega }_{}^2x_{}`$ $`=`$ $`{\displaystyle \frac{k}{2}}x_+`$ (42)
where $`\mathrm{}_\pm `$ indicates an average over the oscillators with the frequency $`\stackrel{~}{\omega }_\pm `$ respectively. This system can be diagonalized analytically in full generality. The eigenmodes of the coupled system are characterized by the complex frequencies

$`\mathrm{\Omega }_\pm ^{(1)}`$ $`=`$ $`i{\displaystyle \frac{\gamma }{2}}\pm \omega _0\sqrt{1+{\displaystyle \frac{k}{2\omega _0^2}}{\displaystyle \frac{\gamma ^2}{4\omega _0^2}}+\sqrt{A^2+\left({\displaystyle \frac{k}{2\omega _0^2}}\right)^2}}`$ (43)
$`\mathrm{\Omega }_\pm ^{(2)}`$ $`=`$ $`i{\displaystyle \frac{\gamma }{2}}\pm \omega _0\sqrt{1+{\displaystyle \frac{k}{2\omega _0^2}}{\displaystyle \frac{\gamma ^2}{4\omega _0^2}}\sqrt{A^2+\left({\displaystyle \frac{k}{2\omega _0^2}}\right)^2}}.`$ (44)
This case clearly illustrates the distinction between what we have called “single oscillator instabilities” and “chain instabilities”. The former refer to the frequency (40) while the latter refer to (44). While the single oscillators would remain stable until $`A=1+k/2\omega _0^2`$ (at which point $`\stackrel{~}{\omega }_{}`$ becomes imaginary), the chain becomes destabilized when $`A`$ reaches the value $`\sqrt{1+k/\omega _0^2}`$, where the imaginary part of $`\mathrm{\Omega }_{}^{(2)}`$ becomes negative (note that the transition point is independent of $`\gamma `$). Beyond that the system is in a saddle-point/unstable-node instability region of non-oscillatory exponential growth due to potential inversion, which was the primary reason of instability in the previous work. In the phase diagrams Figs. 12 and 14 this translates to a stability boundary that settles at $`A=3`$ as $`T_p\mathrm{}`$. Note that in particular the boundary $`A=\sqrt{1+k/\omega _0^2}`$ remains valid in the overdamped limit $`\gamma \mathrm{}`$, consistent with our early work.

The origin of the instabilities presented as narrow blue and pink tongues in the low-$`A`$ region of Fig. 14 is entirely different from the mechanism based on the temporarily inverted potential and unique to the underdamped case. In the following sections we will explain the cause of these instabilities using a dimer of parametric oscillators.

## VII Comparison With Parametric Oscillator Dimer

In a recent paper  we reported results for a model of two coupled oscillators subject to parametric modulations with a phase difference $`\theta `$. The equations of motion for this system are just the $`N=2`$ version of Eq. (1):

$`\ddot{x}_1`$ $`=`$ $`\omega _0^2[1+\xi _1(t)]x_1{\displaystyle \frac{k}{2}}(x_1x_2)\gamma \dot{x}_1`$ (45)
$`\ddot{x}_2`$ $`=`$ $`\omega _0^2[1+\xi _2(t)]x_2{\displaystyle \frac{k}{2}}(x_2x_1)\gamma \dot{x}_2.`$ (46)
The piecewise-constant periodic modulations of the two oscillators differ by a constant phase $`\theta `$, so that we can write Eq. (2) for this case as

$`\xi _1(t)`$ $`=`$ $`A\text{sgn}\left(\mathrm{sin}(\omega _pt)\right)`$ (47)
$`\xi _2(t)`$ $`=`$ $`A\text{sgn}\left(\mathrm{sin}(\omega _pt+\theta )\right).`$ (48)
We want to investigate whether the mean position $`x(x_1+x_2)/2`$ reproduces the macroscopic behavior of the mean in the globally coupled model.

In the absence of parametric modulations, a dimer has two eigenmodes: symmetric (or mutually synchronized, $`x_1(t)=x_2(t)`$) and antisymmetric (or mutually antisynchronized, $`x_1(t)=x_2(t)`$), with the former having a lower energy. When the parametric modulations are applied, these modes are in general no longer the eigenmodes of the dimer (except for $`\theta =0`$). However, the motions can be expressed as linear combinations of these modes and, in particular, the behavior of the mean of interest is reflected in the excitation of the symmetric mode by the parametric modulations.

In previous work  we focused on the behavior of the system as a function of the parameters $`A`$, $`T_p=2\pi /\omega _p`$, $`\omega _0`$, $`\gamma `$, $`k`$, and $`\theta `$. Among our conclusions is the fact that the regions of parametric instability are sensitively dependent on the phase difference $`\theta `$. Of particular interest for the analysis in this paper is the behavior of the anti-phased dimer ($`\theta =\pi `$). This particular dimer captures many of the features of the mean field coupled system with unexpected accuracy. This assertion, which was originally based on our previous comparison of the regions of parametric resonance , is reinforced when the dimer bifurcation diagrams are further refined to take into account details of qualitatively different trajectories, as we shall see below.

For the piecewise constant parametric modulation (47) the solution of the dimer problem is formally simple. The stability analysis is based on the eigenvalues $`\{\lambda _j\}`$ ($`|\lambda _1||\lambda _2||\lambda _3||\lambda _4|`$) of the Floquet operator : parametric instability occurs when $`|\lambda _1|>1`$. To characterize different types of parametric instability we present bifurcation diagrams using the same color conventions as in the phase diagrams of the mean field model (Figs. 12 and 14). If $`\text{Im }\lambda _10`$, then clearly one has an oscillatory instability (pink regions). If $`\text{Im }\lambda _1=0`$, the instability can be either oscillatory or monotonic. Using the eigenvector corresponding to $`\lambda _1`$ as initial condition, we have determined whether or not $`x`$ crosses zero during one period of the modulation. If it does, the point is assigned to a pink region. If it does not, the second largest eigenvalue $`\lambda _2`$ will determine whether the point belongs to a blue ($`|\lambda _2|1`$) or green ($`|\lambda _2|>1`$) region. The pink, blue and green regions are all caused by the instability of the symmetric mode. The yellow region, on the other hand, requires the instability of the antisymmetric mode and the decay of the symmetric mode. However, such a purely antisynchronous solution is forbidden in the $`\theta =\pi `$ case , which means that yellow regions cannot appear at all in the anti-phased dimer .

Results for relatively large damping ($`\gamma /\omega _0=0.4`$) are presented in the bifurcation diagram of Fig. 15, which should be compared to Fig. 12. The similarity between the two figures is remarkable. Despite some extra green regions and the absence of the yellow tongue, one notices that the principal resonance regions of the dimer bifurcation diagram ($`1/2T_p/T_k3/2`$) fits the same region in the mean field model almost exactly. The green regions connected to pink regions in the mean field model are also well mimicked by the dimer.

Figure 16 illustrates the bifurcation diagram for a small value of damping ($`\gamma /\omega _0=0.01`$). Comparison with Fig. 14 shows that although the agreement between the two models is not as good as for higher values of $`\gamma /\omega _0`$, the basic structure and similarities of the phase diagram and the bifurcation diagram are nonetheless preserved. The main differences are the complex pink patterns in region $`T_p/T_k3`$ of the dimer. Also, as in the high-gamma case, the dimer has larger green regions than the mean field model, suggesting that the coupling in the latter plays a stronger role in stabilizing the system. But in spite of these differences and most importantly, the principal resonance region ($`1/2T_p/T_k3/2`$) again shows an almost perfect match.

In the dimer, competition between two kinds of synchronization plays a key role in the destabilization of the system: on the one hand, synchronization between each oscillator and its modulation; on the other hand, synchronization between the two oscillators. This competition is essentially governed by the values of $`A`$ and $`k`$. Larger values of $`A`$ favor the former, while larger values of $`k`$ favor the latter. When the coupling is weak, the energy difference between symmetric and antisymmetric modes is small and both can be excited. In this case, the individual oscillators are nearly independent and the stability diagram of the dimer is similar to that of a single oscillator. As the coupling strength increases, the energy of antisymmetric oscillations increases until eventually only in-phase oscillations are energetically accessible. This mutually synchronized motion brings the system out of synchronization with the modulation.

The stability diagram of the anti-phased dimer in the $`(T_p/T_k,k)`$ plane shown in Fig. 17 illustrates the above story. First, consider $`T_p/T_k0.5`$, where a single oscillator is in the main parametric instability region ($`R<1`$, the first light grey region in the lower panel of Fig. 17). For small $`k`$, the dimer is also unstable and still dominated by the antisymmetric mode (even though the symmetric mode cannot disappear, as mentioned above). However, as $`k`$ increases, excitation of the antisymmetric mode becomes more difficult and the symmetric mode becomes dominant, which stabilizes the system. When $`0.5T_p/T_k1`$ a single oscillator is stable ($`|R|<1`$, between the dark and light grey regions). In this parameter region, individual oscillators do not have to be synchronous to the modulation (see Fig. 1). They are free to become mutually synchronized and the system becomes unstable above a certain coupling strength. Since the symmetric mode dominates, this instability persists even for large $`k`$. Finally, when $`T_p/T_k1`$, the single oscillator is again in an unstable region ($`R>1`$, the dark grey region). Although the situation is similar to the first case, the individual oscillators now have a zero-frequency mode which is not subject to synchronization with the modulation. Therefore, the zero-frequency mode of the two oscillators can be mutually synchronized, which produces a monotonic growth of the mean. Again, this instability does not have strong dependence on the coupling strength and persists even for large $`k`$.

The good agreement between the dimer and mean field models is not merely a coincidence. Consider one particular test oscillator in the globally coupled system. Now define the group of its “friends” as the set of all oscillators whose modulation phase lies within an interval of $`\pm \pi /2`$ around its phase. The group of its “enemies” comprises the set of all other oscillators, whose average modulation phase is opposite to that of the test oscillator. When the oscillators are synchronized to the modulations, “friends” are also mutually synchronized to one another, regardless of the coupling. When the coupling increases, there is a competition between two kinds of synchronization: synchronization between the “enemies” and “friends”, and synchronization with their own modulations. This situation is similar to the anti-phased dimer, which helps explain the remarkable similarity between the instability diagrams of the dimer and mean field models. Notice that this description becomes exact in the quenched limit, where the “friends” and “enemies” can be represented respectively by $`x_+`$ and $`x_{}`$. From this perspective, the choice $`\theta =\pi `$ for the dimer appears as a natural one, being not only a particularly symmetric case in the general dimer problem  but also the effective phase difference between the two groups.

## VIII Discussion and Conclusions

We have investigated collective instabilities of an infinite set of globally coupled linear oscillators driven by time-periodic piecewise linear modulations with random initial phases. These instabilities occur in certain parameter regimes (and not in others), and we have produced phase diagrams as a function of system parameters indicating detailed stability boundaries and types of instabilities.

Instabilities arise from phase synchronization (although not all phase synchronization leads to instability), and there are two possible competing synchronization mechanisms: synchronization of individual oscillators with the external modulation (“modulation synchronization”), and mutual synchronization between oscillators (“mutual synchronization”). In the absence of the external modulation, only mutual synchronization is possible. On the other hand, when the coupling is absent, a single oscillator in parametric resonance condition ($`|R|>1`$) synchronizes with the external modulation. Note that parametric instability of a single oscillator requires modulation synchronization, and in general any instability requires an external modulation to pump energy into the system.

In the presence of external modulation with random phases and coupling, even in a parameter regime where either alone would lead to synchronization, it is not possible for both types of synchronization to occur simultaneously. Intuition might lead to the conclusion that sufficiently strong coupling would favor mutual synchronization thereby destroying modulation synchronization (and hence stabilizing the system). Intuition might also lead one to conclude that weak coupling necessarily results in a dominance of modulation synchronization and hence to instabilities of individual oscillators not related to one another. However, we have shown that this intuitive picture would be incomplete because it does not account for the existence of collective parametric instabilities due to the coupling and entirely absent in a single parametric oscillator. Here we discuss our results in terms of the two competing synchronization mechanisms.

In a globally coupled model, the individual oscillators are coupled to the mean $`x`$, and thus mutual synchronization can be thought of as synchronization between the individuals and the mean. Note that there is no mutual synchronization when $`x=0`$ (because if there were, the mean would not be zero).

Now suppose we are in the parameter regime $`R<1`$ where individual uncoupled oscillators are parametrically unstable. When these oscillators are coupled, one can imagine one of two possible scenarios. If the coupling leads to mutual synchronization, the individual oscillators can no longer be synchronous with the external modulation and therefore the coupled system has been stabilized by the coupling. On the other hand, if the coupling does not lead to mutual synchronization but instead there is modulation synchronization, then the oscillators may be individually unstable but with $`x=0`$. Our results show that the second scenario is the correct one for sufficiently small values of $`k`$, as shown in the yellow regions of incoherent instability in Fig. 14. Modulation synchronization has “won”. On the other hand, there is a coupling energy cost to the lack of mutual synchronization, which slows down the instability of individual oscillators relative to their uncoupled amplitude growth. For larger $`k`$ the first scenario takes over, and the yellow region disappears above a certain value of the coupling.

Next suppose that we are in the other parametric resonance regime, $`R>1`$. The situation is in some ways similar to the previous case but there is a major difference: there is now a mode, the $`j=0`$ mode, whose frequency is zero and which therefore (contrary to the other modes) need not (indeed cannot) synchronize to the modulation. The amplitude of this mode can grow monotonically in either direction, and the coupling among oscillators leads to a tendency for the zero-frequency mode of all the individual oscillators to move in the same direction. Thus while the growth rate of the $`j0`$ modes is reduced by coupling due to the lack of mutual synchronization, that of the $`j=0`$ mode is enhanced because the coupling fosters mutual synchronization of this mode.

If the individual oscillators are not in regimes of parametric instability ($`|R|<1`$), there is no synchronization to the modulation and the oscillators are free to synchronize with one another. Mutual synchronization is thus established and the mean $`x`$ becomes oscillatory with the same frequency as that of individual oscillators in the coupled system. The oscillatory mean drives the system into unstable states via the mean field coupling.

We have found that the instability bifurcation diagram for a simple anti-phased dimer model reproduces the phase diagram for the mean field system with surprising accuracy.

In the absence of external modulation, a dimer has two eigenmodes: symmetric or mutually synchronized (lower energy), and antisymmetric or mutually anti-synchronized (higher energy). In the presence of time-periodic piecewise linear modulations which are exactly out of phase on the two coupled oscillators ($`\theta =\pi `$), a competition between these two modes (which are no longer eigenmodes) ensues. This competition is in many ways similar to the competition between modulation synchronization and mutual synchronization described for the globally coupled system, and here again it determines the stability of a dimer. When the coupling is weak, the energy difference between symmetric and antisymmetric modes is small and both can be excited. In this case, the individual oscillators are nearly independent and the stability diagram of the dimer is similar to that of a single oscillator. The synchronization of each oscillator with its modulation dominates the behavior, and instabilities thus represent boundless excitation of the antisymmetric mode. With increasing coupling the energy of the antisymmetric mode increases until it is too high to be excited. Only the symmetric mode can be excited, i.e., the oscillators become mutually synchronized. The synchronization with the modulation is thus destroyed and the associated parametric instability is suppressed.

Although the similarity between the dimer and globally coupled models is remarkable, they also exhibit various important differences. In the dimer model, mutual synchronization involves only two oscillators. On the other hand, in the global coupling model an oscillator must be synchronous with essentially all others to create collective motion. Therefore, in the thermodynamic limit $`N\mathrm{}`$, the collective instability in the globally coupled system is a genuine phase transition, whereas the instabilities in the dimer are simple bifurcations. Nevertheless, the stability boundaries and dynamics of the mean amplitudes in both cases show impressive similarities.

An interesting case that is in some sense “in between” these two and that promises interesting new features is that of a one-dimensional chain of oscillators with nearest-neighbor coupling. When the phase of the modulations of the oscillators in the chain is chosen at random, there is a significant chance that both neighbors of any given oscillator have a phase “similar” (suitably defined within some range) to its phase. In this case, the middle oscillator can easily establish simultaneously both mutual synchronization with its neighbors and synchronization with the modulation. Therefore, locally this oscillator could become unstable. On the other hand, if the neighbors of a given oscillator are modulated with phases opposite to its own modulation phase, that oscillator may be stabilized. Therefore, the spatial pattern of the modulation phase is expected to play an important role and instability may become wavelength dependent, suggesting spatial pattern formation. Such patterns can of course not be observed in either a dimer or a globally coupled model. A detailed analysis of such systems will be presented elsewhere .

## IX Acknowledgments

This work was supported in part by the National Science Foundation under grant No. PHY-9970699. Support of the Program on Inter-University Attraction Poles of the Belgian Government and the F. W. O. Vlaanderen (CVdB) is also gratefully acknowledged.

## A Time-evolution operator of a single parametric oscillator

We solve the equation of motion (6) using a standard Floquet method. The damping term can be eliminated by introducing a new variable $`y`$ defined as

$$x(t)=e^{\gamma t/2}y(t).$$
(A1)
Equation (6) becomes

$$\ddot{y}+\omega ^2(t)y=0$$
(A2)
with the time-dependent frequency

$$\omega ^2(t)\omega _0^2[1+\xi _\tau (t)]\gamma ^2/4.$$
(A3)
The solution of the undamped frequency-modulated oscillator (A2) can be expressed in terms of the time-evolution operator from $`t=0`$ to $`t`$, $`𝗀_\tau (t)`$, as

$$\left(\begin{array}{c}y(t)\\ \dot{y}(t)\end{array}\right)=𝗀_\tau (t)\left(\begin{array}{c}y(0)\\ \dot{y}(0)\end{array}\right)$$
(A4)
For a piecewise constant modulation such as Eq. (2), the explicit form of the time-evolution operator is known. Using its periodicity and composition property we note that if $`t=nT+u`$ then

$$𝗀_\tau (t)=𝗀_\tau (nT+u)=𝗀_\tau (u)𝗀_\tau (nT)=𝗀_\tau (u)\left[𝗀_\tau (T)\right]^n.$$
(A5)
It is thus sufficient to find $`𝗀_\tau (t)`$ for $`t[0,T]`$

When the phase is $`\tau [0,T/2]`$, the frequency varies as

$$\omega (t)=\{\begin{array}{cc}\omega _+\hfill & \text{for}t[0,\frac{T}{2}\tau )\hfill \\ \omega _{}\hfill & \text{for}t[\frac{T}{2}\tau ,T\tau )\hfill \\ \omega _+\hfill & \text{for}t[T\tau ,T]\hfill \end{array}$$
(A6)
while for $`\tau [T/2,T]`$,

$$\omega (t)=\{\begin{array}{cc}\omega _{}\hfill & \text{for}t[0,T\tau )\hfill \\ \omega _+\hfill & \text{for}t[T\tau ,\frac{3T}{2}\tau )\hfill \\ \omega _{}\hfill & \text{for}t[\frac{3T}{2}\tau ,T]\hfill \end{array}$$
(A7)
where the $`\omega _\pm `$ are defined in Eq. (14). During each constant frequency time window the system evolves according to the propagator of a simple harmonic oscillator of the appropriate frequency. This propagator is well known:

$$𝗀_\pm (t,t^{})=\left(\begin{array}{cc}\mathrm{cos}[\omega _\pm (tt^{})]& \frac{1}{\omega _\pm }\mathrm{sin}[\omega _\pm (tt^{})]\\ \omega _\pm \mathrm{sin}[\omega _\pm (tt^{})]& \mathrm{cos}[\omega _\pm (tt^{})]\end{array}\right).$$
(A8)
The full operator $`𝗀_\tau (t)`$ can be expressed as products of the $`g_\pm `$. For the case (A6),

$$𝗀_\tau (t)=\{\begin{array}{cc}𝗀_+(t,0)\hfill & \text{for}t[0,\frac{T}{2}\tau )\hfill \\ 𝗀_{}(t,\frac{T}{2}\tau )𝗀_+(\frac{T}{2}\tau ,0)\hfill & \text{for}t[\frac{T}{2}\tau ,T\tau )\hfill \\ 𝗀_+(t,T\tau )𝗀_{}(T\tau ,\frac{T}{2}\tau )𝗀_+(\frac{T}{2}\tau ,0)\hfill & \text{for}t[T\tau ,T]\hfill \end{array}$$
(A9)
and for the case (A7),

$$𝗀_\tau (t)=\{\begin{array}{cc}𝗀_{}(t,0)\hfill & \text{for}t[0,T\tau )\hfill \\ 𝗀_+(t,T\tau )𝗀_{}(T\tau ,0)\hfill & \text{for}t[T\tau ,\frac{3T}{2}\tau )\hfill \\ 𝗀_{}(t,\frac{3T}{2}\tau )𝗀_+(\frac{3T}{2}\tau ,T\tau )𝗀_{}(T\tau ,0)\hfill & \text{for}t[\frac{3T}{2}\tau ,T].\hfill \end{array}$$
(A10)
These expressions can be further simplifed by using time-translation symmetry, $`𝗀_\pm (t+u,t^{}+u)=𝗀_\pm (t,t^{})=𝗀_\pm (tt^{},0)`$. We can thus simplify our notation and set $`𝗀_\pm (t,0)𝗀_\pm (t)`$ with no loss of information.

In particular, for $`t=T`$

$$𝗀_\tau (T)=𝖲_\tau 𝗀_0(T)𝖲_\tau ^1$$
(A11)
where

$$\begin{array}{ccc}𝗀_0(T)𝗀_{}\left(\frac{T}{2}\right)𝗀_+\left(\frac{T}{2}\right),\hfill & 𝖲_\tau 𝗀_+(\tau ),\hfill & \text{for}\tau [0,\frac{T}{2}]\hfill \\ 𝗀_0(T)𝗀_+\left(\frac{T}{2}\right)𝗀_{}\left(\frac{T}{2}\right),\hfill & 𝖲_\tau 𝗀_{}\left(\tau \frac{T}{2}\right),\hfill & \text{for}\tau [\frac{T}{2},T].\hfill \end{array}$$
(A12)
Since only $`𝖲_\tau `$ depends on $`\tau `$, the trace, determinant, and eigenvalues of $`𝗀_\tau (T)`$ do not depend on the random phase.

Transforming back to the original variables, we finally obtain the time-evolution operator for $`x(t)`$ and $`\dot{x}(t)`$,

$$𝖦_\tau (t)=e^{\gamma t/2}𝗀_\tau (t).$$
(A13)
Next we derive an explicit expresion for the Laplace transform of the evolution operator with the help of Eq. (A5):

$`\stackrel{~}{𝖦}_\tau (s)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}e^{st}𝖦_\tau (t)𝑑t`$ (A14)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}\left[e^{sT}𝖦_\tau (T)\right]^m{\displaystyle _0^T}e^{st}𝖦_\tau (t)𝑑t`$ (A15)
$`=`$ $`\left(Ie^{sT}𝖦_\tau (T)\right)^1{\displaystyle _0^T}e^{st}𝖦_\tau (t)𝑑t`$ (A16)
$`=`$ $`{\displaystyle \frac{\frac{1}{2}e^{\gamma T/2}\left[e^{(s+\gamma )T}I𝖦_\tau (T)\right]_0^Te^{st}𝖦_\tau (t)𝑑t}{\mathrm{cosh}[\left(s+\frac{\gamma }{2}\right)T]R}}`$ (A17)
where $`I`$ is an identity matrix and

$`R`$ $`=`$ $`{\displaystyle \frac{1}{2}}e^{\gamma T/2}\text{Tr }𝖦_0(T)`$ (A18)
$`=`$ $`{\displaystyle \frac{1}{2}}\text{Tr }𝗀_0(T)`$ (A19)
$`=`$ $`\mathrm{cos}({\displaystyle \frac{\omega _+T}{2}})\mathrm{cos}({\displaystyle \frac{\omega _{}T}{2}}){\displaystyle \frac{\omega _+^2+\omega _{}^2}{2\omega _+\omega _{}}}\mathrm{sin}({\displaystyle \frac{\omega _+T}{2}})\mathrm{sin}({\displaystyle \frac{\omega _{}T}{2}})`$ (A20)
The geometric series in Eq. (A17) converges only when

$$e^{sT}𝖦_0(T)<1.$$
(A21)


format: text

text:
# Escape time in anomalous diffusive media

## I Introduction

The old problem of surmounting a potential barrier, known as Kramers’ problem, is doubtlessly relevant in connection with many topics, in fields ranging from physics to finance. It is a key ingredient to understanding phase-transitions in complex systems, both in and far-from thermal equilibrium. In particular, the quantity known as the escape time (or mean first passage time) from one stable state to another has found numerous applications in a variety of interesting and novel problems. For example, it plays a key role in stochastic resonance , in describing fluctuation-induced transport such as occurs in kink motion  and ratchets . Even the extent of chaos in Hamiltonian systems has been shown to have connections with this quantity . A nice collection of these and other stochastically driven processes can be found in Ref .

However, all of the above examples have been formulated within a standard Brownian framework, for which diffusion properties are normal. In this paper we look at the problem of calculating the escape time for systems exhibiting anomalous diffusion of the correlated type (in contrast to Levy type diffusion, which we do not discuss here). An understanding of escape time properties in such systems could open the door for understanding new stochastically driven phenomena. To our knowledge there has yet been little work done along these lines, although we are aware of some studies relating the anomalous transport properties on a random comb to the distribution of mean first passage times .

The systems we are interested in are such that the diffusion is dependent on the density of particles $`\rho `$, resulting in a diffusion coefficient which is proportional to a power ($`\nu 1`$) of $`\rho `$. Many physical systems are well-described by this class of processes. Let us mention, amongst other examples, percolation of gases through porous media ($`\nu 2`$), thin saturated regions in porous media ($`\nu =2`$) , gravitational spreading of thin liquid films ($`\nu =4`$) , heat transfer by Marshak waves ($`\nu =7`$) , surface growth ($`\nu =3`$) , spatial diffusion of biological populations ($`\nu 2`$), plasma flows ($`\nu <1`$). Explicitly, these processes are ruled by an equation of the type known in the literature as porous media equation

$$_t\rho (x,t)=D_x^2[\rho (x,t)]^\nu ,$$
(1)
where $`x`$ is a dimensionless coordinate representing a bond-length, angle or any other chemical or physical state variable, $`t`$ is the dimensionless time and $`\nu D>0`$. Rewritting the nonlinear term as $`_x(D\nu \rho ^{\nu 1}_x\rho )`$, it becomes evident that the restriction $`D\nu >0`$ guarantees that the flux will be from more dense to less dense regions.

Since the non-linearity in $`\rho `$ is known to lead to anomalous diffusion if $`\nu 1`$ (namely superdiffusion for $`\nu <1`$ and subdiffusion for $`\nu >1`$, as $`<x^2(t)>t^{\frac{2}{\nu +1}}`$) important anomalies are also expected when crossing over a barrier is involved. Precisely, we want to unveil here how escape properties are altered when $`\nu 1`$.

The paper is organized as follows. In Sec. II we present the systems of interest and discuss some of their general features. Because fluctuations are determined by $`\rho (x,t)`$ for $`\nu 1`$, the escape behavior will depend on the initial condition $`\rho (x,0)`$. Therefore we first consider systems in the vicinity of the steady state, a condition which allows analytical treatment. Numerical and analytical results for this case are presented in Secs. III and IV, respectively. In Sec. V we study numerically the escape behavior of systems far from the steady state, comparing the results with the previous ones. Finally, section VI contains concluding remarks.

## II The system

Let us consider a set of identical particles immersed into a thermal environment such as that described by the porous media equation (1). Under the influence of an external bistable potential $`U(x)`$, introduced in order to probe the escape behavior, the density of particles evolves following the nonlinear Fokker-Planck (FP) equation:

$$_t\rho (x,t)=_x[_xU(x)\rho (x,t)]+D_x^2[\rho (x,t)]^\nu .$$
(2)
This class of equations has been the object of diverse previous studies.

The stationary solution of Eq. (2) is

$$\rho _s(x)=[1(\nu 1)\beta V(x)]_+^{\frac{1}{\nu 1}}/Z,$$
(3)
where $`[f]_+=max\{f,0\}`$, $`Z`$ is a (positive) normalization constant, $`\beta =Z^{\nu 1}/(\nu D)`$ and $`V(x)=U(x)U_o`$, with $`U_o`$ the absolute minimum of the potential. In the limit $`\nu 1`$ the standard linear Fokker-Planck equation is obtained. In such a case, the steady state characterized by the Boltzmann-Gibbs distribution $`\rho _s(x)\mathrm{exp}(U(x)/D)`$, is recovered. However, for $`\nu 1`$, the stationary solutions of Eq. (2) have the form of the Maximum Tsallis Entropy probability distributions, as already discussed previously , even in the absence of external drift . It is worth recalling that phemonena such as full developed turbulence , the hadronic transverse moment distribution in high energy scattering process $`\text{e}^+\text{e}^{}\text{hadrons}`$ , among others, have been satisfactorily described in terms of distributions similar to (3) instead of the canonical stationary one.

Steady state solutions are illustrated in Fig. 1 for a quartic potential. Note that a cut-off condition (Tsallis cut-off), yielding regions with null probability, arises in the $`\nu >1`$ case (see Fig. 1b). For a quartic potential the condition $`\nu >3`$ must hold so that the solutions can be normalized. However, the free-particle case requires $`\nu >1`$ so we restrict our discussion to this regime.

The nonlinearity in the diffusion term of Eq. (2) accounts for the fact that the environment presents some kind of disorder or long range correlations in space-time leading to diffusion anomalies. The expression $`\beta =Z^{\nu 1}/(\nu D)`$ can be interpreted as a generalized Einstein relation for this scenario. Note that in disordered or correlated systems such as those discussed here, the standard Einstein relation is expected to be recovered in the absence of disorder . This corresponds to the case of $`\nu =1`$ yielding the well-known result $`D=1/\beta `$. Also, as was shown in , the time-dependent form of these Einstein relations can be used to demonstrate the anomalous scaling properties of these nonlinear diffusion systems. For the free particle one obtains $`<x^2(t)>1/\beta (t)Z^2(t)t^{\frac{2}{\nu +1}}`$.

The Ito-Langevin (IL) counterpart of Eq. (2) reads

$$\dot{x}=_xU(x)+\sqrt{|D|}[\rho (x,t)]^{\frac{\nu 1}{2}}\eta (t),$$
(4)
where $`\eta (t)`$ is a delta correlated Gaussian noise with zero mean and variance 2. In the particular case $`\nu =1`$, the standard Langevin equation for constant noise is recovered. It is noteworthy that this is a phenomenological description, in which the microscopic trajectories are determined by the macroscopic quantity $`\rho `$ when $`\nu 1`$. Physically, this represents a kind of statistical feedback. As with state-dependent noise, it is to be seen as the influence of the environment, which is otherwise not explicitly taken into account by the equations of motion. As a particle evolves, it interacts with the environment such that it reacts to the collective density of states around it. We can think of the subdiffusive case as a kind of ”attraction” to the other particles: Particles tend to stay close to the other particles, fluctuating not far from them. Conversely, we can think of super-diffusive cases as a kind of reaction to the sparseness: If the particle is in a highly populated region then it is in a sense confined by the other particles, and fluctuations are not so large, but as soon as it gets into less dense regions it does not experience this confinement and fluctuations can get very large.

## III Numerical results in the vicinity of the steady state

For numerical experiments we chose as prototype of double-well potential the quartic polynomial $`V(x)=ax^4+bx^3+cx^2+d`$. The coefficients were chosen as in Fig. 1, for which $`(x_L,x_O,x_R)=(0,1,3)`$, with $`x_L`$, $`x_O`$ and $`x_R`$ corresponding to the bottom of the left-hand well, the top of the barrier and the bottom of the right-hand well, respectively. We studied the escape behavior close to the steady state. That is, once a population of a large number of particles has already attained the steady state described by Eq. (3), a probe was injected at $`x_L`$. Then its trajectory was obtained by solving, following the numerical scheme in ref. , the IL (4) for $`\rho (x,t)=\rho _s(x)`$, starting from $`x(t=0)=x_L`$. Typical trajectories are displayed in Fig. 2. For $`\nu >1`$, fluctuations are reduced and trajectories result confined to the region within the cut-off boundaries (see also Fig. 1b); moreover, when the diffusion constant $`D`$ is smaller than a critical value $`D_c`$ (here $`D_c0.17`$ for $`\nu =2`$), the state space becomes disconnected and crossings become forbidden. For $`\nu <1`$, the amplitude of noise is enhanced in the regions of low density and the entire space tends to be populated.

We measured the mean first passage time, i.e., the average time interval $`T(x_Lx)`$ that a particle at $`x_L`$ takes to reach for the first time a given state $`x>x_L`$. In Fig. 3 we present plots of $`T(x)T(x_Lx)`$ vs. $`x`$. For $`\nu 1`$ (Fig. 3a), plateaux become evident as $`D`$ approaches $`D_c`$ indicating that most of the time is spent overcoming the barrier around $`x_O`$. On the other hand, for $`\nu <1`$ (Fig. 3b), the passage time is sensitive to the exact final state and there is not a well defined plateau, even in the small-$`D`$ regime. Moreover, as $`D`$ decreases, the curves collapse to a limiting one for states below $`x_R`$, but grow faster above $`x_R`$, diverging in the limit $`D0`$. The escape behavior seems to be discontinuous at $`D=0`$. In fact, for $`D=0`$ there is no diffusion, however, for finite $`D`$ the particle is attracted towards the deepest valley at $`x_R`$ and becomes trapped within a typical time interval which is bounded from above. This effect can be understood having in mind that fluctuations depend on $`D`$ not only through the factor $`\sqrt{|D|}`$ but also by means of the density through a factor that, for $`\nu <1`$, becomes very large outside the neighborhood of the absolute minimum where particles tend to concentrate as $`D0`$. In other words, the deterministic case is not recovered when $`D0`$ since the effective diffusion coefficient $`D\rho ^{\nu 1}`$ does not vanish in that limit due to the singularity at $`\rho =0`$.

## IV Analytical considerations

Let us show that these results can be understood analytically. For a system in the vicinity of the steady state, we can consider the following approximation for Eq. (2)

$$_t\rho (x,t)_x[_xU(x)\rho (x,t)]+D_x^2[\{\rho _s(x)\}^{\nu 1}\rho (x,t)].$$
(5)
Once the FP equation is linear, the problem of escape from a well can be treated directly, following the same lines as for homogeneous processes characterized by time independent drift and diffusion coefficients . Basically an equation for the probability that the particle is still within a given interval of state space at time $`t`$ is found using the corresponding backward Fokker-Planck equation and solved under appropriate boundary conditions. In this way, one finds that the mean first passage time $`T(x_1x_2)`$, for $`x_1<x_2`$, is given by

$`T(x_1x_2)`$ $`=`$ $`|\nu |\beta {\displaystyle _{x_1}^{x_2}}[1(\nu 1)\beta V(y)]_+^{\frac{|\nu |}{1\nu }}𝑑y`$ (6)
$`\times `$ $`{\displaystyle _{\mathrm{}}^y}[1(\nu 1)\beta V(z)]_+^{\frac{\mu }{\nu 1}}𝑑z,`$ (7)
where $`\mu =1`$ if $`\nu >0`$ and $`\mu =12\nu `$ if $`\nu <0`$. Expression (7) reproduces numerical experiments with excellent agreement as illustrated in Fig. 3.

In Fig. 4 we show $`TT(x_R)T(x_Lx_R)`$ as a function of $`1/D`$ (full lines), for different values of $`\nu >0`$, as calculated from Eq. (7). $`T`$ represents a measure of the escape time from the left to the right-hand well, even in the $`\nu <1`$ cases where plateaux are not well defined. In the range $`\nu >1`$, $`T`$ diverges at a value $`D_c`$, defined by the cut-off prescription, below which the right-hand well becomes inaccessible. In the $`0<\nu <1`$ case, $`T`$ saturates as $`1/D`$ increases. The hyperdiffusive regime $`\nu <0`$ (hence $`D<0`$), where spreading is faster than ballistic, demonstrates the same general features discussed for the region $`0<\nu <1`$ but $`|D|`$ must be considered instead of $`D`$. For any $`\nu `$ and small $`1/|D|`$ the escape time follows the power law $`T\beta ^{\frac{3}{4}}1/|D|^{\frac{3}{\nu +3}}`$.

If $`x_1x_L`$ and $`x_2x_R`$, then, it is possible to find an approximate expression for the escape time $`T`$ when $`|D|`$ (hence $`1/\beta `$) is sufficiently small, noting that the integrands in Eq. (7) present sharp peaks at $`x_O`$ and $`x_L`$ respectively. In that case the integrals can be evaluated by a saddle-point approximation extending the integration limits to the whole space. Following this procedure we arrive at

$$T\frac{2\pi }{\sqrt{\omega _L\omega _O}}\frac{2|\nu |}{|\nu |+\mu }\left(\frac{1(\nu 1)\beta V(x_O)}{1(\nu 1)\beta V(x_L)}\right)^{\frac{|\nu |+\mu }{2(1\nu )}},$$
(8)
where $`\omega _L`$ and $`\omega _O`$ are the frequencies at the bottom of the left well and at the top of the barrier, respectively. Expression (8) is a generalization of the Arrhenius law, which, as expected, is recovered in the limit $`\nu 1`$. In fact, in that limit, $`T(2\pi /\sqrt{\omega _L\omega _O})\mathrm{exp}(\mathrm{\Delta }V/D)`$, where $`\mathrm{\Delta }VV(x_O)V(x_L)`$ is the barrier height.

For comparison, the approximation given by Eq. (8) is also exhibited in Fig. 4 (dashed lines). The approximation is good for large $`1/|D|`$, as expected. It works better for $`\nu >1`$. Let us comment the main features revealed by this expression. When $`\nu >1`$, it foresees the divergence of $`T`$ at finite $`D`$. In fact, $`D_c`$ is obtained from $`1/\beta _c(\nu 1)V(x_O)`$. When $`\nu <1`$, saturation of $`T`$ for large $`1/|D|`$ is also predicted (unless $`V(x_L)=0`$) since $`\beta `$ is an unbounded increasing function of $`1/|D|`$. If $`V(x_L)=0`$, then Eq. (8) indicates that $`T`$ diverges for vanishing $`|D|`$. In particular, if $`0<\nu <1`$, $`T\beta ^{\frac{\nu +1}{2(1\nu )}}1/D^{\frac{1}{1\nu }}`$ and the deterministic limit is achieved. In the limit $`\nu 1`$ the exponential growth of $`T`$ with $`1/D`$ is always recovered.

## V Numerical results far from the steady state

The problem in the vicinity of the steady state actually corresponds to a linear one with a state dependent diffusion coefficient. However, it allows an analytical treatment which can be had in mind as a reference when studying more general cases. In order to test how the previous results compare to those of a more general situation, we also performed numerical studies of the escape properties far from the steady state. Particularly, we studied the case where particles are injected all at the same time at $`x_L`$. This instance requires simultaneous integration of the FP equation, in order to follow the evolution of $`\rho (x,t)`$ starting from $`\rho (x,0)=\delta (xx_L)`$, together with integration of the IL equation (4), starting from $`x(t=0)=x_L=0`$. Now, the parameter $`\nu `$ must lie in the region $`\nu >0`$ due to the divergence in Eq. (2). An implicit finite-difference scheme with centered space differences was employed for numerical integration of the nonlinear FP equation. The time evolution of the density is illustrated in Fig. 5.

The escape time $`T`$ as a function of $`1/D`$ (symbols) obtained for different values of $`\nu `$ was included in Fig. 4. Let us compare this case to the precedent steady one. For sufficiently large $`D`$, $`T`$ is not sensitively dependent on the initial distribution and Eq. (7) fits well to the numerical results for any $`\nu >0`$, following the power law $`T1/D^{\frac{3}{\nu +3}}`$ derived above. On the other hand, for small $`D`$, crossing times become closer to those of the standard case $`\nu =1`$ for any $`\nu `$. This can be understood as follows. For $`\nu >1`$, passage times are smaller than those given by Eq. (7) since, as the distribution evolves, there is an initial passage even between regions disconnected at the steady state (see Fig. 5(a)). However, our results suggest that the divergence of $`T`$ for a finite critical $`D`$, close to $`D_c`$, still occurs. On the other hand, in the range $`\nu <1`$, crossing times are larger than those given by Eq. (7) since now the density of particles is initially unfavorable for surmounting the barrier (see Fig. 5(b)). Saturation is not observed and the escape time increases with $`1/D`$ apparently following a power-law. It is worth noting that, as derived above, a power-law with exponent $`1/(1\nu )`$ is the one expected if the average effective potential felt by crossing particles has the absolute minimum at $`x_L`$ which is consistent with the observed density evolution (see Fig. 5(b)).

## VI Final remarks

Summarizing, we have obtained the escape time for systems exhibiting anomalous diffusion due to a stochastic nonlinear dependency on the density. For steady-state conditions, we obtain an analytical expression for the mean first passage time whose predictions are in excellent agreement with numerical results (Fig. 3). This analytical expression yields a generalization of Arrhenius law. A behavior quite different from that of the standard Brownian case $`\nu =1`$ is depicted. Under close to stationary conditions, two regimes are detected: In the region $`\nu <1`$ (superdiffusion), the escape time $`T`$ saturates for vanishing $`D`$, if $`V(x_L)0`$, and grows with $`1/D`$ following a power-law, otherwise. In the region $`\nu >1`$ (subdiffusion), $`T`$ diverges at $`D_c`$ (Fig. 4).

These results give hints on what should be expected in more general cases. For systems far from the steady state, $`T`$ grows with $`1/D`$ apparently following a power-law in the superdiffusive cases while $`T`$ diverges at finite $`D`$ in the subdiffusive ones.

Acknowledgements: We want to thank R. O. Vallejos for enlightening discussions. We acknowledge Brazilian Agencies, CNPq and FAPERJ for partial financial support.

## Captions for figures

Figure 1: The cut-off condition. (a) Dimensionless double-well potential $`V(x)=ax^4+bx^3+cx^2+d`$, with $`a=1/48,b=1/9,c=1/8,d=3/16`$. The stationary distribution $`\rho _s(x)`$ is shown for $`\nu =2`$ (b) and $`0.5`$ (c), for different values of $`D`$ as indicated on the figure. For $`\nu 1`$ the full state space is covered with power-law tails. For $`\nu >1`$ a cut-off restricts the attainable space. Observe in (b) that as $`D`$ decreases particles become more confined until only the neighborhood of the deepest valley is allowed. The horizontal lines in (a) represent the cut-off condition $`V(x)=1/\beta `$ which defines the allowed regions for $`\nu =2`$ and the same values of $`D`$ as in (b). All quantities are dimensionless.

Figure 2: Typical trajectories $`x`$ vs. $`t`$ for $`(\nu ,D)=(0.5,0.5)`$ (dark gray), $`(2,0.5)`$ (black) and (2,0.15) (light gray).

Figure 3: $`T(x)T(x_Lx)`$ vs. $`x`$ for different values of $`D`$ indicated on the figure and $`\nu =2`$ (a), $`0.5`$ (b). Circles correspond to numerical experiments (mean value over 1000 realizations) and full lines to theoretical prediction given by Eq. (7).

Figure 4: Escape time $`TT(x_R)`$ as a function of $`1/D`$, for different values of $`\nu >0`$ indicated on the figure. Full lines are generated from Eq. (7). Dashed lines correspond to the low-$`D`$ approximation given by Eq. (8). Symbols correspond to the initial condition where all the particles (at least 1000) are injected at the same time at $`x_L`$. Dotted lines are guides for symbols. Insert: Detail (semi-log) of the low-$`D`$ region for $`\nu 1`$.

Figure 5: Time evolution of the density of particles obtained by numerical integration of Eq. (2) with $`\rho (x,0)=\delta (x)`$ for ($`\nu `$,$`D`$)=(4.0,2.5) (a) and (0.5,0.1) (b). The profiles correspond to times $`t`$ indicated on the figure.