Brownian Behavior in Coupled Chaotic Oscillators

: Since the dynamical behavior of chaotic and stochastic systems is very similar, it is some-times difﬁcult to determine the nature of the movement. One of the best-studied stochastic processes is Brownian motion, a random walk that accurately describes many phenomena that occur in nature, including quantum mechanics. In this paper, we propose an approach that allows us to analyze chaotic dynamics using the Langevin equation describing dynamics of the phase difference between identical coupled chaotic oscillators. The time evolution of this phase difference can be explained by the biased Brownian motion, which is accepted in quantum mechanics for modeling thermal phenomena. Using a deterministic model based on chaotic Rössler oscillators, we are able to reproduce a similar time evolution for the phase difference. We show how the phenomenon of intermittent phase synchronization can be explained in terms of both stochastic and deterministic models. In addition, the existence of phase multistability in the phase synchronization regime is demonstrated.


Introduction
Brownian motion is the random motion of particles in a fluid medium as a result of collisions with fluid molecules. This transport phenomenon was named after the Scottish biologist Robert Brown, who in 1827 observed through a microscope how particles are retained in cavities inside a pollen grain in water. Although this discovery of Brownian motion was remarkable, the mechanisms causing this motion had not been identified. Later, Einstein and Smolowski, in their famous papers of 1905 and 1906, respectively [1,2], linked molecular kinetic theory to Brownian motion. Their theory was experimentally verified by French physicist Jean Baptiste Perrin, who established that Brownian motion in liquids is caused by the movement of molecules, and thus convincingly proved the existence of molecules and atoms.
Brownian motion is among the simplest random processes and is akin to two other simpler and more complex stochastic processes: random walk and Donsker's theorem [3]. This universality is closely related to the versality of the normal distribution. In both cases, the use of normal distribution is often associated with mathematical convenience rather than model accuracy. Despite the great progress in understanding the nature of Brownian motion, many aspects remain unclear, and questions remain unanswered, for example, how this process depends on the interaction of particles with the environment, whether it is stochastic or deterministic, etc.
In mathematics, Brownian motion is interpreted as a Wiener process named after Norbert Wiener. This continuous random process is often used in various differential equations to describe physical phenomena that are affected by noise. In many cases, the quantum nature of equations requires a proper understanding of this process [4]. In chaos theory, the question of whether Brownian motion is stochastic or deterministic is rather philosophical. Although the equations describing chaotic phenomena are deterministic, it is inherently difficult to determine the initial conditions or parameters that govern the chaotic dynamical system with sufficient accuracy to replicate the exact same experiment in the laboratory without considering environmental noise. In recent decades, with the advent of computers, the scientific community has discovered that chaotic systems behave like stochastic ones when we perform simulations in which we can control the accuracy of parameters and initial conditions [5]. This brings us back to the philosophical question of what behavior in nature is truly random.
It is now well known that coupled chaotic oscillators, similar to limit cycles, can exhibit coherence phenomena [6]. One such phenomenon is phase synchronization, which is usually characterized by the phase difference between oscillators. Over the past few years, in articles devoted to chaotic phase synchronization, the authors reported that the phase difference between chaotic oscillator fluctuations in time, is like Brownian behavior [7,8].
In this paper, we explore the possibility of simulating the synchronous dynamics of coupled chaotic oscillators with stochastic differential equations, although it is not at all obvious that random behavior can be reproduced using deterministic equations. To do this, we consider a stochastic model and describe the biased Brownian motion in the same way as it was done in the framework of quantum mechanics. We show that the parameters governing such a system can be derived by studying the random motion trajectories. We then demonstrate that the phase difference between two coupled chaotic oscillators behaves like a trajectory derived from the stochastic model, and that parameters obtained with this model can be used to describe the chaotic behavior of the phase difference.

Stochastic Model
The simplest chaotic dynamical systems are usually described by ordinary differential equations that model nonlinear oscillators with noticeable fluctuations in amplitude and frequency. Systems of this type often have periodic solutions under the influence of external noise. Hence, from a nonchaotic deterministic model, we can come to a stochastic model that exhibits generalized dynamics of simple chaotic oscillators. This is the case we discuss below.

Modified Kuramoto Model
One of the most common deterministic models which represents limit cycle synchronization is the Kuramoto model [9], defined by the phase time derivative of each limit cycle as where ω i is the natural frequency of i-th limit cycle. Note that Equation (1) describes phase dynamics of a symmetric all-to-all network of N limit cycle oscillators coupled with strength K > 0. In an asymmetric network, K ij is under the sum in Equation (1) and represents the coupling strength between each pair of limit cycles. If a pair of oscillators are not connected, then K ij = 0. An important property of a chaotic oscillator is the unpredictability of amplitude fluctuations. To describe these fluctuations in the Kuramoto model, we introduce a driving term f (A i ), which takes into account the frequency dependence of the amplitude. We can now derive the modified Kuramoto model for an asymmetric network with a frequencydependent amplitude, as follows

Langevin Dynamics of Two Coupled Chaotic Oscillators
Synchronization of coupled oscillators is often characterized by their phase difference θ = φ 2 − φ 1 . Consider now a pair of unidirectionally coupled chaotic oscillators such that K 12 = 0 and K 21 = κ. Then, two differential equations which govern their dynamics can be written as follows where Equations (2a) and (2b) represent respectively master and slave oscillators. Hence, the phase difference is given by where In chaotic oscillators, the mean frequency of each isolated oscillator differs from the natural frequency due to nonlinearity. The mismatch between these two frequencies can be measured from the power spectrum of each oscillator, determining a maximum known as the dominant frequency (this effect will be discussed in the next section). We can now rewrite Equation (3) to obtain the following Langevin equation where ∆Ω = Ω 2 − Ω 1 is the average mismatch between dominant frequencies of two chaotic oscillators [10]. The averaging is performed over the power spectra for different initial conditions, since we expect diffusion of the phase difference, and therefore it is nonergodic. It is easy to see that when the oscillators are identical and uncoupled, the average mismatch vanishes.

Mean Angular Velocity and Diffusion Coefficient
The time average of Equation (4) is calculated as Since standard white noise has a null expected value over time, in the ensemble sense the average of Equation (5) holds where · means the ensemble average. Then, ∆Ω is the mean angular velocity of the Langevin Equation (4). Taking into account Equation (6), the variance of θ(t) can be defined as where · is also the ensemble average. For a normal diffusion process [11], when t → ∞ the variance is expected to be linear matching σ 2 = 2Dt α with α = 1, thus

Deterministic Model
Consider two identical Rössler oscillators with diffusive unidirectional coupling. The system dynamics is modelled by the following equationṡ where a = 0.165, b = 0.2, c = 10, and ω 1 = ω 2 = 1 are fixed parameters resulting in a chaotic regime. Equations (8a) and (8b) describe respectively the chaotic behavior of the master and slave oscillators. Since the oscillators are identical, they have the same natural frequency. Thus, when the coupling strength parameter κ vanishes, they oscillate with the same dominant frequency. Due to nonlinearity, the natural and dominant frequencies do not match. If a small mismatch exists between natural frequencies, the phenomenon of phase synchronization can be observed [7,[12][13][14][15]. When the coupling strength is very small, the phase difference fluctuates as a random walk [8].

Chaotic Phase Difference
Since the Rössler attractor has a spiral shape, the phase of each oscillator can be calculated using the time series [16]. In projection on the (x, y) plane (see Figure 1), the position vector forms with the positive part of the x axis the following relative phase where arctan(·) ∈ [−π, π) and (x , y ) are the Cartesian axes parallel to the (x, y) axes, in which the origin (x * , y * ) is the attractor center. The relative phase is bounded and exhibits a jump every time a revolution is completed. To calculate the phase difference, we must first compute the absolute phase of each oscillator. To do this, we unwrap radian phases by replacing absolute jumps greater than π to their 2π complement. Therefore, the absolute phase is and t k occurs when y(t k ) = y * and x(t k ) < x * with k ∈ N 0 . Since the Rössler oscillator rotates counterclockwise, we consider φ(t) as a continuous increasing function. Linearizing the entire system in Equation (8) allows us to find the coordinates of the center of each attractor. Since the locations of two centers coincide and are independent of the coupling strength κ, we can now calculate the phase difference as θ(t) = φ 2 (t) − φ 1 (t).

Brownian-like Phase Diffusion
With the intention of studying the phase difference diffusion, we carry out numerical simulations using the fifth-order Runge-Kutta method to obtain the time series of Equation (8) for various random initial conditions. We can obtain phase synchronization when the coupling strength exceeds a certain threshold value (κ > κ th ). In the phase synchronization regime, the phase difference fluctuates around a certain fixed value. Whereas for the coupling strength below the threshold, the phase difference changes linearly with time or switches to 2π radians.
The phase difference fluctuation looks like a random process similar to Brownian motion. In order to characterize this behavior, we compute probability density functions (PDF) [17] and mean power spectral densities (PSD) for the set of trajectories initiated from random initial conditions. Now, we consider the phase difference evolution for different values of the coupling strength κ.
When the oscillators are uncoupled (κ = 0), there is no interaction between them. As seen from Figure 2, both the phase difference and the corresponding PDF increase or decrease in time depending on the initial conditions. For weakly coupled oscillators, when the coupling strength is below the phase synchronization threshold (κ < κ th ), all phase difference trajectories gradually increase in time regardless of initial conditions, as seen from Figure 3a. At the same time, the PDF reveals intermittent phase locking, represented by high probability bands, multiples of 2π (Figure 3b). This phenomenon is detected for all values of the coupling strength below the phase synchronization threshold.  Brownian motion is characterized by a β ≈ −2 power low of the mean PSD [18]. To check how close the observed phase difference is to Brownian motion PSD, we calculate log S(ω) (S(ω) being the PSD of θ(t) for each random initial condition) versus log ω. The results displayed in Figure 4 confirm the Brownian-like behavior of our deterministic system.

Matching Both Models
The phase drift shown in Figures 2 and 3 can be computed with Equation (6). Specifically, the drift in the deterministic model is the mean angular velocity ∆Ω which can be calculated as a slope of the time evolution of the average mismatch between dominant frequencies θ(t) − θ(0) (see Equation (6) and Figure 5). To find the diffusion coefficient D, we plot in Figure 6 the variance calculated with Equation (7) and fit this time dependence by a straight line. The slope of this dependence in a log scale is always near 1, that means normal diffusion, and the ordinate at the origin is the diffusion coefficient.

Itô-Langevin Form
Equation (4) can be transformed into the Itô-Langevin [19] form as follows where W t is the standard Wiener process such that dW t /dt = w(t).
The physical interpretation of Equation (4) is a random walk of a particle in a onedimensional potential [20] and the deterministic force can be written as So, the stochastic differential Equation (10) is the Itô diffusion for a biased Brownian motion Θ t in the tilted periodic potential Some potentials are plotted in Figure 8 employing the values of ∆Ω(κ) calculated in Section 4.1.

Numerical Simulations
Searching analytical solutions of stochastic equations is not an easy task. Although there are some approaches to find them (see, for example, [21]), it is not the aim of the present work. In this paper, we mainly focus on numerical simulations of Equation (10).
The variation in the diffusion coefficient D(κ) shown in Figure 7 is due to the variation of the potential V(θ; κ). So, in Equation (10) we set D = D(κ = 0) because it yields the value of the diffusion coefficient when there is no interaction, that is, V(θ) = 0.
To solve Equation (10), we apply the Euler-Maruyama method as follows where γ > 0 is the step size and δW j j∈N 0 ∼ N (0, 1). In our simulations, we set Θ 0 = θ(0), that is, we choose the initial conditions as in the last section, and the step size γ is the mean time between local extrema from deterministic trajectories.
In Figures 9-11 we compare the probability density functions between simulations performed with the deterministic and stochastic models. One can see that the biased Brownian motion is a faithful representation of the deterministic model not only for the asynchronous regime, but even for coupling strength values corresponding to the phase synchronization regime.
Another interesting phenomenon which exhibits this system is the coexistence of three attractors, as seen from Figure 11, where three branches correspond to different initial conditions. This is the manifestation of phase multistability [22,23].   . Probability density function of θ(t) over time for (left) deterministic and (right) stochastic (Θ t ) models for κ = 0.01 and the same initial conditions. In the phase synchronization regime, both models exhibit confined Brownian-like motion. The three coexisting branches indicate underlying phase multistability.

Conclusions
In this paper, we have analyzed the dynamics of two unidirectionally coupled chaotic oscillators with identical natural frequencies using a quantum mechanics approach. Nonlinearity in the basic equations leads to an increase in the oscillation frequency of the slave oscillator if the oscillators are not synchronized. This, in turn, increases the system randomness, as seen in Figure 7. In addition, the diffusion coefficient increases with increasing mean angular frequency. This phenomenon is contrary to common sense, since an increase in coherence between oscillators is usually expected when the oscillators are coupled. We suspect that when the mean angular velocity is higher, there is some kind of resonance between the harmonics of both oscillators at low frequencies, so that the trajectory of one of the oscillators will act as a "chaotic stationary wave" with respect to the trajectory of the other oscillator. Thus, the phase difference would be greater or less than at lower values of the mean angular velocity.
In addition, as shown in Figures 9 and 10, intermittent phase synchronization occurs in the asynchronous regime. This phenomenon was revealed when studying the probability density of the phase difference. The stochastic model showed that some kind of force acts in this direction, which becomes stronger as the coupling strength is increased. Thus, the biased Brownian motion in an inclined potential, observed in quantum mechanics, reproduces the dynamics of the chaotic deterministic system.
Finally, we observed phase multistability in the weakly phase synchronized chaotic system, which manifested itself in the coexistence of attractors with different phase differences, shifted by 2π.