Synchronization of chaotic systems

Because of its positive Lyapunov exponents, an isolated chaotic system synchronizes with nothing else. The synchronization of two or more chaotic systems is actually another version of the usual thought experiment where one imagines two initial conditions starting close to each other in phase space. For two separated chaotic systems, trajectories starting at close initial condition will diverge at first. It was counterintuitive to have the two initial conditions case converge to the same trajectory. We now realize that two identical uncoupled chaotic systems have twice as many positive Lyapunov exponents as either of the systems by themselves. If, on the other hand, we add some sort of coupling between the isolated systems, we can alter the Lyapunov spectrum-depending on the coupling and the specific systems, we could decrease or increase the number of positive Lyapunov exponents. In our original synchronizing system, we went from a total of two positive Lyapunov exponents in the combination of two uncoupled systems to one positive Lyapunov exponent in the pair of coupled systems. We show below how we originally handled this situation mathematically so we could calculate the stability of the synchronized state. As we go through this article, we will also show how this approach matured and became general with a geometric way to view the synchronization in large networks of coupled oscillators.

By 1988, we had entered the world of nonlinear dynamics and published several papers on nonlinear of spin waves in magnetic materials.^1,2 While the materials we were studying were important for radiofrequency electronics, the difficulty of the theory of spin waves prevented us from translating nonlinear dynamics into useful applications. While the Navy was supportive of basic research, at some point we needed to be able to show how our basic science might affect future applications. We needed to come up with a more immediate application of chaos. One idea that came to mind was some type of message masking or hiding, where a chaotic signal covers or hides a real message in a signal. Had we known a lot about cryptography and communications we might not have been so naive to think that this could be done simply. But sometimes ignorance is a blessing in another way. We talked for months on and off about ways to do this, but came up with no good, testable version.

Early in 1989 after a dynamics conference, we decided that combining two identical chaotic systems in a particular way would be worth a try. The idea was to take a signal from a component of one (the transmitter) and send it to a duplicate system (the receiver) where the receiver was missing the part of the system that the signal came from in the transmitter, but this missing part was compensated for by using the received signal at points in the receiver where the missing part would have supplied the signal. This is a bit confusing in words but clear in a diagram in Fig. 1.

If we have chosen the correct variable to drive the response, then y′ and z′ subsystems will converge to y and z as the systems evolve together on identical trajectories, i.e., the response variables will synchronize to the drive variables. We can think of this situation as the differences $| y−y′ |$ and $| z−z′ |$ going to zero as time progresses. Below we will make this rigorous and quantitative.

A type of system like Fig. 1 that we first tried was actually a 2D map system constructed from logistic maps in a very simple way

where (1) is the drive and (2) is the response. For α = 3.9 and β = 1.1, we see the chaotic attractor of the drive in Fig. 2(a). If we start y and z at different values, we see the plot of ln(|z_n − y_n|) in Fig. 2(b) which shows convergence of z_n → y_n. This is not completely trivial since setting β = 1.6 gives Fig. 2(c) for the difference in y and z, which shows there are regions of the attractor where there are instabilities repelling y and z although the systems still synchronize.

Although we were encouraged by this little calculation, we realized that we needed a more realistic drive-response system to display convergence and hence synchronization of one subsystem with another. That meant we had to set up differential equations that would more realistically mimic actual physical systems. And after making that last statement, we realized we eventually would make this whole approach more convincing if we could find a physical system that could be built to show synchronization.

We set up the usual Runge-Kutta (4,5) numerical solver with some unusual (at the time) vector fields made from an original vector field parts of the original (a subsystem) linked to the original by a drive signal from the missing part. We did this for the Lorenz and Rössler systems along with several others. The equations more clearly display what the word description says (note drives and responses always have the same parameter values)

(a) is the Lorenz drive- response system using x driving;

(b) is the Rössler drive-response system using x driving;

(c) is the Rössler drive- response system using y driving.

Immediate observations were that for parameters that gave chaotic behavior x-driving worked to synchronize Lorenz, y-driving worked to synchronize Rössler, but x-driving did not work to synchronize Rössler. For example, if we start a Rössler drive and a y-driven response from different x and z values then let them evolve trajectories as we see in Fig. 3(a). For a Lorenz system with x-driven response and different y and z values, we get Fig. 3(b). This synchronization has come to be called identical synchronization since the result (at infinite time) is that the two systems (drive and response) have exactly equal dynamical variables.

All this is not obvious, but it emphasized that we needed a stability theory to explain why the phenomenon worked or did not work for all systems we tried. The usual Lyapunov exponents for the original Lorenz or Rössler systems give no information for the drive- response systems. But we can follow a similar process to develop a variational equation for our cases. For the Lorenz with x driving, we set y′ = y + δy and z′ = z + δz and subtract the y-z subsystems of the drive and response. This generates a variational equation for the y-z response subsystem

In a manner analogous to finding the usual Lyapunov exponents, we evolve Eq. (6) along with the associated drive and calculate stability exponents.⁵ Because these are different from the usual Lyapunov exponents and they depend (are conditional) on the x signal we call these exponents conditional Lyapunov exponents.⁶ For the Lorenz and Rössler systems, we get the exponents as shown in Table I.

We can see immediately that the Rössler x-driving y-z subsystems is unstable and will not synchronize since one of the conditional Lyapunov exponents is positive.

It is easy to see that this can apply to dynamical equations and systems of any dimension. And one can also generalize to systems using multiple drive signals and maps.⁴ 

As the above was developed, we worked on implementing this approach to synchronization in a real system. Thomas Carroll had long been interested in analog computers and circuits after seeing one in operation by his undergraduate professor. The circuit used operational amplifiers, capacitors, and resistors to simulate the equations for a bouncing ball with damping, and the circuit output was displayed on an oscilloscope. There was now a good reason to attempt to build an analog computer circuit: chaotic synchronization. Although other physical systems might have worked, a circuit can often be split into sub-circuits easily just like in the synchronization technique that was developed above. Tom found a chaotic circuit in the literature that had been developed by Newcomb and Sathyan⁷ and Newcomb and El-Leithy⁸ of the University of Maryland. The circuit has two lobes of unstable foci analogous to a Lorenz system. The circuit hopped between the lobes, but the hopping occurred at values of the dynamical variables (the voltages and currents) that made the system hysteretic. Tom built a pair of circuits: a drive circuit and a response circuit that synchronized to the drive circuit. This final experimental confirmation was important. It showed that the technique could be implemented in a real system with noise and parameter mismatch.

As for keeping the Navy interested in our work, we did a live demonstration of chaotic synchronization and signal masking in front of our lab director. As a plasma physicist, he was skeptical of this new field of chaos, as he believed plasma physicists were already studying nonlinear dynamics, but he was intrigued enough by our new work to allow us to continue. Tom Carroll even believes that chaotic synchronization directly caused him to be hired as a permanent employee.

Most of the above (except the section on the 2D logistic map synchronization) appears in our paper³ as well as our second paper⁴ along with a lot of other information on this approach to synchronization.

In the world of science and mathematics, new proofs, interesting phenomena, physical behavior, etc., are often discovered more than once, often without knowledge of the previous work beforehand. This happened to us. It turns out that synchronization of chaotic systems of the kind we presented above (identical synchronization) was discovered several times, often using different approaches. We list in reverse historical order publications and other information that we discovered over the years after our initial work. The reverse order here is not the order in which we actually discovered them.

1. In 1989, Aranson and Rul'kov published a paper analyzing synchronization zones in multidimensional dynamical systems.⁹ This paper studied the parameters and bifurcation diagrams of synchronized, driven microwave oscillators. This is an early study of synchronization in nonlinear multidimensional systems.

2. In 1989, Volkovskii and Rul'kov¹⁰ studied the bifurcations which lead to stochastic locking in coupled, self oscillating systems with piecewise nonlinear vector fields. They studied both mutual and unidirectional coupling. As in so many Russian papers on this subject, the authors also built experiments to test the theory.

3. Afraimovich et al. published a now famous paper in 1986 (Ref. 11) on the mutual synchronization of two time nonautonomous, chaotic oscillators. When we first came across this paper, we were confused by the label of “stochastic” and only later realized that it was identical synchronization through the dynamics similar to our system. At a conference, Rabinovich assured us that their approach was actually a general version of what we had done. Their equations of motion in Ref. 12 were

The mutual couplings c₁ and c₂ provided a dissipative dynamics pulling the two systems into synchrony despite their chaotic behavior if the coupling was large enough. If we set c₁ = 0 and let c₂ → ∞, we turn the first line in Eq. (7) into a drive and the second line into a response with y₂ slaved to y₁ just like our original approach. So Rabinovich was right. However, he was gracious in saying that they missed the possible applications of such synchronization, which is what our paper provided.

4. Pikovskii published two papers that had synchronization as their theme in 1984.^13,14 By coupling chaotic systems, Pikovskii found he could synchronize two chaotic attractors.

5. Starting in 1983, Fujisaka and Yamada published a series of four papers on synchronization in chaotic systems in Progress in Theoretical Physics.^15–18 These papers investigated in much detail several aspects of symmetrically coupled identical chaotic systems and their synchronization. The systems were thought of as locations along spatial axes with a diffusive type coupling through a rdiscrete analog of the usual diffusive coupling in continuous spatial systems, that is the spatial Laplacian operator $∇2$⁠. This yields a set of discrete, coupled ordinary differential equations (ODEs)

where D is the diffusion constant which acts like a coupling constant linking the dynamics of different spatial locations and, because of its dissipative nature, the coupling terms make the synchronization of different spatial dynamical locations possible for certain values of D. Fujisaka and Yamada's work is the earliest we know of for showing chaotic synchronization.

Once one see the interesting phenomena emerging from complex systems constructed from simpler components (e.g., drive- response or coupling systems), one can easily imagine constructing many new and interesting ensembles of simple, chaotic, or nonchaotic components. As the field of synchronization of chaotic systems began to grow, we would characterize the use of various coupling schemes and constructed systems as becoming a “cottage industry” in the world of nonlinear dynamics. The literature on all the possible approaches is huge and we will not attempt to cover it here. We will give only some highlights from our own work and some others that introduced new interesting phenomena discovered by constructing various coupled systems.

One of the first extensions we did on the drive-response approach is to cascade the systems so as to reproduce the incoming drive signal when the drive and response synchronize.¹⁹ This is shown in Fig. 4 below. Basically, we form two drive-response systems so that the second one in the cascade has the original drive component as part of the response, thus reproducing the incoming drive signal.

We also showed that we could synchronize systems in such a way that both the drive and response filled a subset of the phase space with non-zero measure and they would synchronize. The goal here was to find more random looking functions which would serve communications better, but still be deterministic so they could be synchronized. We accomplished this by starting with an n-dimensional phase space and a linear, expanding map L (this provides the positive Lyapunov exponents) that acts on the phase space points and couples the components, i.e., the output is a linear combination of the input components. Thus, if x(m) is the n-dimensional vector of components, where m is the time step, then we write that the linearly transformed vector is just L x(m). However, we must add another function acting on this expanded vector to fold it back into the original volume so the system will remain finite. We do this by using a piecewise linear folding function F. Fig. 5 shows a simple folding function analogous to a Bernoulli map. Another version is just the map F(x) = (x₁ mod k, x₂ mod k, …, x_n mod k).

Then the iterated coupled maps are just written as

where $xd$ is the drive vector, $xr$ is the response vector, and $xr=(x1d,x2r,...,xnr)$⁠, i.e., we drive only with the 1st component of $xd$⁠. Obviously, this can all be generalized (see Ref. 20). The important point is that a very random attractor can be synchronized with its driven partner, the response. Fig. 6 shows a drive attractor for a 3D system. Note that the attractor fills a 3D cubic volume. The synchronized response is exactly the same.

An advantage of this synchronized system is that the autocorrelations drop off to near zero after one step. The synchronization decay is very fast with differences between the drive and response down by a factor of 1000 in four steps. However, the spectrum of the system is broad-band and much like white noise. We built a circuit based on these ideas to show proof of principle. In addition we could mix in a message with the drive signal using the more robust XOR between the signal and the drive signal. This was easily recovered on the response side (see Ref. 20).

An important finding by Peng et al. was that systems with more than one positive Lyapunov exponent could still be synchronized using one drive (scalar) signal.²¹ There was speculation that one would need as many drive signals as there were positive Lyapunov exponents in the original system. The paper by Peng et al. cleared that controversy up. Using their formulation as shown below in Eq. (10), we found another interesting behavior of synchronized systems.

A postdoc working with us (Johnson) discovered that we could optimize the linear coupling between the drive and response by Peng et al. by tuning the coupling matrix to have the best properties to rapidly and stably synchronize the response.²² This resulted in the ability to synchronize the response so well that even when the drive's parameters changed by as much as 50% and it underwent a bifurcation (e.g., chaotic to periodic behavior), the response was carried with it undergoing the same bifurcation despite having different parameters. The general scheme looks like this

where B and K are vectors which can be tuned to guarantee the greatest stability of the response following Johnson et al.²² We developed a four-dimensional experimental system (a circuit) which has the dynamics of the following ODE system:

where ρ is the bifurcation parameter varied between 0.05 and 0.91 and Θ is the Heaviside step function. With Eq. (11) as the drive and the same with the coupling term as the response, we ran experiments in which ρ is varied by 50% of its value through a set of bifurcations. With the coupling term tuned to give maximum stability for the parameter matching case, we get the progression of attractors shown in Fig. 7 going from periodic, to multiperiodic, and on to chaotic dynamics. The response (with ρ held constant) follows the drive through the bifurcations. We called this behavior imposed bifurcations.

We note briefly here that for the last 10–15 years, many control theory approaches have been proposed to ensure the synchronization of a drive-response pair. This work is a bit far afield for CHAOS and is certainly too large to include in any concise fashion. We simply note that most control theory techniques have been used to successfully synchronize two chaotic systems, such as sliding mode control, linear matrix inequalities (LMI), simple coupling, impulse stabilization, and many others along with several applications of observer theory to chaotic synchronization. One important contribution of control theory approaches is that for some systems, it has been possible to prove global synchronization. This is a much stronger result than linear stability theory, such as conditional Lyapunov exponents above.

A few years after working on synchronization of two chaotic systems, we ventured into the realm of synchronization in networks of oscillators. We began with the case of oscillators coupled using general forms of linear diffusive coupling, e.g., nearest neighbor linear diffusive coupling. We examined both periodic and non periodic boundary conditions. These systems are straightforward to analyze because we can diagonalize the coupling matrix and reduce the equations to easily interpreted forms. However, the ingredients for much more general scenarios are in these simple cases and we find some dynamical surprises along with the seed for a very general result in cases of more complex networks of oscillators.

In the case of linear diffusive coupling, we can write the coupling matrix as a shift-invariant or circulant matrix. Here are the equations of motion for the near-neighbor coupling case along with the coupling matrix

where $x=(x1,...,xN)$⁠, $F(x)=(F(x1),...,F(xN))$⁠, σ is the coupling strength, and each x_i is an n- dimensional dynamical variable at one node/oscillator in the network and F is an n-dimensional vector field, identical for each x_i. The coupling is through the nN-dimensional matrix G = $A⊗B$⁠, where B is an n × n matrix picking out which components in each x_i are coupled to its neighbor and the coupling matrix G is an N × N matrix given by

If s(t) is the synchronized motion of the system, that is, when all oscillators in the network are synchronized and follow the trajectory s(t), then we can view this as all the motion settles onto a synchronization manifold of n dimensions in an nN dimensional space. The synchronization manifold is a flat hyperspace defined by all the equations of the type x_i = x_j for i,j = 1,.,N. What is important is whether s(t) is a stable motion at each value of σ. This can be calculated from the variation of the equation of motion, Eq. (12), along the s(t) trajectory given by the linear stability equations of motion

The Lyapunov exponents of Eq. (14) determine the stability of the system, but as it stands, Eq. (14) does not give a clear picture of the stability of the synchronous state. This is because the perturbations which leave the state in synchrony, i.e., perturbations which have all equal coordinates, keeping the motion on the synchronization manifold, are mixed in with the perturbations which desynchronize the system. The former do not matter in the calculation of the Lyapunov exponents for the stability of the synchronous state. Only the perturbations which desynchronize the system determine the stability of the synchronous state. To separate the two sets of parallel and transverse perturbations effects on stability, we need to make a transformation of the variational equations which separate the two. That transformation is the diagonalization of the coupling term since the individual vector field part (F(x)) is already diagonal as a multiple of the identity matrix.

The transformation which diagonalizes G is the discrete Fourier transformation. This is the case for all shift invariant (circulant) matrices. Recall the first term in Eq. (14) is really $1N⊗DF(s(t))$⁠, where 1_N is the N × N identity matrix. Since we are diagonalizing the first terms in the direct products, this term is unchanged. Writing the new variational dynamical variables as $ζ=Tξ$ with $T=S⊗1n$⁠, where S diagonalizes the circulant matrix A, we have a set of uncoupled ODEs as one would get with any mode decomposition in a linear equation

that is, λ_k's are the eigenvalues of matrix A.

At this point, Heagy, a postdoc working with us on this problem who derived Eq. (15) made several important observations. One is that the eigenvalue λ₀ = 0 is associated with the variations in the synchronization state, i.e., variations that do not desynchronize the oscillators. The eigenvector for this state is (1,1,.,1)/√N. All other eigenvalues are associated with transverse variations which do desynchronize the oscillators. It is the stability of the latter variations that matter for the study of the stability of the fully synchronized state. If all the variations of the k ≠ 0 eigenmodes have negative Lyapunov exponents, then the synchronized state is stable. If one or more are positive, the synchronized state is unstable. And the most important observation is that the variational equations (Eq. (15)) are of the same form for all modes, except for the factor of λ_k multiplying the coupling strength σ. So if we solve for the Lyapunov exponents of one mode as a function of coupling σ, we automatically have the exponents for all modes by just scaling their couplings accordingly. It turns out that all these observations generalize to the case of arbitrary networks, including nonsymmetric and weighted coupling cases; more on that later.

When we applied this theory to Rössler oscillators coupled to nearest neighbors through the x variable, we discovered a new phenomenon called a short wavelength bifurcation.²³ The origin of this name and the phenomenon is easily seen by considering the point above about each mode having the same functional form of Lyapunov exponent vs. σ except for the argument being scaled by the eigenvalues λ_k. Suppose we have a ring of four Rössler oscillators coupled through their x variables to nearest neighbors, then we will have two modes present (one of which is two-times degenerate) given by λ₁ = λ₃ = $4 sin2(π4)=2$ for the longer wave length mode and the other given by λ₂ = $4 sin2(π2)=4$ for the shorter wavelength. By wavelength here we are thinking of the form of the perturbations associated with each mode around the ring. Hence, the shorter wavelength mode acts as though the coupling is twice as large. This leads to a plot of the maximum Lyapunov exponents μ for these transverse modes vs. coupling as in Fig. 8.

The curves are the same, but are scaled differently in their argument. Note that for the x-coupled Rössler oscillators increasing the coupling can desynchronize the oscillators. The shortest wavelength mode goes unstable first. This is somewhat counterintuitive since most instabilities in diffusive coupling come from decreasing coupling strength until the longest wavelength mode becomes unstable. In a sense, this is an extreme form of a Turing bifurcation.²⁴ The loss of synchronization stability for this case at high coupling is consistent with Table I which shows the x-driven response is unstable since the drive- response systems is the same as a diffusive, one-way coupling in which the coupling becomes infinite. A short wavelength bifurcation or instability also implies that some networks of the oscillators never have stable chaotic synchronization. In the case of the ring configuration above for the parameters used for the Rössler there is a size effect²⁵ for which adding more oscillators to the ring will eventually prohibit any synchronization because the highest mode will undergo a short wavelength bifurcation as coupling increases before the lowest mode becomes stable.

In 1991, Pikovsky and Grassberger²⁶ found that symmetry breaking bifurcations, as the coupling strength is changed between initially identical, synchronized chaotic one-dimensional systems, can cause certain orbits on the attractor to become unstable causing excursions out of the synchronous state. They exhibited the local strange invariant set which caused the excursions. This was the first evidence noted in the literature for a phenomenon latter called attractor bubbling. In 1994, Ashwin et al.,^27,28 published papers which also showed that as one changes a parameter (e.g., coupling strength) that decreases the stability of the synchronous state of coupled, chaotic oscillators the dynamics and associated attractor geometry go through a series of changes. The first occurs when the synchronized attractor loses asymptotic stability. For example, certain periodic orbits in the synchronized chaotic attractor become unstable. At this stage, if there is a coexisting attractor that the system will eventually be attracted to as the parameter changes, the basin of attraction for the chaotic attractor becomes riddled with points from the basin of attraction of the coexisting attractor.^29–31 The second change occurs at another parameter value the attractor becomes a chaotic saddle, no longer a true attractor. And the third change in behavior as the parameter continues to be varied is that the saddle becomes a chaotic repeller.

During the first stage if there is a little noise or parameter mismatch on the synchronized chaotic systems, then the systems will occasionally diverge from the synchronized state, then return to it. This constant fluctuation is what is usually called “Attractor Bubbling.” It occurs when the synchronized trajectory wanders near one of the unstable sets (e.g., a periodic orbit) and the system is temporarily kicked out of sync. If there is no coexisting attractor, the system will just bubble. Ashwin et al.,^27,28 demonstrated the bubbling behavior in two coupled electronic circuits.

In studying the short-wavelength bifurcation/loss of synchronization stability, we noted that there are two coexisting attractors near the desynchronization threshold for the highest mode. Ott suggested that this situation might display riddled basins of attraction.³² They are arranged spatially like the wave of instability of the highest mode and a symmetric partner to it. In Fig. 9(a), we show the chaotic attractors for the original ring of four Rössler and the two symmetrically arranged periodic attractors which coexist with the synchronized chaotic attractors, and the attractors left after the coupling parameter is set at a value where the synchronized state is unstable. We automated our original four-Rössler experiment so we could scan a part of the phase space which showed the highest spatial frequency of deviations from the chaotic synchronized state (the spatial frequency of the periodic oscillator placement). This gave the first experimental picture of a riddled basin of attraction.³³ 

An important concept came out in the early 1990s which captures possible relationships between a drive and a synchronized response or, more generally, a stable response. This was called generalized synchronization³⁴ because the response is stable (conditional Lyapunov exponents are all positive), but is a different system from the drive, this is not identical synchronization as above, but a generalized form of it, hence the name. A criterion of generalized synchronization is that there is a function, say, φ, from the drive to the response so each point on the drive attractor x(t) is mapped to a point y(t) = φ(x(t)) on the response. The stability of the response is necessary for such a function to exist, but is not sufficient (see below). In a followup paper Abarbanel et al.³⁵ provided a tool to check easily for the possibility of generalized synchronization, the auxiliary system method. This consists of starting with two identical responses and driving both with the same system. If the systems are stable and there is only one basin of attraction, then the two responses will synchronize and plotting the variables of one against their counterparts on the other response should yield the signature straight, diagonal line of synchronization. There have been many followup studies and discoveries of generalized synchronization and to this day the concept remains interesting.

We mentioned above that there are times when the stability of the response is not enough to guarantee that the function φ exists. An example is a periodic system driving a nonlinear system that has gone through a period-doubling bifurcation. This means that the drive goes around its attractor twice while the response goes around once. Hence, there must be two points on the response which correspond to each point on the drive. Of course, mathematically, a function cannot do that. This situation was found to exist (surprisingly) in a chaotic response by Parlitz et al.³⁶ In the latter case, the drive and response were chaotic, but the response was entrained in a 2:1 way. For example, an unstable periodic orbit in the drive had a period-doubled unstable orbit in the response associated with it.

We mentioned above that we were initially somewhat naïve about the difficulties of actually using chaos for a communication system. Identical synchronization between two chaotic systems in the presence of noise becomes impossible for noise whose amplitude is greater than approximately 10% of the signal amplitude. The chaotic response system mixes noise and signal in a nonlinear fashion, making it impossible to separate the two by conventional methods, such as filtering or correlation.

Sterling made one attempt to reduce the effect of noise on a chaotic shift map.³⁷ Sterling considered maps of the form

where x is a driving signal. A noise signal added to x_n could cause y_n to exceed the threshold for the modulus function, which would show up as a sudden change in the value of y_n; lack of a corresponding change in x_n indicated that the change was due to noise.

We realized that conventional noise reduction schemes depend on some form of time integration over a long time, the signal should integrate to some finite values, while the noise integrates to zero. Even simple filters sum values of a signal from a range of different times. We designed a chaotic system that operated on two different time scales to take advantage of time integration.

The two-frequency chaotic system was based on a piecewise linear version of the Rossler system.³⁸ Denoting the piecewise linear Rossler vector field as R_pl allowed the two frequency Rossler system to be represented schematically as

where f = fast, s = slow, and pl = piecewise linear.

The slow vector field, R_s, was similar to R_pl, but it included extra damping so it did not oscillate on its own. The signal that was transmitted was a linear combination of the fast and slow signals.³⁹ The transmitted signal was coupled into an identical response system using linear diffusive coupling. The output of the response system was the first component of the slow variables. As the difference between the time constants τ_fast and τ_slow increased, the low frequency part of the response system synchronized more and more closely to the low frequency part of the drive system, even for noise larger than the drive signal.

The reason the low frequency synchronization was so stable was traced to the presence of long period unstable periodic orbits with purely real Floquet multipliers. The presence of such orbits may be more common than is realized in other physical systems- in this same paper on chaotic systems, it was also shown that low period unstable periodic orbits (UPOs) with real Floquet multipliers were also present in a published neuron model.⁴⁰ Such UPOs may have been responsible for the noise robustness of information transmission in the neuron model.

Above, we mentioned our first venture into synchronization of oscillators in a network for a shift-invariant coupling matrix. As we will show, the analysis of that case leads to a very general analysis of synchronization in arbitrarily coupled networks of oscillators. But first we should mention two other early and important works on synchronization in networks of coupled, chaotic oscillators.

In 1995 and 1996, Gade et al.⁴¹ and Gade,⁴² respectively, showed that synchronization in certain types of networks could have a general form for the analysis of the system's stability. In the 1995 paper, they showed that coupled maps on tree-structured networks would have better synchronization robustness than the same chaotic maps on a regular (e.g., nearest neighbor coupled) lattice. They also exhibited the results for synchronization in the limit of an infinite lattice. They also showed that solving the eigenvalue problem for the coupling matrix enabled a very general solution for the tree network coupling synchronization problem. The 1996 paper analyzed synchronization of chaotic maps in networks which had nonlocal coupling, for example, smallworld networks. Again the eigenvalue spectrum of the coupling matrix was the key to finding relations between the coupling constant and synchronization stability.

In 1998, Hu et al.⁴³ showed that diffusively, linearly coupled oscillators on lattices with periodic boundary conditions can be analyzed for stability of the synchronous state using an eigenvalue approach. They introduced the idea of a matrix that picks out which components of each oscillator will be coupled to others and that the regions of the complex plane encompassing the eigenvalues of the coupling matrix will change drastically with choice of oscillator coupling matrix. They showed that the regions of stability for a specific type of coupling matrix remained the same as the lattice of oscillators is increased.

In 1998, we returned to the topic of synchronization of chaotic systems in networks. We started with the above case of circulant coupling matrices and quickly decided that we could do stability analyses for an arbitrary network for which the coupling matrix was diagonalizable. We started with the very general set of equations of motion for the coupled system

where everything has the same meaning except that we have a coupling function H: R^m—> R^m which makes this formulation of the problem very general. We started by assuming the coupling G is bidirectional, i.e., G is a symmetric, real matrix and has zero row sum, which enables global synchronization of all oscillators.

Out starting point is the important observation of scaling as above in Sec. III D. By diagonalizing G, we would end up using a scaling argument to generate the stability of each eigenmode from a calculation of only one, as in Eq. (15), but where λ_k is an eigenvalue of the coupling matrix under study. However, we realized that it was better to turn this argument around. Instead of generating a separate curve for each eigenmode, we could have a function which was the maximum Lyapunov exponent for the variational equation

Then this function would represent the maximum Lyapunov exponent as a function of the “coupling” α. To find the stability of a particular eigenmode (λ_k) at a particular value of coupling σ, we locate the point on the α axis where α = σ λ_k. An example of this function for the Rössler system where H merely picks out the x component to use in the coupling is shown in Fig. 10(a). The important observation is that Eq. (17) does not depend on the coupling matrix G, but only on the dynamics F and the coupling function H. This means that for symmetric coupling G, Fig. 10(a) shows a stability function that applies for any bidirectionally coupling network with F and H the same. We named the stability function the Master Stability Function.

We realized that having a symmetric coupling matrix G is not necessary so long as G is diagonalizable. In general, the eigenvalues will be complex so we can solve for the MSF over the complex plane using a complex extension of Eq. (17), $ζ̇=[ DF(s(t))+(α+iβ) DH(s(t)) ]⋅ζ$⁠. To determine the stability of an eigenmode, locate the position of the complex eigenvalue (λ_k) and determine if it is in the stable or unstable region of the MSF (see Fig. 10(b)). The value of the MSF is that is separates the dynamics (the stability) from the network (the eigenmodes and eigenvalues) so each can be calculated separately and one calculation of the MSF can be used for the stability analysis of any network configuration with the same dynamics.

There is much more that can be done with the MSF, but we stop at this point. There is a separate article in this issue that explores many facets of synchronization in networks of oscillators. We recommend any reader interested in synchronization or collective phenomena in networks explore that article.

We complete this overview of synchronization with the brief mention of two areas of future work using synchronization. One is parameter estimation, which is related to a similar problem in control theory, the other is a continuation of synchronization into the larger topic of collective behavior and dynamical patterns in networks.

An early exploration of how synchronization can be used to estimate parameters was the paper by Parlitz.⁴⁴ In this paper, he starts with a dynamical system (chaotic or not), $ẋ=f(x,p)$⁠, where p is a vector of parameters for the system, and then shows how to construct a system driven by a component of the x system, $ẏ=g(s,y,q)$ so that y synchronizes with x when q = p. Thus, we can get an estimate of p by synchronizing y to x. The trick is to choose s as a function of x properly so we get that synchronization.

More recently, a more sophisticated approach has emerged which takes into consideration the type of stability of the synchronizing system (y). This is work by Abarbanel et al.^45–47 The basic idea is similar to Parlitz' above, but more care is taken to avoid mathematical difficulties and the method is generalized to more complex situations, including those involving time series and data using sophisticated regularization approaches. We do not have sufficient space here to even crudely cover all the ideas in this work, but we note that it sets the tone for much more development of the use of synchronization for parameter estimation and system analysis. We suggest that this is an area that has much potential for more development.

With the exception of parameter estimation and data assimilation (above) the synchronization of two systems by driving and coupling has been well explored. However, over the last decade, a new area of more general types of collective behavior has developed out of the synchronization studies. This is the area of synchronization patterns in networks of oscillators. By this, we mean that a network of identical oscillators may not be fully synchronized as in the above case of the MSF, but subsets of oscillators can synchronize to each other within the subset. This has been described in various terms, but the phrase cluster synchronization (CS) seems to be most appropriate and common. The term cluster is used to describe the collection of nodes (oscillators) in a network that are synchronized. They are not necessarily in close proximity it is just that their dynamics are the same. Many versions of CS have been displayed.^48–56 These involve many different scenarios (unidirectional coupling, time delays, special network structures, etc.), many of which are engineered to yield certain cluster synchronization patterns. Recently, the use of symmetries^57,58 and group theory has been extended to complex networks to discover patterns of cluster synchronization which are not obvious or which are too complicated to be discerned by eye or even human means.^59–61 

The area of cluster synchronization of other collective patterns is large and still growing. It opens up the older global synchronization (where all oscillators are synchronized) to a much vaster range of possible types of dynamical order. This area should be fruitful for the future as larger and more complex networks are explored and more general approaches to finding and/or engineering synchronization patterns are discovered.