# Coalescence of limit cycles in the presence of noise

Sergei Shmakov<EMAIL_ADDRESS>Peter B. Littlewood
<EMAIL_ADDRESS>James Franck Institute and Department of Physics, The
University of Chicago, Chicago, Illinois 60637, United States

###### Abstract

Complex dynamical systems may exhibit multiple steady states, including time-
periodic limit cycles, where the final trajectory depends on initial
conditions. With tuning of parameters, limit cycles can proliferate or merge
at an exceptional point. Here we ask how dynamics in the vicinity of such a
bifurcation are influenced by noise. A pitchfork bifurcation can be used to
induce bifurcation behavior. We model a limit cycle with the normal form of
the Hopf oscillator, couple it to the pitchfork, and investigate the resulting
dynamical system in the presence of noise. We show that the generating
functional for the averages of the dynamical variables factorizes between the
pitchfork and the oscillator. The statistical properties of the pitchfork in
the presence of noise in its various regimes are investigated and a scaling
theory is developed for the correlation and response functions. The analysis
is done by perturbative calculations as well as numerical means. Finally,
observables illustrating the coupling of a system with a limit cycle to a
pitchfork are discussed and the phase-phase correlations are shown to exhibit
non-diffusive behavior with universal scaling.

††preprint: APS/123-QED

##  I Introduction

Dynamical systems exhibiting bifurcations - where there is a sudden change of
state - are of importance in many fields of study, ranging from biological
sciences to quantum systems [1, 2, 3]. For example (relevant to the present
work) the onset of polariton condensation and lasing corresponds to collective
synchronization of oscillators [4, 5, 6, 7], and the development of
superconductivity and superfluidity can be mathematically described in the
same way, perhaps most easily understood in the context of exactly solvable
mean field Richardson-Gaudin models [8]. In biology the proximity to an
instability can be used to tune the sensitivity of a sensor to an external
stimulus [9], seen in the tuning of hair bundles in the cochlea [10]. Many
physiological functions depend on the maintenance of synchronised patterns in
time (limit cycles), including breathing [11], cardiac function, and circadian
rhythms [12]. Photonic oscillators such as optically injected lasers also
exhibit limit cycles, allowing for unique tunability and response properties
[13, 14].

Bifurcations in spatially extended dynamical systems may fall into generalized
classes, and the onset of oscillatory behavior from a stationary state bears
analogy to a phase transition in a thermodynamic system. There is a particular
subclass of such transitions which is characterized by an exceptional point of
the dynamics, where two modes merge with identical eigenvalues [15, 16]. This
critical point will be the focus of our paper.

All physical systems operate in the presence of environmental fluctuations and
noise. This can mean that dynamical transitions in finite systems (or infinite
systems below the lower critical dimension) are blurred. However, dynamical
responses can still be tuned - the classic case is usually called stochastic
resonance [17, 18], where a bistable system in the presence of noise is found
to be maximally sensitive to a.c. fields whose frequency matches the typical
rate of hopping between minima. As an example of a bistable system (at least
in one regime), we take a pitchfork bifurcation, which can be coupled to more
general dynamical systems in order to imbue those with bifurcation behavior.
Unlike the pitchfork itself, these coupled systems can exhibit much broader
classes of dynamics, where the coalescing states can be fixed points, various
stable orbits such as limit cycles, and other attractors [19]. In this paper
we study such a dynamical system: noiseless dynamics can merge two limit
cycles (by parameter tuning) at an exceptional point. What are the response
and correlation functions of this system in the presence of noise?

Noisy limit cycles and oscillators are of considerable interest in many fields
of study, and have been explored in various contexts [20, 21, 22, 23]. In the
vicinity of the bifurcation, particular dynamical systems can be expressed in
their normal forms that are generic to a class of models. These normal forms
can therefore be studied to infer general information on how coupled dynamical
systems with bifurcations behave. In this paper we will investigate a Hopf
oscillator coupled to the pitchfork, which can describe a class of systems
with coalescing limit cycles. We will show that there are characteristic
features in the dynamics that mix the established scaling of stochastic
resonance with additional time scales associated with the recurrent limit
cycle oscillations.

One reason to study critical points is that they may have universal dynamical
signatures. The existence of a scaling theory for a model system (which we
will derive in this case) should, near the critical point, provide a
universality class for dynamical correlations of systems with the same
underlying symmetry. Consequently experimental systems that are
microscopically much more complicated than our model may nonetheless have
universal features in their correlation functions which reveal underlying
dynamical transitions. Our work will suggest methods to analyze data to
explore whether the bifurcation model is appropriate for the system at hand.

The nature of the noise depends on the details of the system, and in this
paper we investigate the simplest case: Gaussian white noise. One of the main
results of this paper, which describes how expectation values of operators
factorize between the pitchfork and the rest of the dynamical system, are
valid for any type of noise, as long as it is uncorrelated between the
pitchfork and the rest of the system. The behavior of the pitchfork with noise
has been investigated in various contexts and in various regimes, such as
perturbative methods [24] and stochastic resonance [17, 18]. However, to
authors’ knowledge, there is no treatment that combines the noisy pitchfork in
it’s various regimes into one theory. Here we investigate and collect results
for it in all regimes, and sew them together with scaling theory. That theory
relies on the observation that the pitchfork has a very limited number of
dimensionless parameters that the perturbative expansions can run over.

Hence, once the factorization result and the scaling theory of the pitchfork
are explored, our goal is to identify observables that can be measured to see
whether an experimental system at hand does indeed exhibit factorization and
whether it exhibits a bifurcation governed by the pitchfork. Other possible
venues of exploration for future study include investigating other types of
orbits that can bifurcate coupled to a pitchfork, other types of noise, and
also noise that is correlated between the pitchfork and the rest of the
dynamical system.

The reason to investigate an abstract model is that it may offer universal
pointers for physical behavior. Slow fluctuations between different conceptual
states are common in complex systems and are often loosely assumed in biology.
It is almost never the case that one will have a complete enough understanding
of the system in order to produce a detailed model. We rely instead on general
principles, which are necessarily mostly static. However, the statistics of
fluctuations near to the critical point can be deterministic of the underlying
model. We will show that fluctuations between limit cycles have a greater
richness than fluctuations between fixed points (traditional stochastic
resonance). Such correlations could be determined from data, agnostic from a
model.

This paper is organized as follows. In Sec. II we introduce the model system
of a Hopf oscillator coupled to a pitchfork bifurcation and demonstrate that
it can be factorised into noisy dynamics of the pitchfork and of the limit
cycle, eqn. (7). Here we also explore what that result implies for various
observables of the limit cycle, eqns. (9), (10), (11), (12). In Sec. III we
perform a scaling analysis of the noisy pitchfork (which is a classic model
for stochastic resonance). In Sec. IV we demonstrate the numerical results and
discuss how physical correlation functions of variables that correlate both
the periodic ”limit cycle” space and the ”pitchfork” space couple and lead to
novel interference features in the dynamics. In Sec. V we close with a
discussion, and in particular talk about potential experimental impacts.

##  II Theoretical Model

We take a Hopf oscillator as a model for a limit cycle and couple it to a
pitchfork bifurcation:

$\displaystyle\dot{\psi}=(a+i\Omega+b|\psi|^{2})\psi+ic\psi w^{n}\;,$ (1)
$\displaystyle\dot{w}=Rw-gw^{3}\;.$ (2)

Here $a,\Omega,b$ are real parameters chosen to produce a limit cycle, $c$ is
the coupling strength between the pitchfork and the oscillator and $R$ and
$g>0$ are parameters of the pitchfork. All parameters are of appropriate units
and are taken to be constant in time. We assume that the Hopf oscillator is in
a regime with a stable limit cycle ($a>0,b<0$), although the particular values
of the parameters will not be relevant for most of the analysis. This system
has been shown to exhibit bifurcation behaviour and for $n=1$ and $n=2$ there
is also an exceptional point at $R=0$ [19]. Notice that in this case the Hopf
oscillator is coupled to the pitchfork in a multiplicative way. Since $\psi$
is complex, we can break it up into its real and imaginary components. We then
add noise with the assumption that it is not correlated between the the
pitchfork and the oscillator.

$\displaystyle\dot{\psi}_{x}=a\psi_{x}-\Omega\psi_{y}+b|\psi|^{2}\psi_{x}-c\psi_{y}w^{n}+\sqrt{B_{\psi_{x}}(\psi)}\xi_{1}(t)\;,$
$\displaystyle\dot{\psi}_{y}=a\psi_{y}+\Omega\psi_{x}+b|\psi|^{2}\psi_{y}+c\psi_{x}w^{n}+\sqrt{B_{\psi_{y}}(\psi)}\xi_{2}(t)\;,$
$\displaystyle\dot{w}=Rw-gw^{3}+\sqrt{B_{w}(w)}\xi_{3}(t)\;,$

where the stochastic differential equations (SDEs) will be taken under Ito
discretization and the noise increments are defined to be independent Gaussian
white noises
$\left<\xi_{i}(t_{1})\xi_{j}(t_{2})\right>=\delta_{ij}\delta(t_{1}-t_{2})$. We
will assume that the noises along the orthogonal directions of the Hopf
oscillator are equal: $B_{\psi_{x}}=B_{\psi_{y}}=B_{\psi}$. Calling
$\psi_{x}=r\cos(\theta)$ and $\psi_{y}=r\sin(\theta)$ we can convert the Hopf
oscillator into polar coordinates, keeping the terms produced due to Ito
discretization:

$\displaystyle\dot{r}=ar+br^{3}+\frac{B_{\psi}}{r}+\sqrt{B_{\psi}}\xi_{r}(t)\;,$
(3)
$\displaystyle\dot{\theta}=\Omega+cw^{n}+\frac{\sqrt{B_{\psi}}}{r}\xi_{\theta}(t)\;,$
(4)

where $\xi_{r}(t)$ and $\xi_{\theta}(t)$ are defined through an orthogonal
transformation from the original $\xi_{1}(t)$ and $\xi_{2}(t)$, and are also
Gaussian white uncorrelated noises [25]. We can see here that the coupling to
the pitchfork only affects the dynamics of $\theta$, as in the corresponding
classical case, and in fact the dynamics of $r$ are completely independent of
both $\theta$ and $w$.

We can now introduce the Martin-Siggia-Rose response fields $w^{q}$, $r^{q}$
and $\theta^{q}$ and apply the Janssen-De Dominicis-Peliti formalism [26, 27,
28, 29] to express the problem in terms of path integrals:

$\begin{split}\mathcal{Z}=\int\mathcal{D}[w,w^{q},\theta,\theta^{q},r,r^{q}]&\exp\left[\int
dt\,iw^{q}(\dot{w}-Rw+gw^{3})-\frac{B_{w}}{2}w^{q2}\right]\exp\left[\int
dt\,i\theta^{q}(\dot{\theta}-\Omega-
cw^{n})-\frac{B_{\psi}}{2r^{2}}\theta^{q2}\right]\\\
&*\exp\left[\int\,dtir^{q}(\dot{r}-f(r))-\frac{B_{\psi}}{2}r^{q2}\right]\;,\end{split}$
(5)

where $\mathcal{Z}=1$ identically, but can be used as a generating functional
for operators comprised of our dynamic variables $r,\theta,w$. Here
$f(r)=ar+br^{3}+B_{\psi}/r$ contains the deterministic part of the evolution
of $r$. We will assume unspecified initial conditions, in other words the
systems were settled into their equilibrium distributions at $t=-\infty$. Here
we also note that we have not taken into account the compactness of $\theta$,
which in principle should be handled. However, the following result still
holds and is confirmed numerically. We can now demonstrate that expectation
values of these operators factorize into expectation values of the
constituents corresponding to field theories of the Hopf oscillator and
pitchfork on their own. That fact relies on the linearity of the equation for
$\theta$. It can be clearly seen that when calculating an expectation value of
any operator $w^{n}(w^{q})^{m}$, we can integrate out the fields corresponding
to the Hopf oscillator and we’ll just be left with the dynamics of the
pitchfork. Therefore, the complexity comes from any parts of the operators of
the form $\theta^{n}(\theta^{q})^{m}r^{k}(r^{q})^{l}$. These operators can be
obtained from a generating functional:

$\begin{split}&\mathcal{Z}[\lambda_{\theta},\lambda^{q}_{\theta},\lambda_{r},\lambda^{q}_{r}]=\\\
&\left<\exp\left[\int
dt\lambda_{\theta}\theta+\lambda_{\theta}^{q}\theta^{q}+\lambda_{r}r+\lambda_{r}^{q}r^{q}\right]\right>\;,\end{split}$
(6)

where all $\lambda$’s are functions of time. Employing the linearity of the
equation for $\theta$ it can be shown by integrating out $\theta$ and
$\theta^{q}$ that this generating functional factorizes into the generating
functional for the Hopf oscillator in the absence of the coupling and an
expectation value of an appropriate operator of the pitchfork only:

$\begin{split}&\mathcal{Z}[\lambda_{\theta},\lambda^{q}_{\theta},\lambda_{r},\lambda^{q}_{r}]=\mathcal{Z}_{\text{Hopf}}[\lambda_{\theta},\lambda^{q}_{\theta},\lambda_{r},\lambda^{q}_{r}]*\\\
&\left<\exp\left[\int
dt\,cw^{n}\int_{t}^{\infty}ds\,\lambda_{\theta}(s)\right]\right>_{w}\;,\end{split}$
(7)

where $\left<...\right>_{w}$ is the expectation value with respect to the
pitchfork action:

$S_{w}[w,w^{q}]=\int
dt\,\left(iw^{q}(\dot{w}-Rw+gw^{3})-\frac{B_{w}}{2}w^{q2}\right)\;.$ (8)

Hopf oscillators in their normal form have been studied in many contexts, with
varying mathematical depths and in application to a variety of subjects, both
with and without noise [30, 31, 32, 33, 34]. Thus, we won’t address part of
the generating functional $Z_{Hopf}$ that comes from the dynamics of the Hopf
oscillator alone. However, we can look at some quantities of interest of the
model and see how the coupling to the pitchfork affects them. For example, for
$n=1$ the connected angle correlation function for the full system can be
shown to be:

$\begin{split}&\left<\theta(t_{1})\theta(t_{2})\right>_{c}=\frac{\delta}{\delta\lambda_{\theta}(t_{1})}\frac{\delta}{\delta\lambda_{\theta}(t_{2})}\ln\mathcal{Z}[\lambda_{\theta},\lambda^{q}_{\theta},\lambda_{r},\lambda^{q}_{r}]\bigg{|}_{\lambda=0}\\\
&=\left<\theta(t_{1})\theta(t_{2})\right>_{c,Hopf}+c^{2}\int^{t_{1}}dt\int^{t_{2}}dt^{\prime}\left<w(t)w(t^{\prime})\right>_{c,w}\;.\end{split}$
(9)

On the other hand, the phase response function of the coupled system is
unchanged by the presence of the pitchfork:

$\begin{split}&\left<\theta(t_{1})\theta^{q}(t_{2})\right>_{c}=\frac{\delta}{\delta\lambda_{\theta}(t_{1})}\frac{\delta}{\delta\lambda^{q}_{\theta}(t_{2})}\ln\mathcal{Z}[\lambda_{\theta},\lambda^{q}_{\theta},\lambda_{r},\lambda^{q}_{r}]\bigg{|}_{\lambda=0}\\\
&=\left<\theta(t_{1})\theta^{q}(t_{2})\right>_{c,Hopf}\;.\end{split}$ (10)

More importantly, since the object of interest is a limit cycle, the
observable of interest will be phase-phase correlations, which, due of the
factorization discussed above, will take the form:

$\begin{split}&\left<\exp\left[i\theta(t_{1})-i\theta(t_{2})\right]\right>=\\\
&\left<\exp\left[i\theta(t_{1})-i\theta(t_{2})\right]\right>_{Hopf}\left<\exp\left[-ic\int_{t_{1}}^{t_{2}}dtw^{n}(t)\right]\right>_{w}\;.\end{split}$
(11)

Another object of interest might be the phase diffusion constant [20], which
we will consider in the absence of external fields, leaving it time-
translation invariant:

$\begin{split}D&\equiv\int_{-\infty}^{\infty}d\tau\left<\left(\dot{\theta}(t)-\left<\dot{\theta}\right>(t)\right)\left(\dot{\theta}(t+\tau)-\left<\dot{\theta}\right>(t+\tau)\right)\right>\\\
&=\int_{-\infty}^{\infty}d\tau
B_{\psi}\left<\frac{\delta(\tau)}{r(t)r(t+\tau)}\right>\\\
&+c^{2}\left(\left<w^{n}(t)w^{n}(t+\tau)\right>-\left<w^{n}(t)\right>\left<w^{n}(t+\tau)\right>\right)\;,\end{split}$
(12)

where we used the independence of statistics of $\xi_{\theta}$, $w$ and $r$ to
arrive at the second equality. For $n=1$, the diffusion constant is directly
dependent of the connected two-point function of the pitchfork. For $n=2$ it
will depend on the variance of the pitchfork distribution (which is just the
two-point function at zero time separation) as well as a particular form of
the four-point function. We can therefore see that to answer general questions
about the coupled system, it is useful to understand the dynamics of the
stochastic pitchfork bifurcation on it’s own, and it’s correlation functions.

##  III Pitchfork Scaling Theory

The noisy, forced pitchfork obeys the stochastic differential equation (SDE):

$\dot{w}=Rw-gw^{3}+h(t)+\sqrt{B_{w}}\xi_{3}(t)\;,$ (13)

and the bifurcation occurs at $R=0$ in the classical system. It’s not possible
to obtain the correlation and response functions for the pitchfork in all
regimes exactly, so we instead employ a scaling theory and patch together
different regimes. In the SDE, we assume that $g>0$ and $B_{w}\geq 0$. The
parameters and the functions of interest all have dimensions restricted by the
dimensions of $w$. We can therefore define:

$\displaystyle\left<w\right>$ $\displaystyle\equiv
M(R,g,B_{w},h)=\left(\frac{B_{w}}{g}\right)^{\frac{1}{4}}F_{1}\left(\frac{R}{B_{w}^{\frac{1}{2}}g^{\frac{1}{2}}},\frac{h}{B_{w}^{\frac{3}{4}}g^{\frac{1}{4}}}\right)\;,$
(14a) $\displaystyle\left<w(\Delta t)w(0)\right>_{c}$ $\displaystyle\equiv
C(R,g,B_{w},h,\Delta
t)=\left(\frac{B_{w}}{g}\right)^{\frac{1}{2}}F_{2}\left(\frac{R}{B_{w}^{\frac{1}{2}}g^{\frac{1}{2}}},\frac{h}{B_{w}^{\frac{3}{4}}g^{\frac{1}{4}}},g^{\frac{1}{2}}B_{w}^{\frac{1}{2}}\Delta
t\right)\;,$ (14b) $\displaystyle\left<w(\Delta t)w^{q}(0)\right>_{c}$
$\displaystyle\equiv\chi(R,g,B_{w},h,\Delta t)=i\Theta(\Delta
t)F_{3}\left(\frac{R}{B_{w}^{\frac{1}{2}}g^{\frac{1}{2}}},\frac{h}{B_{w}^{\frac{3}{4}}g^{\frac{1}{4}}},g^{\frac{1}{2}}B_{w}^{\frac{1}{2}}\Delta
t\right)\;,$ (14c)

where $F_{1},F_{2},F_{3}$ are dimensionless functions of dimensionless
variables. Similar scaling theories have been defined for the pitchfork,
especially in the context of stochastic resonance [18], but here we proceed to
define it for completeness. By $\left<...\right>_{c}$ we mean the connected
correlation function, or cumulant. Several assumptions were made. First, we
assume unconstrained initial conditions, so the system is time-translation
invariant (at least in the absence of forcing). Second, all other
dimensionless parameters that might arise from parameters implicit in the
definition of $h(t)$ are present in these functions, but are not explicitly
written for brevity. Third, causality is assumed for the response function
$\chi$ [24, 35], and the Heaviside function $\Theta(t)$ is defined in the
usual manner with the prescription that $\Theta(0)=0$ to be consistent with
Ito discretization [35].

These functions are universal for any collection of parameters. Suppose some
crossover is found that is intrinsic to these functions, for example in
$F_{2}(x,0,z)$ for some $x$. The universality of these functions then
guarantees that in physical parameters this crossover or transition has to
behave as $R\sim\sqrt{gB_{w}}$, with the exact proportionality constant to be
determined. The only way for this relationship to be broken is through the
prefactors that give these functions proper dimensions to match the physical
quantities of interest. However, these prefactors are also defined
universally, and thus their presence can be integrated generically into
understanding of such crossovers.

We can now extract general information about the functions
$F_{1},F_{2},F_{3}$. No attempt was made to study these functions for all
values of the parameters of the system. However, various regimes allow for use
of different techniques to compute these functions. We define the effective
bifurcation parameter $\tilde{R}=R/\sqrt{gB_{w}}$. For very large and negative
$\tilde{R}$, it is possible to apply the usual field-theoretical perturbative
techniques [35, 24, 36, 29, 28, 37]. For large and positive $\tilde{R}$, the
system is going to hop between the two wells so that we cannot just expand the
field theory around one of the wells, but rarely enough so that the system has
enough time to settle into either well before hopping again. In other words,
the relaxation around either stable state will be negligible for long-term
dynamics. Here it is possible to approach the study of the functions of
interest from the perspective of Kramer’s problem and stochastic resonance
[29, 17]. There will be two relevant time scales. The first one is responsible
for the decay within either well, and can be obtained by expanding the action
around either well. The lowest order approximation of this time scale is
$\left(2R\right)^{-1}$. The other one is responsible for the rate of hopping
between the wells and can be obtained through standard Kramers problem
calculation to be $\tau\approx\sqrt{2}\pi
R^{-1}\exp\left[R^{2}/(2gB_{w})\right]$. The crossover between the two is
therefore governed by the competition of the two time scales. Unsurprisingly,
perhaps, the ratio of these time scales is also a function of the effective
bifurcation parameter: $\pi/\sqrt{2}\exp\left(\tilde{R}^{2}/2\right)$. The
final regime centered around $R=0$ is mainly studied numerically. We also
define the effective time $\tilde{t}=\sqrt{gB_{w}}t$ and the effective
external field $\tilde{h}=h/(B_{w}^{\frac{3}{4}}g^{\frac{1}{4}})$.

##  IV Numerical results

###  IV.1 Pitchfork

To demonstrate and support the results of the scaling theory, we turn to
numerical simulations of our random process. First, we consider the pitchfork
itself. In Fig. 1, we simulate the process and plot the non-dimensionalized
correlation function $F_{2}(\tilde{R},0,\tilde{t})$ for different values of
$\tilde{R}$. Each line on the plot is an overlap of five results of
simulations (in different colors) at different values of $R$, $B_{w}$ and $g$
that leave $\tilde{R}$ invariant. Their complete overlap demonstrates the
validity of the developed scaling theory.

Figure 1:  Numerical simulation of $F_{2}(\tilde{R},0,\tilde{t})$ as a
function of $\tilde{t}$ at various $\tilde{R}$. Each line is an overlap of 5
simulations, in different colors, with different values of the individual
parameters that leave $\tilde{R}$ the same. Figure 2:  Numerical simulation of
$F_{2}(\tilde{R},0,\tilde{t})$ as a function of $\tilde{R}$ at various
$\tilde{t}$.

We expect that the decay of the correlation function will be governed by time-
scales that are described above. In Fig. 2 we plot the time slices of
$F_{2}(\tilde{R},0,\tilde{t})$ as a function of $\tilde{R}$ for different
values of $\tilde{t}$. As can be seen, there is a smaller time scale
associated with the immediate decay of the correlation function, which
perturbatively we expect to be $~{}R^{-1}$ in the $\tilde{R}<<0$ regime and
$~{}(2R)^{-1}$ in the $\tilde{R}>>0$ regime. In the $R\leq 0$ regime that is
the only time scale. In $R>0$ regime there is also Kramers time scale which is
exponential and controls the decay starting with values of smaller and smaller
$\tilde{t}$ as $\tilde{R}$ increases.

###  IV.2 Phase-phase correlations

Now we can return back to the original system. Since the information about the
period and the diffusion due to noise of the limit cycle is contained in the
phase-phase correlations, the operator of interest is the one in Eq. (11).
This operator can also be subjected to scaling theory, and with an additional
parameter - the coupling constant $c$, there is an additional dimensionless
parameter:
$\tilde{c}=cB_{w}^{\frac{n}{4}-\frac{1}{2}}g^{-\frac{n}{4}-\frac{1}{2}}$. The
phase correlations factor coming from the pitchfork is dimensionless,
therefore we call this factor
$F_{p}(\tilde{R},\tilde{h},t_{2,eff}-t_{1,eff},\tilde{c})$. The scaling is
demonstrated in Fig. 3 for $n=1$ and in Fig. 4 for $n=2$, where we plot
$F_{p}$ for various values of $\tilde{R}$ and $\tilde{c}$ as a function of
$\tilde{t}$ (and there is no external field). With this, the two cases of
$n=1$ and $n=2$ can be studied separately.

Figure 3:  Numerical simulation of
$Re[F_{p}(\tilde{R},0,\tilde{t},\tilde{c})]$ for $n=1$ as a function of
$\tilde{t}$ at various $\tilde{R}$ and $\tilde{c}$. Each line is an overlap of
5 simulations, in different colors, with different values of individual
parameters that leave the dimensionless variables the same. Figure 4:
Numerical simulation of $Re[F_{p}(\tilde{R},0,\tilde{t},\tilde{c})]$ and
$Im[F_{p}(\tilde{R},0,\tilde{t},\tilde{c})]$ for $n=2$ as a function of
$\tilde{t}$ at various $\tilde{R}$ and $\tilde{c}$. Each line is an overlap of
5 simulations, in different colors, with different values of individual
parameters that leave the dimensionless variables the same.

###  IV.3 $n=1$

First, by formal expansion it can be seen that for $n=1$, all the imaginary
terms will have an odd power of $w$. In the absence of external fields the
pitchfork SDE has a $Z_{2}$ symmetry, and thus they will all be zero
identically, and only the real part of the phase-phase correlations will
survive. To formally see that the dimensionless parameter above (which is
$cB_{w}^{-\frac{1}{4}}g^{-\frac{3}{4}}$ for $n=1$) is indeed what controls the
expansion, we can employ the scaling theory to describe the lowest order
terms:

$\begin{split}&F_{p}=\left<\exp\left[-ic\int_{t_{1}}^{t_{2}}dtw(t)\right]\right>_{w}\approx
1-\\\
&\frac{c^{2}}{2g^{\frac{3}{2}}B_{w}^{\frac{1}{2}}}\int_{0}^{\tilde{t}_{2}-\tilde{t}_{1}}dz^{\prime}\int_{-z^{\prime}}^{\tilde{t}_{2}-\tilde{t}_{1}-z^{\prime}}dzF_{2}(\tilde{R},\tilde{h},z)\;.\end{split}$
(15)

As an example, the following limit of this function can be calculated:

$\lim_{\tilde{R}\rightarrow-\infty}F_{p}(\tilde{R},0,\tilde{t},\tilde{c})=\exp\left[-\tilde{c}^{2}\left(\frac{\tilde{t}}{2\tilde{R}^{2}}-\frac{e^{\tilde{R}\tilde{t}}}{2\tilde{R}^{3}}+\frac{1}{2\tilde{R}^{3}}\right)\right]\;.$
(16)

The leading corrections to $\tilde{R}$ can be calculated perturbatively [24].

For $n=1$ the contribution of the pitchfork to the diffusion constant is also
conceptually quite simple:

$\begin{split}D_{p,n=1}&=c^{2}\int_{-\infty}^{\infty}d\tau\left(\left<w(t)w(t+\tau)\right>-\left<w(t)\right>\left<w(t+\tau)\right>\right)\\\
&=\frac{c^{2}}{g}\int_{-\infty}^{\infty}dz\,F_{2}(\tilde{R},\tilde{h},z)\;,\end{split}\;,$
(17)

Where we used the definition of the dimensionless correlation function $F_{2}$
and non-dimensionalized the integral measure $z=\sqrt{gB_{w}}\tau$. Although
we assumed time translation invariance, this expression demonstrates the
scaling of the diffusion constant in the presence of forcing as well. The
numerical results for the integral are demonstrated in Fig. 5 for various
values of $\tilde{R}$.

Figure 5:  Numerical simulation of $\int_{-\infty}^{\infty}dzF_{2}$ as a
function of $\tilde{R}$ on a logarithmic scale at $\tilde{h}=0$.

###  IV.4 $n=2$

The $n=2$ case is simpler conceptually. The coupling to the pitchfork changes
the period of the limit cycle governed, to lowest order, by the variance of
the pitchfork distribution. Since it is always positive, there are always
oscillations in the phase-phase correlations, as well as modulation of the
period. The oscillations can be understood intuitively: whatever the state of
the pitchfork, it always drives the limit cycle along the positive phase
direction, so at some points it is likely that on average the phase will be
opposite to the starting phase, hence the correlations get a negative real
part.

The effective coupling constant for the $n=2$ case is $\tilde{c}=cg^{-1}$,
which is independent of pitchfork noise. Therefore, changing the noise
strength of the pitchfork amounts to changing $\tilde{R}$ and thus, with
corresponding tuning of $R$, can be scaled away (with the time axis being
stretched as well).

###  IV.5 Frequency space

A useful representation of the phase-phase correlations that is experimentally
viable given the access to the oscillator and demonstrates the crossover is
the Fourier representation. In the absence of external fields acting on the
pitchfork (so that the time translation invariance is preserved), the Fourier
transform of phase-phase correlations takes the form:

$\begin{split}&\int
dt\left<\exp\left[-ic\int_{0}^{t}dtw^{n}(t)\right]\right>_{w}\exp(-i\omega
t)\\\ &=\int
dt\left<\exp\left[-ic\int_{0}^{t}dt\left(w^{n}(t)+\frac{\omega}{c}\right)\right]\right>_{w}\;.\end{split}$
(18)

An example of Fourier transforms for $n=1$ is given in Fig. 8. The bifurcation
in the pitchfork is clearly exhibited by the phase-phase correlations in
frequency space. It is clear that the dominant contributions to the transform
will come from the most likely states that the pitchfork can be found in. To
demonstrate it further, we plot in Fig. 6 the frequencies of the peaks of the
Fourier distribution as a function of $\tilde{R}$, and see that they have the
same bifurcation diagram as the pitchfork.

Figure 6:  Peak frequencies (blue dots) of the absolute value of the Fourier
transform of the phase-phase correlations as a function of $\tilde{R}$ for
$n=1$ and $\tilde{c}=1$. The bars represent the widths of half maxima of the
peaks. In red are the mean field stable states of the pitchfork. The
frequencies take discrete values due to the finite number of time steps in the
simulations of the Ito processes. The frequency is offset by the natural
frequency $\Omega$ of the oscillator.

For $n=1$, the mean field stable states of the pitchfork for $R>0$ are
$w=\pm\sqrt{R/g}$. Non-dimensionalizing the frequency of the Fourier transform
as $\tilde{\omega}=\omega/\sqrt{gB_{w}}$, we get the scaling for the peak
locations for $R>0$: $\tilde{\omega}_{peaks}=\pm\tilde{c}\sqrt{\tilde{R}}$.

For the case $n=2$, the expression for the mean field positions of the peaks
takes the form $\tilde{\omega}_{peaks}=-\tilde{c}\tilde{R}$, with the
replacement of $\tilde{c}$ by the corresponding definition at $n=2$ (see
above).

Figure 7:  Peak frequencies (blue dots) of the absolute value of the Fourier
transform of the phase-phase correlations as a function of $\tilde{R}$ for
$n=2$ and $\tilde{c}=1$. The bars represent the widths of half maxima of the
peaks. The dashed red line are the mean field stable states of the pitchfork.
The black line is $-\tilde{c}F_{2}(\tilde{R},0,0)$. As before, the frequencies
take discrete values due to the finite number of time steps in the simulations
of the Ito processes. The frequency is offset by the natural frequency
$\Omega$ of the oscillator.

Notice that the peaks will only appear in the negative frequencies. However,
the dominant contribution to the Fourier transform will come from the factor
$\left<w(t)^{2}\right>$, which means the peak locations will approximately
obey $\tilde{\omega}_{peaks}=-\tilde{c}F_{2}(\tilde{R},0,0)$. These
relationships are demonstrated in Fig. 7.

The crossover in this case is not as clear as in the $n=1$ case. However, in
both cases the widths of the peaks undergo a broadening around $R=0$ and are
controlled by $\tilde{c}$. This can be used as an indication of a crossover.

As a final remark, the phase-phase correlations also have a term coming from
the Hopf oscillator itself that would go into the Fourier transform. However,
in the limit cycle regime, the dynamics of the oscillator are simple, in that
the Fourier transform of the phase-phase oscillations has a peak at $-\Omega$,
the width of which depends on the noise strength. For the coupled system, it
means that in the above discussion the zero frequency will be shifted to
$-\Omega$, and the peaks widths will in reality be broader. However, their
locations and widths will undergo the same behavior described above, with the
same scaling laws.

Figure 8:  Example of absolute values of Fourier transforms of the phase-phase
correlations for $n=1$ and $\tilde{c}=1$ versus dimensionless frequency
$\tilde{\omega}=\omega/\sqrt{gB_{w}}$. The vertical axis is rescaled for each
plot, and is thus measured in arbitrary units (a.u.). a) $\tilde{R}=-3$, b)
$\tilde{R}=0$, c) $\tilde{R}=4$

##  V Conclusion

We have taken a model of bifurcating limit cycles that we believe to be
generally applicable to a wide variety of systems near a bifurcation and
explored it in the presence of noise. A pitchfork bifurcation is used to imbue
a limit cycle with a bifurcating behaviour and the main result of the study is
that this model system can be factorized into corresponding noisy dynamics of
the pitchfork and the oscillator.

The factorization of the generating functional depends only on the statistical
independence between the noise acting on the pitchfork and the noise acting on
the limit cycle. Therefore, specifics of the noise spectrum of either part of
the system only determine the details of their individual effects. For
example, a temporally correlated noise added to the pitchfork will add an
additional dimensionless time scale to the scaling theory, therefore changing
the above analysis. However, as the noise strength decreases, the scaling
functions are driven further into their asymptotic regimes, where their
behavior is easier to obtain analytically, and is also believed to be less
sensitive to changes in the nature of noise and to presence of external
fields.

The effect of the pitchfork with Gaussian white noise on the oscillator was
then explored, and a scaling theory was applied to demonstrate how the
parameters of the pitchfork must come into the functions of interest of the
limit cycle. The scaling theory is supported numerically, and measurable
effects of the pitchfork on the limit cycle, such as on the phase diffusion
constant and on the spectral resolution of the phase-phase correlations, are
demonstrated.

Our model is explicitly low-dimensional but we believe it may be a starting
point for a more general approach to extended systems. The transition between
silent (but excitable) and spontaneous activity has long been studied in
neural cultures (e.g. [38, 39, 40]), and the transition may be accompanied by
large scale activity that has been characterised as “avalanches” with a
1/f-like power spectrum. This is well-known as an outcome for dynamical
systems that are driven very slowly close to a threshold, named self-organised
criticality [41], and also for extended disordered nonlinear systems under
periodic driving [42, 43]. Such states are highly correlated in space and
time, but not fully periodic. Such systems are based directly on
generalisations of the low dimensional model we discuss here, though
extensions are not simple.

Close to the bifurcation, our model could potentially be applicable to a wide
variety of other systems, such as visual perception bistability [44] and
auditory perception bistability [45, 46], with suitable choices of bifurcation
parameters. Scaling ideas similar to our discussion have also been explored in
the context of multistable perception [47]. Although direct access to the
limit cycle dynamics of such systems might not be readily available, the
processing of stimuli is often inherently oscillatory (i.e. neural activity
[48]), thus driving the choice to explore limit cycles. Alternatively, in
perception research, the stimulus itself can be oscillatory [49]. Bifurcations
of limit cycles and bistable regimes also occur in queuing theory [50], among
others.

The choice of the type of bifurcation other than the pitchfork that can be
coupled to an oscillator would be based on the phenomenological knowledge of
each particular system. Same also applies to the type of noise that should be
used in the model, as well as it’s spectrum. The exploration of these ideas,
in addition to increasing the dimensionality of these models is reserved for
future work.

## Acknowledgments

The authors would like to thank Cheyne Weis and Ed Awh for helpful
discussions.

## References

  * Strogatz [2015] S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With applications to physics, biology, chemistry, and engineering_ (CRC Press, 2015).
  * Guckenheimer and Holmes [2002] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of Vector Fields_ (Springer, 2002).
  * Kuramoto [2003] Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Dover Publ, 2003).
  * Dicke [1954] R. H. Dicke, Coherence in Spontaneous Radiation Processes, Phys. Rev. 93, 99 (1954).
  * Rice and Carmichael [1994] P. R. Rice and H. J. Carmichael, Photon statistics of a cavity-QED laser: A comment on the laser–phase-transition analogy, Phys. Rev. A 50, 4318 (1994).
  * Eastham and Littlewood [2001] P. R. Eastham and P. B. Littlewood, Bose condensation of cavity polaritons beyond the linear regime: The thermal equilibrium of a model microcavity, Phys. Rev. B 64, 235101 (2001).
  * Kirton and Keeling [2018] P. Kirton and J. Keeling, Superradiant and lasing states in driven-dissipative Dicke models, New Journal of Physics 20, 015009 (2018).
  * Dukelsky  _et al._ [2004] J. Dukelsky, S. Pittel, and G. Sierra, Colloquium: Exactly solvable Richardson-Gaudin models for many-body quantum systems, Rev. Mod. Phys. 76, 643 (2004).
  * Winfree [2001] A. T. Winfree,  _The Geometry of Biological Time_ (Springer New York, NY, 2001).
  * Camalet  _et al._ [2000] S. Camalet, T. Duke, F. Jülicher, and J. Prost, Auditory sensitivity provided by self-tuned critical oscillations of hair cells., Proc. Natl. Acad. Sci. 97, 3183 (2000).
  * Smith  _et al._ [1991] J. C. Smith, H. H. Ellenberger, K. Ballanyi, D. W. Richter, and J. L. Feldman, Pre-Bötzinger complex: a brainstem region that may generate respiratory rhythm in mammals, Science 254, 726 (1991).
  * Rijo-Ferreira and Takahashi [2019] F. Rijo-Ferreira and J. S. Takahashi, Genomics of circadian rhythms in health and disease, Genome Med. 11, 82 (2019).
  * Himona  _et al._ [2022] G. Himona, V. Kovanis, and Y. Kominis, Isochrons, phase response and synchronization dynamics of tunable photonic oscillators, Phys. Rev. Res. 4, 10.1103/physrevresearch.4.l012039 (2022).
  * Valagiannopoulos and Kovanis [2021] C. Valagiannopoulos and V. Kovanis, Injection-Locked Photonic Oscillators: Legacy results and future applications, IEEE Antennas Propag. Mag. 63, 51 (2021).
  * Hanai  _et al._ [2019] R. Hanai, A. Edelman, Y. Ohashi, and P. B. Littlewood, Non-Hermitian Phase Transition from a Polariton Bose-Einstein Condensate to a Photon Laser, Phys. Rev. Lett. 122, 185301 (2019).
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Non-reciprocal phase transitions, Nature 592, 363 (2021).
  * McNamara and Wiesenfeld [1989] B. McNamara and K. Wiesenfeld, Theory of Stochastic Resonance, Phys. Rev. A 39, 4854 (1989).
  * Gammaitoni  _et al._ [1998] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Stochastic resonance, Rev. Mod. Phys. 70, 223–287 (1998).
  * Weis  _et al._ [2022] C. Weis, M. Fruchart, R. Hanai, K. Kawagoe, P. B. Littlewood, and V. Vitelli, Coalescence of Attractors: Exceptional points in non-linear dynamical systems, arXiv e-prints , arXiv:2207.11667 (2022).
  * Goldobin  _et al._ [2010] D. S. Goldobin, J.-n. Teramae, H. Nakao, and G. B. Ermentrout, Dynamics of Limit-Cycle Oscillators Subject to General Noise, Phys. Rev. Lett. 105, 154101 (2010).
  * Rosenblum and Pikovsky [2004] M. G. Rosenblum and A. S. Pikovsky, Controlling Synchronization in an Ensemble of Globally Coupled Oscillators, Phys. Rev. Lett. 92, 114102 (2004).
  * Grenfell  _et al._ [1998] B. T. Grenfell, K. Wilson, B. F. Finkenstädt, T. N. Coulson, S. Murray, S. D. Albon, J. M. Pemberton, T. H. Clutton-Brock, and M. J. Crawley, Noise and determinism in Synchronized Sheep Dynamics, Nature 394, 674–677 (1998).
  * Yoshimura and Arai [2008] K. Yoshimura and K. Arai, Phase Reduction of Stochastic Limit Cycle Oscillators, Phys. Rev. Lett. 101, 154101 (2008).
  * Hertz  _et al._ [2016] J. Hertz, Y. Roudi, and P. Sollich, Path Integral Methods for the Dynamics of Stochastic and Disordered Systems, J. Phys. A: Math. Theor. 50, 033001 (2016).
  * Gardiner [1997] C. W. Gardiner,  _Handbook of Stochastic Methods_ (Springer-Verlag, 1997).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical Dynamics of Classical Systems, Phys. Rev. A 8, 423–437 (1973).
  * De Dominicis and Peliti [1978] C. De Dominicis and L. Peliti, Field-theory renormalization and critical dynamics above $T_{c}$: Helium, antiferromagnets, and liquid-gas systems, Phys. Rev. B 18, 353 (1978).
  * Atland and Simons [2010] A. Atland and B. Simons,  _Condensed Matter Field Theory_ (Cambridge University Press, 2010).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Faber and Bozovic [2019a] J. Faber and D. Bozovic, Noise-induced chaos and signal detection by the nonisochronous Hopf oscillator, Chaos 29, 043132 (2019a).
  * Juel  _et al._ [1997] A. Juel, A. G. Darbyshire, and T. Mullin, The effect of noise on Pitchfork and Hopf bifurcations, Proc. R. Soc. A 453, 2627–2647 (1997).
  * Doan  _et al._ [2018] T. S. Doan, M. Engel, J. S. Lamb, and M. Rasmussen, Hopf bifurcation with additive noise, Nonlinearity 31, 4567–4601 (2018).
  * Risler  _et al._ [2004] T. Risler, J. Prost, and F. Jülicher, Universal critical behavior of noisy coupled oscillators, Phys. Rev. Lett. 93, 175702 (2004).
  * Faber and Bozovic [2019b] J. Faber and D. Bozovic, Chaotic dynamics enhance the sensitivity of inner ear hair cells, Sci. Rep. 9, 18394 (2019b).
  * Chow and Buice [2015] C. Chow and M. Buice, Path Integral Methods for Stochastic Differential Equations, J. Math. Neurosc. 5 (2015).
  * Barabási and Stanley [1995] A.-L. Barabási and H. E. Stanley,  _Fractal Concepts in Surface Growth_ (Cambridge University Press, 1995).
  * Binney  _et al._ [1992] J. Binney, N. Dowrick, A. Fisher, and M. Newman,  _The Theory of Critical Phenomena: An Introduction to the Renormalization Group_ (Clarendon Press, Oxford, 1992).
  * Beggs and Plenz [2003] J. M. Beggs and D. Plenz, Neuronal avalanches in neocortical circuits, J Neurosci. 23, 11167 (2003).
  * Mazzoni  _et al._ [2007] A. Mazzoni, F. D. Broccard, E. Garcia-Perez, P. Bonifazi, M. E. Ruaro, and V. Torre, On the Dynamics of the Spontaneous Activity in Neuronal Networks, PLoS ONE 2, e439 (2007).
  * Benayoun  _et al._ [2010] M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, Avalanches in a Stochastic Model of Spiking Neurons, PLoS Comput Biol 6, 1000846 (2010).
  * Bak  _et al._ [1987] P. Bak, C. Tang, and K. Wiesenfeld, Self-organized criticality: An explanation of the 1/f noise, Phys. Rev. Lett. 59, 381 (1987).
  * Coppersmith and Littlewood [1987] S. N. Coppersmith and P. B. Littlewood, Pulse-duration memory effect and deformable charge-density waves, Phys. Rev. B 36, 311 (1987).
  * Tang  _et al._ [1987] C. Tang, K. Wiesenfeld, P. Bak, S. Coppersmith, and P. Littlewood, Phase organization, Phys. Rev. Lett. 58, 1161 (1987).
  * Rodríguez Martínez and Castillo Parra [2018] G. A. Rodríguez Martínez and H. Castillo Parra, Bistable perception: Neural bases and usefulness in psychological research, Int. J. Psychol. Res. 11, 63–76 (2018).
  * Rankin  _et al._ [2015] J. Rankin, E. Sussman, and J. Rinzel, Neuromechanistic model of auditory bistability, PLoS Comput. Biol. 11, 10.1371/journal.pcbi.1004555 (2015).
  * Curtu  _et al._ [2019] R. Curtu, X. Wang, B. W. Brunton, and K. V. Nourski, Neural signatures of auditory perceptual bistability revealed by large-scale human intracranial recordings, J. Neurosc. 39, 6482–6497 (2019).
  * Cao  _et al._ [2016] R. Cao, A. Pastukhov, M. Mattia, and J. Braun, Collective activity of many bistable assemblies reproduces characteristic dynamics of multistable perception, J. Neurosc. 36, 6957–6972 (2016).
  * Stiefel and Ermentrout [2016] K. M. Stiefel and G. B. Ermentrout, Neurons as oscillators, J. Neurophysiol. 116, 2950–2960 (2016).
  * Liu  _et al._ [2012] C.-H. Liu, O. J. Tzeng, D. L. Hung, P. Tseng, and C.-H. Juan, Investigation of bistable perception with the “Silhouette spinner”: Sit still, spin the dancer with your will, Vision Res. 60, 34–39 (2012).
  * Collera [2022] J. A. Collera, Bubbling, bistable limit cycles and quasi-periodic oscillations in queues with delayed information, Symmetry 14, 376 (2022).