Extracting order parameters from global measurements with application to coupled electrochemical oscillators

Physica D 205 (2005) 57–69

Extracting order parameters from global measurements with application to coupled electrochemical oscillators夽 Yumei Zhai a , Istv´an Z. Kiss a , Hiroaki Daido b , John L. Hudson a,∗ b

a Department of Chemical Engineering, 102 Engineers’ Way, University of Virginia, Charlottesville, VA 22904-4741, USA Department of Mathematical Sciences, Graduate School of Engineering, University of Osaka Prefecture, Sakai 599-8531, Japan

Available online 25 March 2005

Abstract The calculation of order parameters to characterize the degree of synchronization usually requires the measurement of time series of observable quantities for all individual dynamical elements. In this article, a new method of extracting the Kuramoto order from global measurements for weakly coupled, nearly identical oscillators is established and successfully tested with electrochemical populations of smooth, relaxational and chaotic oscillators. An experimental study on a larger array of coupled smooth oscillators was performed and the results confirmed the theoretical predictions on size effects including enhanced fluctuations around the synchronization transition. The proposed method opens up the possibility for the study of synchronization in large populations of interacting oscillators in many fields including biology where measurements are limited to global signals and few local sites. © 2004 Elsevier B.V. All rights reserved. PACS: 82.40.Np; 05.45.Xt; 87.10.+e Keywords: Order parameter; Synchronization; Populations of oscillators; Metal electrodissolution

1. Introduction Large populations of interacting oscillatory entities are commonly observed in physical, chemical, and biological systems [1,2]. Mutual synchronization among

DOI of original article:10.1016/j.physd.2004.09.014.

夽 This paper was originally and erroneously published in

Physica D, vol. 199 (3,4), pp. 387–399. ∗ Corresponding author. Tel.: +1 434 924 6275; fax: +1 434 982 2658. E-mail address: [email protected] (J.L. Hudson). 0167-2789/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2004.09.017

elements in an oscillator population has been shown to be important in interpreting glycolitic oscillations [3,4], alpha rhythms in the brain [5], and aggregate beating of heart cells [1,6]. Engineering applications include microwave communications [7], high-power laser devices [8], and superconducting electronic systems [9]. Such a system can be modeled by a large assembly of coupled limit-cycle oscillators with intrinsic frequencies different from one another. Kuramoto, based on Winfree’s intuition about phase models [10], made a major advance in the theory of the onset of synchronization in populations with weak global coupling

58

Y. Zhai et al. / Physica D 205 (2005) 57–69

[11]. In the Kuramoto model, each oscillator of an infinite set is described by a single variable, the phase, φ, and is coupled to all other elements with equal strength K:

one another and coupling among them is sufficiently weak, we can express each observable quantity as

N dφj K
= ωj + sin(φl − φj ), dt N

where A(θ) is a periodic function with period 2π and θ j the phase of the jth oscillator (mod 2π) such that in the absence of coupling, it varies as θ j = Ωj t + constant, where Ωj is the angular frequency of the jth oscillator [2]. A(θ) may be Fourier-expanded as

(1)

l=1

where N is the system size, j = 1, . . ., N and ωj the frequency of the jth oscillator. The theory predicts a transition with increasing global coupling strength K at which some of the oscillators with originally different frequencies become coherent. Order √ parameters such  as Z = 1/N N j=1 exp(iφj ), (i ≡ −1), were introduced to characterize the extent of the synchronization and an analytical solution of order parameter as a function of coupling strength was obtained. Daido expanded the study to generic coupling functions and introduced generalized order parameters [12,13]. However, by definition, these order parameters require the knowledge of the phase of each oscillator. In experiments, very often measurements of all individual oscillators are impossible but global/mean quantities are readily available. In this paper, the measurement of order from global information for an interacting population of oscillators is established and tested with an electrochemical system experimentally and numerically. The new method was used to calculate orders of weakly coupled smooth, relaxational and phase coherent chaotic oscillators in laboratory experiments; the results are compared with the ones we obtained earlier from individual measurements where Kuramoto’s theory on the relationship of the order and coupling strength was verified [14]. In addition, a study on a larger array of globally coupled smooth electrochemical oscillators was performed based on the new method and effects of size on synchronization process of a finite-size population are discussed. Finally, the limitations of the new method are presented and suggestions for its proper application to laboratory systems are made.

2. Order from a global signal Consider first a population of periodic oscillators. Let the time series of an observable variable from the jth element be Ij (t). If all the oscillators are similar to

Ij (t) ∼ = A(θj ),

A(θ) =

∞


Ak eikθ ,

(2)

(3)

k=−∞

where Ak are Fourier coefficients (in general complex) which are determined by the waveforms of the uncoupled oscillators:  1 Tj Ak ∼ Ij (t) e−ikΩj t dt, (4) = Tj 0 where j is a representative oscillator and Tj ≡ 2π/Ωj . Let us now consider the behavior of the mean (javeraged) oscillation denoted below by I(t), i.e., I(t) ≡

N 1
Ij (t). N

(5)

j=1

Using Eq. (3) and noting that in the synchronized region, the phase variables can be written θj (t) = Ωe t + ψj ,

(6)

with the angular frequency of entrainment Ωe as well as slowly varying variables ψj , we obtain the following results I(t) =

∞


Ak Zk eikΩe t ,

(7)

k=−∞

where Zk ≡

N 1
ikψj e , N

(8)

j=1

are the generalized order parameters introduced in [15,12,13] and Z1 is the Kuramoto order parameter Z. The analysis of Z2 , Z3 , etc., in addition to Z1 , is crucial to the understanding of the asymptotic dynamics of oscillator populations through the order function that is defined using Zk [13]. These parameters are constants

Y. Zhai et al. / Physica D 205 (2005) 57–69

in the thermodynamic limit N → ∞, but exhibit persistent fluctuations for N < ∞ [16,17]. Since Zk are all slowly varying by definition, if Ak = 0, then we obtain Zk (t) ∼ =

1 Ak Te



t+Te/2 t−Te/2

I(u) e−ikΩe u du,

(9)

where Te ≡ 2π/Ωe . The coefficients Ak can be calculated from Eq. (4). The frequency Ωe can be experimentally measured. It should be noted that without appropriate symmetries of the frequency distribution f(Ω) and the coupling function h(θ), Ωe generically differs from the average of natural angular frequencies and depends on the control parameter [13]. In the nonsynchronized region, Ωe cannot be defined, but the average of natural frequencies can be used because the value of |Zk | as well as the slowly varying nature of |Zk | are not affected by this. In the right-hand-side of Eq. (9), the coefficients of the Fourier expansion of the mean signal I appear; if we denote them as Bk , then Bk Zk (t) ∼ , = Ak

(10)

where Bk =

1 Te



t+Te/2 t−Te/2

I(u) e−ikΩe u du.

(11)

Hence, using Eq. (9) or Eq. (10), the temporal variation of the order parameters can be reproduced exclusively from the mean signal.

59

3. Experiment and simulation 3.1. Experiments The experimental system is an array of 64 or 100 nickel electrodes in an 8 × 8 or 10 × 10 geometry, respectively. A standard three-compartment electrochemical cell consisting of the array of nickel working electrodes, a Hg/Hg2 SO4 /K2 SO4 reference electrode and a Pt mesh counter electrode was used. A schematic of the experimental setup is shown in Fig. 1a. Nickel dissolution on the surface of the electrode occurred at an applied potential (V) held constant by a potentiostat (EG&G Princeton Applied Research, model 273). The working electrodes are embedded in epoxy and reaction takes place only at the ends. Experiments were carried out in H2 SO4 solution at a temperature of 11 ◦ C. The electrolyte was stirred with a magnetic stirrer at a speed of 250 rpm. Current, proportional to the rate of reaction, was measured at a sampling rate of 100 Hz on each electrode to provide dynamical information of the elements of the population. (The sum of the individual currents is within 3% of the independently measured total current obtained from the potentiostat.) Zero resistance ammeters were used to convert each individual current into a potential that was recorded by a transputer (Adwin Pro, Keithley). Periodic or chaotic oscillations were observed, depending on conditions such as applied potential (V), acid concentration (c), and added external resistance [18,19]. The electrodes are 1 mm in diameter and equally spaced with 2 mm; it has been shown that the short-range coupling through the effect of stirring can be ignored with this spacing [20].

Fig. 1. Schematics of the experiments: (a) experimental setup; (b) equivalent circuit.

60

Y. Zhai et al. / Physica D 205 (2005) 57–69

The electrodes were connected to the potentiostat through one series (collective) resistor (Rs ) and through N parallel (individual) resistors (RP ), where N is the number of electrodes. Heterogeneity due to variations in the surface of the electrodes and transport always exists in the experiments and causes a frequency distribution in the individual oscillating currents [21]. In addition, we used small random resistors (with a standard deviation of 21 ) to increase the variance of the frequency distribution of the periodic populations; experiments without these resistors gave similar results to those shown but with a slight loss in reproducibility. In the studies of periodic oscillators, mean parallel resistors RP = 652  and c = 3 mol/L were used; V = 1.075 V (versus Hg/Hg2 SO4 /K2 SO4 electrode) for smooth oscillators and V = 1.22 V for relaxational oscillators. Population of chaotic oscillators were obtained with identical individual resistors of RP = 906 , c = 4.5 mol/L and V = 1.310 V. The collective resistor couples the electrodes globally [21,22]. To study the effect of various global coupling strengths on the behavior of the system, we held the total external resistance (Rtot = Rs + RP /N) constant while the fraction dedicated to individual currents, as opposed to the total current, was varied. Thus all other parameters of the system are constant (i.e., the kinetics and the parameters of the individual oscillators are not changed) while the strength of global coupling can be varied. A global coupling parameter K is defined as K=

Rs /Rtot . 1 − Rs /Rtot

(12)

For K = 0, the external resistance furnishes no additional global coupling; for K = ∞, maximal external global coupling is achieved. 3.2. Simulations A simple equivalent circuit of an electrode array is shown in Fig. 1b. For clarity, only two electrodes are shown in the figure. The electrical double layer at the metal–solution interface is described as a double layer capacitor (Cd ) and a nonlinear Faradaic resistor in parallel. The current passing through the Faradaic resistor (iF ) is proportional to the rate of reaction of each electrode, which is determined by the reaction kinetics. The individual resistors (rj ) and collective re-

sistor (Rs ) are inserted into the circuitry in parallel and in series, respectively. The individual current through each electrode (Ij ) can be directly measured in experiments and thus serves as observable variable. However, the direct dynamical variable of this system, i.e., that which appears in the governing ODE model below, is the potential drop through the double layer (ej ) of the electrode. These two quantities are related by Ij =

(v − ej ) , rj

(13)

where v is the potential drop cross the  individual resistor and the double layer (v = V − Rs Ij ). Simulations on coupled electrochemical oscillators with population size of 64 or 512 were carried out for comparison with the experimental results. We used a model of anodic electrodissolution of a single nickel electrode proposed by Haim et al. [23]. The model in a dimensionless form involves two variables: the dimensionless double layer potential drop (e) and the surface coverage of NiO + NiOH (θ). Based on this kinetic model of a single electrode, the following dimensionless equations were derived for coupled nickel electrodes [24]: dej V − ej 1 = − iF,j (θj , ej ) + K(emean − ej ), dt Rj R0 (14) Γj

dθj exp(0.5ej ) = (1 − θj ) dt 1 + Ch exp(ej ) −

bCh exp(2ej ) θj , cCh + exp(ej )

(15)

 where the subscripts j = 1, . . ., N, emean = (1/N) ej is the mean potential drop, and K the same coupling parameter as that used in the experiments. V is the dimensionless applied potential; Rj the dimensionless equivalent individual resistance; Γ j the surface capacity. Rj ’s and Γ j ’s were taken slightly different values for each electrode to simulate the heterogeneity. R0 = (1/N) Rj · iF,j is the Faradaic current: iF,j (θj , ej )   Ch exp(0.5ej ) = + a exp(ej ) (1 − θj ). 1 + Ch exp(ej )

(16)

Eq. (14) is for the charge balance in the equivalent circuit (see Fig. 1b) of the electrochemical reaction; Eq.

Y. Zhai et al. / Physica D 205 (2005) 57–69

(15) is obtained from the simplified mass balance and kinetics. Since the coupling is electrical (through resistors), it only appears in the equation for the variable e. K = 0 represents uncoupled oscillators; K → ∞ yields maximum global coupling. The parameter values Ch = 1600, a = 0.3, b = 6 × 10−5 , c = 1 × 10−3 were used [23] to obtain dynamical features similar to experiments. Additional details of the simulations can be found in Ref. [24]. 3.3. Order parameter from individual signals If individual signals for all oscillators are available, an order parameter similar to Z1 (the Kuramoto order) can be found. An analytical signal ξ j can be obtained for each oscillator with the help of the Hilbert transform H(t) [25], ξj = (Ij (t) − Ij ) + iH(t), where H(t) =

1 PV π



∞ −∞

Ij (τ) − Ij  dτ. t−τ

(17)

(18)

The order parameter is the normalized sum of the analytical signals: N j=1 ξj . (19) r(t) = N j=1 |ξj | This order parameter has been previously used to characterize the degree of synchronization in experiments and simulations on populations of electrochemical oscillators [14,24].

4. Populations of nearly smooth oscillators In this section, we first calculate the order parameters with both Eq. (19) and the new proposed method (Eq. (9)) for populations of 64 smooth oscillators and compare them at different coupling strengths. Then we show the results of an experimental study on a larger size array of 100 elements. 4.1. Order measurement for array of 64 elements The time series of an individual current for the population of 64 smooth electrochemical oscillators

61

without added global coupling is illustrated in Fig. 2a. The oscillation can be considered approximately as ‘smooth’ because its angular velocity is almost constant, although there was a slight slowing down at minimum values of the current. Correspondingly, a sharp peak of the magnitude of fast Fourier transform (FFT) dominated around 0.45 Hz as shown in Fig. 2b while other peaks were much smaller (only the first two leading peaks were shown here). Each peak in this plot corresponds to a Fourier coefficient |Ak |; the peak of the lowest frequency (ωl ) is |A1 |, the one of the second lowest frequency (≈2ωl ) is |A2 |, and so on. The waveform of all individual oscillations looked similar to the one shown here and the Fourier expansion of the currents gave similar values of coefficients |Ak |. |Ak | exhibited a small change as experiments proceeded because of slight, unavoidable drift in our experiments. Therefore, for more appropriate representation of a ‘typical’ individual current at all coupling strengths, we used average |Ak (K)| obtained at each coupling strength K instead of |Ak (K = 0)| from any element randomly selected at the base case of K = 0, as the theoretical method implies. Without added global coupling, the frequency distribution was approximately unimodal with a mean of 0.4526 Hz and a standard deviation of 6.54 mHz. (The frequency was obtained by taking the slope of the phase versus time curve; the phase was calculated from the angle of analytical signal. For further details, see [24].) As a result, the mean current of 64 elements fluctuated irregularly with small amplitude (Fig. 2c) and in the plot of magnitude of FFT, only a small bump near the mean frequency appeared, i.e., |Bk | were quite small. Obviously, the order, |Z1 | = |B1 |/|A1 | = 0.005/0.030 = 0.17 was very low although nonzero due to the finite size effect. With some finite but still weak global coupling of K = 0.085, a more coherent state was reached. The waveform of the individual current did not change qualitatively as shown in Fig. 2e–f but the mean current now oscillated with large amplitude equivalent to that of the individual one. The latter can be better seen in the FFT plots (Fig. 2f and h) where |B1 | = 0.0295 and |A1 | = 0.033 indicates a high order |Z1 | = 0.0295/0.033 = 0.89. As mentioned above, the first mode of Zk , Z1 , should be the same as Kuramoto order r in definition, although Z1 is calculated with Eq. (9) while r is obtained with

62

Y. Zhai et al. / Physica D 205 (2005) 57–69

Fig. 2. Experiments of 64 smooth electrochemical oscillators. Left column: time series of electrode currents. (a) An individual current without added global coupling (K = 0). (c) The mean current at K = 0. (e) An individual current at phase synchronized state (K = 0.085). (g) The mean current at K = 0.085. Right column: time series of the magnitudes of the FFT of the corresponding currents shown in the right panels. (b) FFT of the individual current at K = 0. (d) FFT of the mean current at K = 0. (f) FFT of the individual current at K = 0.085. (h) FFT of the mean current at K = 0.085.

Eq. (9). Detailed comparisons of r and Z1 therefore are made at three different coupling strengths as shown in Fig. 3. The time series of order magnitudes in the left column shows that evolution of |Z1 | was very similar to |r| in all shown coupling strengths except for the fasttime-scale variation of |Z1 | was greatly suppressed. The elimination of the in-cycle fluctuation of the order for finite size population is one benefit of using Z1 instead of r because of the cycle-average feature in Eq. (9). The middle and right columns show r and Z1 in the complex plane, respectively. We can see that the two methods of order calculation result in differences in the phase of the order. Z1 should eliminate the spread of the phase of the order parameter according to the theory in Section 2; however, some small changes in the phase of Z1 still remained in the application due to a very small inaccuracy of the mean frequency Ωe we used in Eq. (9). The first three modes of |Zk | were calculated with the mean current obtained by taking the mean of the individual currents in the experiments of 64 smooth oscillators. The dependence of the temporal mean of the orders |Zk | (k = 1, 2, 3) on the coupling strength

K is shown in Fig. 4a in comparison with the mean Kuramoto order |r|. We can see that |Z1 |, the dots, have practically the same values as those of |r|, the line-connected circles, at all coupling strengths; slight deviations were present in a narrow region of K right after the critical coupling strength Kc . In Fig. 4b, the variance of |Z1 | shows enhanced fluctuation of the order near the Kc due to the finite size effect as that of |r|. The variances of |Z1 | were generally smaller than those of |r| because of the suppression of the in-cycle fluctuations. The results for a population of 64 simulated smooth electrochemical oscillators (Fig. 4c and d) further confirm the exchangeability of |Z1 | and |r| for the study on the relationship of order and coupling strength. In Fig. 4d, the multiple close peaks around Kc are a feature in the simulation of smaller oscillator population and reveals asymmetries in the synchronization process [24]. The temporal mean of |Z2 | and |Z3 | in the experiments and simulations are also plotted as a function of K in Fig. 4a and c, respectively. The three |Ak |’s had following relationships: |A1 | = 3 × |A2 | = 10 × |A3 | for the experimental data and |A1 | = 5.5 × |A2 | = 20 × |A3 |

Y. Zhai et al. / Physica D 205 (2005) 57–69

63

Fig. 3. Comparison of r and Z1 at three different coupling strengths in experiments of smooth population. Top row: K = 0. Middle row: K = 0.035 (Kc = 0.03). Bottom row: K = 0.085. Left column: time series of the magnitude of the orders. Thin curve: |r|. Bold curve: |Z1 |. Middle column: r in the complex plane. Right column: Z1 in the complex plane.

Fig. 4. Dependence of orders (both |r| and |Zk |, k = 1, 2, 3) on coupling strength (K) in experimental and simulated 64 smooth electrochemical oscillators. (a) The mean order (|r| or |Zk |) vs. K in experiments. Open circles with lines: |r|. Dots: |Z1 |. Stars: |Z2 |. Triangles: |Z3 |. (b) The variance of the order (|r| or |Z1 |) vs. K in the same experiments as in (a). (c) The mean order vs. K in simulations. (d) The variance of the order (|r| or |Z1 |) vs. K in the same simulations as in (c).

64

Y. Zhai et al. / Physica D 205 (2005) 57–69

for the simulation data. Both proved that the conditions we chose for smooth oscillators in the experiments and simulations were appropriate for categorizing the obtained oscillators as ‘smooth’; the first mode of Fourier transform of the individual current was evidently more dominant than the rest ones. In both experiments and simulations somewhat higher Kc ’s for |Z2 | and |Z3 | were observed. Furthermore, the increase of |Z2 | and |Z3 | after Kc were much slower than that of |Z1 | [13], especially in the experiments where the population was more heterogeneous ([max(f) − min(f)]/f = 0.04 for the experiments, 0.025 for the simulations). This is not surprising because the factor k in the exponential in Eq. (8) effectively increases the width of the phase distribution and this effect is stronger for larger values of k. 4.2. Effect of the array size on the synchronization process The behavior of an array of 100 elements were experimentally explored in which frequency distribution was similar to that of the array of 64 electrodes (both approximately normal distribution with the same standard deviations of 0.65 mHz; the distribution of the 100-

electrode array was somewhat smoother and more symmetric). Fifty six individual currents were recorded and the mean current was obtained from the potentiostat. Consequently, only Zk were calculated to characterize the degree of the order of the studied larger population. Simulations of 512 electrochemical oscillators in the smooth oscillation region were carried out to further investigate finite size effects. In simulations, r was used as the order parameter because all the individual information can be easily obtained by integrating the equations. The mean order parameters |Z1 | for arrays of 64 and 100 elements in experiments as a function of coupling strength are shown in Fig. 5a. At K = 0, |Z1 | was 0.12 and 0.2 for 100 and 64 oscillators, respectively. The Kc ’s and the increase of the order after the transitions in both cases were similar; only a very slight increase of Kc in the 100-array experiments was observed. The variance of the mean order for 100-array was again enhanced around Kc but with a smaller peak value compared with that of the 64-array. However, if we examine N × var |Z1 | as a function of K, as shown in Fig. 5b, the peak value is somewhat higher for the larger array, in agreement with the theoretical work of Daido on studies of finite size populations of oscillators [16,17]. All these size effects can be better seen

Fig. 5. Effect of system size on the relationship between the order (|r| or |Z1 |) and coupling strength (K) for smooth populations in experiments and simulations. (a) The mean order |Z1 | vs. K in experiments. Solid circles: array of 64 elements. Open circles: array of 100 elements. (b) N × var |Z1 | vs. K in the same experiments as in (a). (c) The mean order |r| vs. K in simulations. Solid circles: population of 64 elements. Open circles: population of 512 elements. (d) N × var |r| vs. K in the same simulations as in (c).

Y. Zhai et al. / Physica D 205 (2005) 57–69

in the simulation results where the population size is increased to 512. In Fig. 5c, the mean order |r| of simulated 512 oscillators at K = 0 is 0.038 while that of the simulated 64 elements is 0.11; both are in accordance with a N−1/2 scaling, i.e., √ assuming a Gaussian distribution of the Z, |Z| = π/N/2, giving 0.039 and 0.111 for N = 512 and 64, respectively. The transition of the order occurs at a little bit higher coupling strength for the larger size population and the increase of the order after Kc is smoother. In the variance plot (Fig. 5d), the finite-size related multi-peak character around the critical coupling strength is less evident for the larger population and a higher peak value around Kc is observed.

5. Populations of relaxational oscillators The applicability of the proposed method of order calculation is now examined with periodic relaxational oscillators both in experiments and in simulations. In experiments, the frequency distribution (f = 0.35 Hz, std f = 39 mHz) was flatter and broader than that of the smooth oscillators considered above. The relaxational character of the individual oscillation was

65

well confirmed by comparing the first three Fourier coefficients, |A1 | = 1.5 × |A2 | = 2.5 × |A3 |; |A2 | and |A3 | had more comparable values relative to |A1 | in the relaxational case. This also implies that the coupling function may include higher harmonic components. Therefore, we expect that the second and third modes of the Fourier expansion of time series of the relaxational oscillators may significantly contribute to the synchronization of the population as well as the first mode. Fig. 6a shows the order parameters (both |r| and first three modes of |Zk |) from the experiments on 64 relaxational oscillators as a function of K. Now the transition of the order can be characterized with a larger Kc and faster increase of the order above the Kc compared with the smooth population. Note that the values for |Z1 | at all coupling strengths were almost the same as |r|, although the minor differences between them were slightly greater than that in the smooth population. The same trend was observed in the plot of the variance of the orders |r| and |Z1 |, as shown in Fig. 6b. Another distinct feature of the relaxational population is that the dependence of |Z2 | on K followed closely to that of |Z1 | before the order was increased above 0.8. The evident deviation of |Z3 | from |Z1 | occurred at

Fig. 6. Dependence of orders (both |r| and |Zk |, k = 1, 2, 3) on coupling strength (K) in experimental and simulated 64 relaxational electrochemical oscillators. (a) The mean order (|r| or |Zk |) vs. K in experiments. Open circles with lines: |r|. Dots: |Z1 |. Stars: |Z2 |. Triangles: |Z3 |. (b) The variance of the order (|r| or |Zk |) vs. K in the same experiments as in (a). (c) The mean order vs. K in simulations. (d) The variance of the order (|r| or |Z1 |) vs. K in the same simulations as in (c).

66

Y. Zhai et al. / Physica D 205 (2005) 57–69

a smaller order value, approximately 0.65, which was still much higher than that in the smooth case. The corresponding simulation results are shown in Fig. 6c and d, where |A1 | = 3 × |A2 | = 6 × |A3 |. |A2 | and |A3 | were slightly smaller relative to |A1 | than those in the experiments. |Z2 | and |Z3 | are lower than |Z1 | above Kc ; however, both are still closer to |Z1 | compared with the simulated smooth case (Fig. 4c). The plot of the variance of |r| and |Z1 | versus K (Fig. 6d), exhibits multiple peaks around Kc because of the enhanced multi-step feature in the synchronization process of simulated finite-size relaxational oscillators.

6. Populations of phase coherent chaotic oscillators We also applied the proposed method to our experimental data of 64 weakly coupled chaotic electrochemical oscillators. The chaotic attractor was low dimensional and phase coherent [26]. There was again a unimodal frequency distribution of the uncoupled oscillators (f = 1.29 Hz, std f = 18 mHz).

Since the chaotic oscillator had time varying amplitudes and corresponding periods, we used Eq. (4) to calculate instantaneous Ak (n) for the nth cycle of each oscillator. The spatial average of the temporal mean of Ak (n) were used to calculate Zk ’s at each coupling strength. The first three modes of Ak followed approximately the relationship: |A1 | = 2.5 × |A2 | = 4.5 × |A3 |. The first mode was larger than the others, as was the case with the smooth oscillators. A sharp highest peak existed in the plot of FFT of the mean current even at K = 0. This further demonstrated that the chaotic oscillator was phase coherent and can be regarded as a smooth oscillator with small noise [27]. The dependence of the mean orders (|r| and |Zk |) on the coupling strength for the population of 64 electrochemical chaotic oscillators is shown in Fig. 7a. Again, |Z1 | had similar values of |r| at most of the weak coupling region before the phase synchronization was achieved at K = 0.11. However, above the critical point, small positive deviation of the |Z1 |’s from the |r|’s gradually developed. This may be caused by the low sampling rate (100 Hz) in the chaotic experiments. The simulation results on the periodic oscillator pop-

Fig. 7. Experiments of 64 chaotic electrochemical oscillators. (a) The dependence of the mean orders (|r| and |Zk |) on coupling strength K. Open circles with lines: |r|. Dots: |Z1 |. Stars: |Z2 |. Triangles: |Z3 |. (b) The variance of the order (|r| or |Zk |) vs. K. (c) Histogram of the magnitudes of instantaneous coefficients A1 (n)’s of an individual current at K = 0. (d) Histogram of the magnitudes of 64 instantaneous coefficients A1 (n) at the 50th cycle (n = 50 or t = 41.635 ± 0.595 s) at K = 0.

Y. Zhai et al. / Physica D 205 (2005) 57–69

67

ulations indicate that the calculation of the order from the global signal can be improved by using a higher sampling rate. As seen with the smooth oscillator population, |Z2 | and |Z3 | had much lower values than the corresponding |Z1 | or |r|. From Fig. 7b, we can see that the variance of |Z1 | had a peak around Kc as that of |r| and it was generally smaller than that of |r| in the entire weak coupling region. Fig. 7c and d shows the histogram of |A1 (n)| of an individual chaotic current and that of 64 instantaneous |A1 (n)| at a given n without added global coupling, respectively. Both the temporal and spatial distribution of |A1 | showed a bimodal feature at K = 0. This feature held true at other weak coupling strengths and in the distributions of |A2 | and |A3 |. This agrees with the finding that many phase coherent chaotic oscillators have two-banded characters [28]. The similarity of the two histograms also implies that the chaotic oscillators are ergodic; the statistical features of the instantaneous amplitudes of one oscillator are the same as those of the spatial distributions at a particular time. So, in a population of chaotic oscillators, the amplitude of each oscillator is different from each other at any specific time, but the variations of amplitudes cancel each other and therefore the spatial mean of Ak is constant in time as that in the periodic oscillators. This may explain why the method originally proposed for smooth periodic oscillators also works for our weakly coupled chaotic oscillators.

distribution, we found that the method can be used for a population of weakly coupled periodic oscillators with bimodal natural frequency distribution as well. On the other hand the method cannot be applied if either assumption (i) or (ii) is not met. For example, the proposed method failed in an application to strongly globally coupled relaxational oscillators, where dynamical clustering occurred [21]; the individual currents were greatly deformed by the strong coupling (K = 0.95) from a period-1 to a period-2 waveform while the mean current was still period-1. Another example in which the new method cannot be applied is that of identical synchronized (strongly coupled) chaotic oscillators, where the spatial mean of Ak varies with time. The success of the application of the method depends on the accuracy of the coefficients Ak . In the study of populations with greater heterogeneities, the mean Ak from a randomly selected small subgroup can be used for better representation of the ‘average’ element. Also, if drift occurs in the application and it changes all the individual elements the same way simultaneously, Ak (K) at each coupling strength can be used instead of Ak (K = 0) in Eq. (9). It should be also noted that high enough sampling rate in the measurement of the signals (the global/mean signal and at least one individual) is required in the application of the new method; otherwise, false in-cycle variation will be introduced.

7. Limitations of the order calculation from global signals

8. Conclusion

In the application of the proposed method of order calculation exclusively from global measurements, it should be noted that Eq. (9) is derived based on two important assumptions which make phase reduction possible [2]: (i) all the oscillators should be similar to one another and (ii) the coupling is sufficiently weak. The method was shown to be successful for populations of weakly coupled smooth, relaxational and phase coherent chaotic oscillators whose elements were similar. The only restriction of the method on the intrinsic frequency distribution is that its width is sufficiently small. Although all the results we showed above were obtained with a unimodal frequency

The generalized order parameters Zk can be obtained from global measurements with Eq. (9). From its modified version, Eq. (10), Zk is the ratio of the Fourier coefficients of the collective signal at the kth multiples of the entrained frequency to those of the individual signal. Therefore, Zk is a measure of order that expresses the resemblance of the collective behavior of the system to that of an individual element; it is done by comparing the coefficients of the each mode of the Fourier expansions of the two waveforms. The method had been successfully applied to the experimental and simulation data of weakly coupled smooth, relaxational and phase coherent chaotic electrochemical oscillators. The new method is capable of obtaining all Zk ; a com-

68

Y. Zhai et al. / Physica D 205 (2005) 57–69

plete examination of the series of Zk may give more detailed information about the synchronization process of the populations [13] or about the process of transient desynchronization [29]. A study on a larger array of globally coupled smooth electrochemical oscillators was performed based on the new method and the results show excellent agreement with our previous conclusions obtained with a smaller array in the study of the relationship of order on coupling strengths. A smoother increase of the order after the critical point was observed for the larger size population and the peak around Kc in the plot of N × var |Z1 | versus K had greater value as N increases, in accordance with the theoretical prediction by Daido [16,17]. The successful application of the proposed method of extracting the order from the global measurements verifies that it can be of great help in the study of synchronization or desynchronization processes of large populations of interacting oscillators in many fields including biology where measurements are limited to the global signals and only few local sites.

Acknowledgements This work was supported by the National Science Foundation (CTS-0317762). JLH, YMZ, and IZK thank Yoshiki Kuramoto for many fruitful discussions and for encouraging us to study populations of weakly coupled limit-cycle oscillators. HD is grateful to him for many enlightening discussions, comments, and criticisms over the last three decades.

References [1] A.T. Winfree, The Geometry of Biological Time, SpringerVerlag, New York, 1980. [2] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer-Verlag, Berlin, 1984. [3] A.K. Ghosh, B. Chance, E.K. Pye, Metabolic coupling and synchronization of nadh oscillations in yeast cell populations, Arch. Biochem. Biophys. 145 (1971) 319–331. [4] S. Danø, P.G. Sørensen, F. Hynne, Sustained oscillations in living cells, Nature 402 (6759) (1999) 320–322. [5] N. Wiener, Nonlinear Problems in Random Theory, MIT Press, Cambridge, MA, 1958.

[6] D.C. Michaels, E.P. Matyas, J. Jalife, Mechanisms of sinoatrial pacemaker synchronization: a new hypothesis, Circ. Res. 61 (5) (1987) 704–714. [7] R.A. York, R.C. Compton, Quasi-optical power combining using mutually synchronized oscillator arrays, IEEE Trans. Microw. Theory Tech. 39 (6) (1991) 1000–1009. [8] R.A. Oliva, S.H. Strogatz, Dynamics of a large array of globally coupled lasers with distributed frequencies, Int. J. Bifurcat. Chaos 11 (9) (2001) 2359–2374. [9] K. Wiesenfeld, P. Colet, S.H. Strogatz, Synchronization transitions in a disordered Josephson series array, Phys. Rev. Lett. 76 (3) (1996) 404–407. [10] A.T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators, J. Theor. Biol. 16 (1967) 15–42. [11] Y. Kuramoto, Cooperative dynamics of oscillator community, Prog. Theor. Phys. Suppl. 79 (1984) 223–240. [12] H. Daido, Generic scaling at the onset of macroscopic mutual entrainment in limit-cycle oscillators with uniform all-to-all coupling, Phys. Rev. Lett. 73 (1994) 760–763. [13] H. Daido, Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions: bifurcation of the order function, Physica D 91 (1996) 24–66. [14] I.Z. Kiss, Y.M. Zhai, J.L. Hudson, Emerging coherence in a population of chemical oscillators, Science 296 (2002) 1676– 1678. [15] H. Daido, Order function and macroscopic mutual entrainment in uniformly coupled limit-cycle oscillators, Prog. Theor. Phys. 88 (1992) 1213–1218. [16] H. Daido, Intrinsic fluctuation and its critical scaling in a class of populations of oscillators with distributed frequencies, Prog. Theor. Phys. 81 (1989) 727–731. [17] H. Daido, Intrinsic fluctuations and a phase-transition in a class of large populations of interacting oscillators, J. Stat. Phys. 60 (1990) 753–800. [18] O. Lev, A. Wolffberg, M. Sheintuch, L. Pismen, Bifurcations to periodic and chaotic motions in anodic nickel dissolution, Chem. Eng. Sci. 43 (1988) 1339. [19] W. Wang, I.Z. Kiss, J.L. Hudson, Clustering of arrays of chaotic chemical oscillators by feedback and forcing, Phys. Rev. Lett. 86 (2001) 4954–4957. [20] I.Z. Kiss, Y. Zhai, J.L. Hudson, Collective dynamics of a weakly coupled electrochemical reaction on an array, Ind. Eng. Chem. Res. 41 (2002) 6363–6374. [21] I.Z. Kiss, W. Wang, J.L. Hudson, Experiments on arrays of globally coupled periodic electrochemical oscillators, J. Phys. Chem. B 103 (1999) 11433–11444. [22] W. Wang, I.Z. Kiss, J.L. Hudson, Experiments on arrays of globally coupled chaotic electrochemical oscillators: synchronization and clustering, Chaos 10 (2000) 248– 256. [23] D. Haim, O. Lev, L.M. Pismen, M. Sheintuch, Modeling periodic and chaotic dynamics in anodic nickel dissolution, J. Phys. Chem. 96 (1992) 2676–2681. [24] Y.M. Zhai, I.Z. Kiss, J.L. Hudson, Emerging coherence of oscillating chemical reactions on arrays: experiments and simulations, Ind. Eng. Chem. Res. 43 (2004) 315–326. [25] D. G´abor, J. Inst. Electron. Eng. Part 3 93 (1946) 429.

Y. Zhai et al. / Physica D 205 (2005) 57–69 [26] I.Z. Kiss, J.L. Hudson, Phase synchronization and suppression of chaos through intermittency in forcing of an electrochemical oscillator, Phys. Rev. E 64 (2001) 046215. [27] A.S. Pikovksy, M.G. Rosenblum, J. Kurths, Synchronization: A Universal Concept in Nonlinear Science, Univeristy Press, Cambridge, 2001.

69

[28] J. Davidsen, I.Z. Kiss, J.L. Hudson, R. Kapral, Rapid convergence of time-averaged frequency in phase synchronized systems, Phys. Rev. E 68 (2003) 026217. [29] P.A. Tass, A model of desynchronizing deep brain stimulation with a demandcontrolled coordinated reset of neural subpopulations, Biol. Cybern. 89 (2003) 81–88.