# Acceleration effect caused by the Onsager reaction term in a frustrated coupled oscillator system

## Abstract

The role of the Onsager reaction term (ORT) is not yet well understood in frustrated coupled oscillator systems, since the Thouless-Anderson-Palmer (TAP) and replica methods cannot be directly applied to these non-equilibrium systems. In this paper, we consider two oscillator associative memory models, one with symmetric and one with asymmetric dilution of coupling. These two systems are ideal for evaluating the effect of the ORT, because, with the exception of the ORT, they have the same order parameter equations. We found that the two systems have identical macroscopic properties, except for the acceleration effect caused by the ORT. This acceleration effect does not exist in any equilibrium system.

Coupled oscillators are of intrinsic interest in many branches of physics, chemistry and biology. Simple coupled-oscillator models involving uniform and global coupling have been investigated in some detail, and it has been found that they can be used to model many types of chemical reactions in solution kuramoto0 . However, in the modeling of more complicated phenomena, including those studied in the fields of neuronal systems, it is more necessary to consider coupled oscillators with frustrated couplings.

The Onsager reaction term (ORT), which describes the effective self-interaction, is of great importance in obtaining a physical understanding of frustrated random systems, because the presence of such an effective self-interaction is one of the characteristics that distinguish frustrated and non-frustrated systems. In the case of equilibrium systems, we can rigorously evaluate the effect of the ORT in the Thouless-Anderson-Palmer (TAP) framework Mezard , and/or using the replica method fukai2 . However, we cannot directly apply these systematic methods to non-equilibrium coupled-oscillator systems. We can define a formal Hamiltonian function on such systems. However, Perez et al. proved that the ground states of this Hamiltonian are not stationary states of the dynamics Perez . Therefore, it is impossible to construct a theory based on a free energy for this class of systems. For this reason, in order to evaluate the macroscopic quantities in such systems that include an ORT, the self-consistent signal to noise analysis (SCSNA) which can be applied to systems without the Hamiltonian function has been used fukai3 . The mathematical treatment of this method is similar to that of the cavity method Mezard . Some results obtained from the SCSNA consist with results from the replica method, but this method includes a few heuristic steps. While the SCSNA has indicated some interesting results, it is not sufficient to give a complete understanding of frustrated systems, and for this reason many theoretically fundamental questions remain in the study of such systems. In fact, even the existence of the type of self-interaction that can be described by the ORT is the subject of some debate aonishi2 ; yoshioka .

In this paper, we discuss an effect of the ORT that exists only in frustrated globally coupled oscillator systems, and in particular cannot be found in equilibrium systems. In order to make this effect clear, it would be ideal for us to compare two frustrated systems that, with the exception of the different quantity of the ORT, have the same order parameter equations. In addition, it is desirable for these systems to have a clear correspondence with an equilibrium system, because the effects of the ORT are well understood in equilibrium systems.

In consideration of the above-mentioned points, a system of the form:

$`{\displaystyle \frac{d\varphi _i}{dt}}=\omega _i+{\displaystyle \underset{ji}{\overset{N}{}}}J_{ij}\mathrm{sin}(\varphi _j\varphi _i+\beta _{ij}+\beta _0),`$ (1)   
is ideal. In fact, such systems are well known as models of coupled oscillator systems kuramoto0 ; Park . Here, $`\varphi _i`$ is the phase of the $`i`$ th oscillator (with a total of $`N`$) and $`\omega _i`$ represents its quenched natural frequency. The natural frequencies are randomly distributed with a density represented by $`g(\omega )`$. We restrict $`g(\omega )`$ to a unimodal symmetric distribution for satisfying the condition that there exists one large cluster of synchronous oscillators. Also in Eq. (1), $`J_{ij}`$ and $`\beta _{ij}`$ denote the amplitude of coupling from unit $`j`$ to unit $`i`$ and its delay, respectively. In the present study, we have selected the following two generalized Hebb learning rules with random dilutions aoyagi2  to determine $`J_{ij}`$ and $`\beta _{ij}`$:

$`J_{ij}\mathrm{exp}(i\beta _{ij})={\displaystyle \frac{c_{ij}}{cN}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}\xi _i^\mu \overline{\xi }_j^\mu ,\xi _i^\mu =\mathrm{exp}(i\theta _i^\mu ),`$ (2)   
$`c_{ij}=\{\begin{array}{cc}1\hfill & \mathrm{with}\mathrm{probability}c\hfill \\ 0\hfill & \mathrm{with}\mathrm{probability}1c\hfill \end{array},`$ (5)   
where $`\overline{}`$ means the complex conjugate. $`\{\theta _i^\mu \}_{i=1,\mathrm{},N,\mu =1,\mathrm{},p}`$ are the phase patterns to be stored in the present model and are assigned to random numbers with a uniform probability on the interval $`[0,2\pi )`$. $`\mu `$ is an index of stored patterns and $`p`$ is the total number of stored patterns. We define a parameter $`\alpha `$ (the loading rate) by $`\alpha =p/N`$. When $`\alpha O(1)`$, the system has frustration. The quantity $`c_{ij}`$ is the dilution coefficient. Let $`c_{ij}=1`$ if there is a non-zero coupling from unit $`j`$ to unit $`i`$ and $`c_{ij}=0`$ otherwise. The number of fan-in (fan-out) is restricted to $`O(N)`$, i.e., $`cO(1)`$. Here, we consider both the cases of symmetric dilution (i.e., $`c_{ij}=c_{ji}`$) and asymmetric dilution (i.e., $`c_{ij}`$ and $`c_{ji}`$ are independent random variables)okada . The quantity $`\beta _0`$ in Eq. (1) represents a uniform bias. Due to the effect of this bias, the mutual interaction between a pair of oscillators is asymmetric, even if $`J_{ij}=J_{ji}`$ and $`\beta _{ij}=\beta _{ji}`$. Such an unbalanced mutual interaction is the essence of the acceleration (deceleration) effect meunier ; Hansel1 ; kurata . In the case of $`g(\omega )=\delta (\omega \omega _0)`$, $`\beta _0=0`$ and $`c_{ij}=c_{ji}`$, this system can be mapped to a XY-spin system aoyagi2 ; cook ; okuda . In this way, we can make a bridge between the frustrated coupled oscillator system and the equilibrium system.

Let us consider steady states of the system in the limit $`t\mathrm{}`$. Our theory is based on the condition that there exists one large cluster of oscillators synchronously locked at frequency $`\mathrm{\Omega }`$ and the number of this cluster scales as $``$ $`O(N)`$. Under such a condition, Daido demonstrated through a scaling plot obtained from numerical simulations that variation of order parameter scales as $`O(1/\sqrt{N})`$ in ferromagnetic systems with one large synchronous cluster daido2 . According to this result, we assume that the self-averaging property holds in our system and order parameters are constant in the limit $`N\mathrm{}`$. These assumptions in our theory were also introduced by Sakaguchi and Kuramoto (SK) kuramoto .

Redefining $`\varphi _i`$ according to $`\varphi _i\varphi _i+\mathrm{\Omega }t`$ and substituting this into Eq. (1), we obtain

$`{\displaystyle \frac{d\varphi _i}{dt}}+\omega _i\mathrm{\Omega }=\mathrm{sin}(\varphi _i)h_i^R\mathrm{cos}(\varphi _i)h_i^I,`$ (6)   
where $`h_i`$ represents the so-called “local field”, which is described as

$`h_i=h_i^R+ih_i^I=`$   
$`e^{i\beta _0}\left({\displaystyle \underset{\mu }{\overset{p}{}}}\xi _i^\mu m^\mu +{\displaystyle \frac{1}{N}}{\displaystyle \underset{\mu }{\overset{p}{}}}{\displaystyle \underset{ji}{\overset{N}{}}}{\displaystyle \frac{c_{ij}c}{c}}\xi _i^\mu \overline{\xi }_j^\mu s_j\alpha s_i\right).`$ (7)   
For convenience, we write $`s_i=\mathrm{exp}(i\varphi _i)`$. The order parameter $`m^\mu `$, which is the overlap between the system state $`\{s_i\}_{i=1,\mathrm{},N}`$ and embedded pattern $`\{\xi _i^\mu \}_{i=1,\mathrm{},N}`$, is defined as

$`m^\mu ={\displaystyle \frac{1}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}\overline{\xi }_j^\mu s_j.`$ (8)   
The effect of the second term of Eq. (7), i.e., $`\frac{1}{N}_\mu ^p_{ji}^N\frac{c_{ij}c}{c}\xi _i^\mu \overline{\xi }_j^\mu s_j`$, is equivalent to that of an effect of additive coupling noise aoyagi2 ; Sompolinsky2 ; okada . In the limit $`c0`$, with $`\alpha /c`$ kept finite, our system is reduced to a glass oscillator. Therefore, our theory proposed here can cover two types of frustrated systems, the oscillator associative memory and the glass oscillator.

In general, the fields $`h_i^R`$ and $`h_i^I`$ involve the ORT corresponding to the effective self-feedback fukai3 . We must eliminate the reaction term from these fields. Here, we assume that the local field splits into a “pure” effective local field, $`\stackrel{~}{h}_i=\stackrel{~}{h}_i^R+i\stackrel{~}{h}_i^I`$, and the ORT, $`\mathrm{\Gamma }s_i`$:

$`h_i=\stackrel{~}{h}_i+\mathrm{\Gamma }s_i.`$ (9)   
Here we have neglected the complex conjugate term of the ORT leading to a higher-harmonic term of the response function aonishi . This can be done in the present model because we employ generalized Hebb learning rules. Hence, by substituting Eq. (9) into Eq. (6), we obtain the equation

$`{\displaystyle \frac{d\varphi _i}{dt}}+\omega _i\stackrel{~}{\mathrm{\Omega }}=\mathrm{sin}(\varphi _i)\stackrel{~}{h}_i^R\mathrm{cos}(\varphi _i)\stackrel{~}{h}_i^I,`$ (10)   
$`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }\left|\mathrm{\Gamma }\right|\mathrm{sin}(\psi ),\psi =\mathrm{Arg}\left(\mathrm{\Gamma }\right),`$ (11)   
which does not contain the reaction term. The quantity $`\stackrel{~}{\mathrm{\Omega }}`$ here represents the effective frequency of synchronous oscillators. We can regard $`\stackrel{~}{\mathrm{\Omega }}`$ as the renormalized version of $`\mathrm{\Omega }`$, of which the ORT has been pulled out, and therefore $`\stackrel{~}{\mathrm{\Omega }}`$ takes a different value from the observable $`\mathrm{\Omega }`$ in general. Thus, $`\mathrm{\Omega }\stackrel{~}{\mathrm{\Omega }}`$ represents the contribution of the ORT to the acceleration (deceleration) effect. $`\stackrel{~}{\mathrm{\Omega }}`$ is one of the order parameters of our theory. In the analysis that follows, $`\stackrel{~}{h}_i`$ and $`\mathrm{\Gamma }`$ are obtained in a self-consistent manner. Under the above assumption expressed by Eq. (9), by applying SK theory kuramoto  to Eq. (10), we obtain the average of $`s_i`$ over $`\omega _i`$:

$`s_i_{\omega _i}`$ $`=`$ $`\stackrel{~}{h_i}{\displaystyle _{\pi /2}^{\pi /2}}𝑑\varphi g\left(\stackrel{~}{\mathrm{\Omega }}+|\stackrel{~}{h}|\mathrm{sin}\varphi \right)\mathrm{cos}\varphi \mathrm{exp}(i\varphi )`$ (12)   
$`+`$ $`i\stackrel{~}{h}_i{\displaystyle _0^{\pi /2}}𝑑\varphi {\displaystyle \frac{\mathrm{cos}\varphi (1\mathrm{cos}\varphi )}{\mathrm{sin}^3\varphi }}`$   
$`\times \left\{g\left(\stackrel{~}{\mathrm{\Omega }}+{\displaystyle \frac{|\stackrel{~}{h}_i|}{\mathrm{sin}\varphi }}\right)g\left(\stackrel{~}{\mathrm{\Omega }}{\displaystyle \frac{|\stackrel{~}{h}_i|}{\mathrm{sin}\varphi }}\right)\right\}.`$   
Equation (6) implies that $`s_i`$ is a function of $`h_i`$, $`\omega _i\mathrm{\Omega }`$ and $`t`$, and to make this explicit we write

$$s_i=X(h_i,\omega _i\mathrm{\Omega },t).$$
(13)   
Note that $`s_i`$ is not a function of renormalized $`\stackrel{~}{h}_i`$ and $`\stackrel{~}{\mathrm{\Omega }}`$ but, rather, the bare $`h_i`$ and $`\mathrm{\Omega }`$ in appearing Eq. (6). We can properly evaluate the ORT with this careful treatment. Here, we assume that microscopic memory effect can be neglected in the $`t\mathrm{}`$ limit. In this analysis, we focus on memory retrieval states in which the configuration has appreciable overlap with the condensed pattern $`𝝃^1`$ ($`m^1O(1)`$) and has tiny overlap with the uncondensed patterns $`𝝃^\mu `$ for $`\mu >1`$ ($`m^\mu O(1/\sqrt{N})`$). Under this assumption, we estimate the contribution of the uncondensed patterns using the SCSNAfukai3 , and determine $`\stackrel{~}{h}_i`$ in a self-consistent manner. Finally, the equations relating the order parameters $`|m^1|`$, $`U`$ and $`\stackrel{~}{\mathrm{\Omega }}`$ are obtained using the self-consistent local field:

$`|m^1|e^{i\beta _0}`$ $`=`$ $`\stackrel{~}{X}(x_1,x_2;\stackrel{~}{\mathrm{\Omega }})_{x_1,x_2},`$ (14)   
$`Ue^{i\beta _0}`$ $`=`$ $`F_1(x_1,x_2;\stackrel{~}{\mathrm{\Omega }})_{x_1,x_2},`$ (15)   
where $`\mathrm{}_{x_1,x_2}`$ is the Gaussian average over $`x_1`$ and $`x_2`$, $`\mathrm{}_{x_1,x_2}=Dx_1Dx_2\mathrm{}`$. The quantity $`U`$ corresponds to the susceptibility, which is the measure of the sensitivity to external fields. Since the present system possesses rotational symmetry with respect to the phase $`\varphi _i`$, we can safely set the condensed pattern as $`\xi _i^1=1(i=1,\mathrm{},N)`$. Now, $`\stackrel{~}{h}`$, $`\stackrel{~}{X}`$, $`F_1`$ and $`Dx_1Dx_2`$ can be expressed as follows:

$`Dx_1Dx_2={\displaystyle \frac{dx_1dx_2}{2\pi \rho ^2}}\mathrm{exp}\left({\displaystyle \frac{x_1^2+x_2^2}{2\rho ^2}}\right),`$ (16)   
$`\rho ^2={\displaystyle \frac{\alpha }{2}}\left({\displaystyle \frac{1}{|1U|^2}}+{\displaystyle \frac{1c}{c}}\right),`$ (17)   
$`\stackrel{~}{h}=|m^1|+x_1+ix_2,`$ (18)   
$`\stackrel{~}{X}(x_1,x_2;\stackrel{~}{\mathrm{\Omega }})=s_\omega ,`$ (19)   
$`F_1(x_1,x_2;\stackrel{~}{\mathrm{\Omega }})={\displaystyle \frac{s_\omega }{\stackrel{~}{h}}},`$ (20)   
where $`s_\omega `$ in Eqs. (19) and (20) is written as Eq. (12).

In the case of the symmetric diluted system, $`\mathrm{\Gamma }`$ can be expressed as

$`\mathrm{\Gamma }e^{i\beta _0}={\displaystyle \frac{\alpha U}{1U}}+{\displaystyle \frac{\alpha (1c)}{c}}U.`$ (21)   
In the case of the asymmetric diluted system, on the other hand, we have

$`\mathrm{\Gamma }e^{i\beta _0}={\displaystyle \frac{\alpha U}{1U}}.`$ (22)   
$`\stackrel{~}{h}`$ and $`\stackrel{~}{\mathrm{\Omega }}`$ are the renormalized versions of $`h`$ and $`\mathrm{\Omega }`$ respectively, of which the ORT has been pulled out, and thus $`\stackrel{~}{h}`$ and $`\stackrel{~}{\mathrm{\Omega }}`$ are independent of the ORT. Therefore, the two models we consider have identical order parameter equations, (14) and (15), written in the term of the renormalized quantities $`\stackrel{~}{h}`$ and $`\stackrel{~}{\mathrm{\Omega }}`$. According to Eq. (11), the difference of the ORT in Eqs. (21) and (22) only leads a different value of the observable $`\mathrm{\Omega }`$ only when $`\beta _00`$. In this way we are able to clearly separate the effect of the ORT, and therefore, by observing the macroscopic parameter $`\mathrm{\Omega }`$ of these two systems, we can analyze the effect of the ORT qualitatively and quantitatively.

The distribution of resultant frequencies $`\overline{\omega }`$ in the memory retrieval state, which is denoted $`p(\overline{\omega })`$, becomes

$`p(\overline{\omega })=r\delta (\overline{\omega }\mathrm{\Omega })`$   
$`+{\displaystyle Dx_1Dx_2\frac{g\left(\stackrel{~}{\mathrm{\Omega }}+(\overline{\omega }\mathrm{\Omega })\sqrt{1+\frac{|\stackrel{~}{h}|^2}{(\overline{\omega }\mathrm{\Omega })^2}}\right)}{\sqrt{1+\frac{|\stackrel{~}{h}|^2}{(\overline{\omega }\mathrm{\Omega })^2}}}},`$ (23)   
$`r={\displaystyle Dx_1Dx_2|\stackrel{~}{h}|_{\pi /2}^{\pi /2}𝑑\varphi g\left(\stackrel{~}{\mathrm{\Omega }}+|\stackrel{~}{h}|\mathrm{sin}\varphi \right)\mathrm{cos}\varphi }.`$ (24)   
The $`\delta `$-function in Eq. (23) indicates the cluster of oscillators synchronously locked at frequency $`\mathrm{\Omega }`$. The value $`r`$ is the ratio between the number of synchronous oscillators and the total number of oscillators $`N`$. The second term in Eq. (23) represent a distribution of asynchronous oscillators.

If $`\beta _0=0`$ and $`g(\omega )`$ is symmetric, our theory reduces to the previously proposed theory aonishi2 . In the case $`g(\omega )=\delta (\omega )`$, $`\beta _0=0`$ and $`c_{ij}=c_{ji}`$, where the present model reduces to XY-spin systems (frustrated equilibrium systems), our theory coincides with the replica theory aoyagi2 ; cook  and the SCSNAokuda . In addition, in the limit of $`\alpha 0`$, where the present model reduces to uniform non-equilibrium systems, our theory reproduces the SK theory kuramoto .

In the numerical simulations we now discuss, we chose the distribution of natural frequencies as $`g(\omega )=(2\pi \sigma ^2)^{1/2}\mathrm{exp}(\omega ^2/2\sigma ^2)`$. Also, we set $`\sigma =0.2`$, $`\beta _0=\pi /20`$, and $`c=0.5`$. Figure 1(a) displays $`\mathrm{\Omega }`$ as a function of $`\alpha `$ in the memory retrieval states, where the solid curves were obtained theoretically, and the data points with error bars represent results obtained by numerical simulation. As seen from Figure 1(a), the oscillator’s rotation for the symmetric diluted system is faster than that for the asymmetric diluted system. As previously discussed, $`\stackrel{~}{\mathrm{\Omega }}`$ in Figure 1(a), which represents the effective frequency of synchronous oscillators, does not depend on the type of dilution, while the observed $`\mathrm{\Omega }`$ depends strongly on it. This dependence is due to the existence of the ORT $`\mathrm{\Gamma }`$. In the case that the local field $`h`$ does not contain the ORT yoshioka , the plots obtained from the numerical simulations for both models should fit the curve of $`\stackrel{~}{\mathrm{\Omega }}`$. Therefore, the dependence of the observed $`\mathrm{\Omega }`$ on the type of dilution is strong evidence for the existence of the ORT in the present system. In this figure, we have shifted the numerical values of $`\mathrm{\Omega }`$ at $`\alpha =0`$ (in the computer simulation) to their corresponding theoretical values at $`\alpha =0`$, in order to cancel fluctuations of the mean value of $`g(\omega )`$ caused by the finite size effect.

Figures 1(b) and (c) display the distributions of the resultant frequencies for the symmetric dilution and asymmetric dilution systems, respectively. As these figures reveal, the theory (solid curve) is in good agreement with the simulation results (histogram). According to the results given in Figures 1(b) and (c), the distributions of the resultant frequencies for the symmetric diluted system is identical to that for the asymmetric diluted system, except for a slight difference between those positions caused by the ORT. From this result we can conclude that the mean field $`\stackrel{~}{h}`$ of the symmetric diluted system is identical to that of the asymmetric diluted system, since $`\stackrel{~}{h}`$ reflects the distribution of resultant frequencies, as represented by Eq. (23). Figures 2(a) and (b) display $`|m^1|`$ as a function of $`\alpha `$ for symmetric and asymmetric diluted systems, respectively. Here, the solid curves were obtained theoretically, and the data plots represent the results obtained from the numerical simulations. According to the results given in Figures 2(a) and (b), the critical memory capacity of the asymmetric diluted systems obtained from numerical simulation is slightly smaller than that of the asymmetric diluted systems, because asymmetric dilutions break the detailed balance, and then, asymmetric dilutions weaken a memory state. At this stage, there is no theory to rigorously treat a system with asymmetric interaction. Almost theoretical studies of asymmetric systems are based on the naive assumption that there exist steady states like symmetric systems okada .

In conclusion, we have found that the symmetric and asymmetric diluted systems have the same macroscopic properties, with the exception of the the acceleration (deceleration) effect caused by the ORT. The quantity of the ORT depends on the type of dilution, and this dependence leads to a difference in the rotation speed of the oscillators for the two cases, as shown in Fig. 1(a). The acceleration (deceleration) effect caused by the ORT is a phenomenon peculiar to non-equilibrium systems, since this effect only exists when $`\beta _00`$. There has been fundamental disagreement regarding the existence of the ORT in a typical system corresponding to our model with $`\beta _0=0`$ and a symmetric $`g(\omega )`$ aonishi2 ; yoshioka . In this work we have reached the following conclusion in this regard. Even if $`\beta _0=0`$ and $`g(\omega )`$ is symmetric, the ORT exists in the bare local field given by Eq. (7). In this case, the effect of the ORT is invisible, since it cancels out of Eq. (10).