# The stability of adaptive synchronization of chaotic systems

###### Abstract

In past works, various schemes for adaptive synchronization of chaotic systems have been proposed. The stability of such schemes is central to their utilization. As an example addressing this issue, we consider a recently proposed adaptive scheme for maintaining the synchronized state of identical coupled chaotic systems in the presence of a priori unknown slow temporal drift in the couplings. For this illustrative example, we develop an extension of the master stability function technique to study synchronization stability with adaptive coupling. Using this formulation, we examine local stability of synchronization for typical chaotic orbits and for unstable periodic orbits within the synchronized chaotic attractor (bubbling). Numerical experiments illustrating the results are presented. We observe that the stable range of synchronism can be sensitively dependent on the adaption parameters, and we discuss the strong implication of bubbling for practically achievable adaptive synchronization.

We consider an adaptive scheme for maintaining the synchronized state in a network of identical coupled chaotic systems in the presence of a priori unknown slow temporal drift in the couplings. Stability of this scheme is addressed through an extension of the master stability function technique to include adaptation. We observe that noise and/or slight nonidenticality between the coupled systems can be responsible for the occurrence of intermittent bursts of large desynchronization events (bubbling). Moreover, our numerical computations show that, for our adaptive synchronization scheme, the parameter space region corresponding to bubbling can be rather substantial. This observation becomes important to experimental realizations of adaptive synchronization, in which small mismatches in the parameters and noise cannot be avoided. We also find that, for our coupled systems with adaptation, bubbling can be caused by a slow drift in the coupling strength.

## I Introduction

It has been shown Fujisaka and Yamada (1983); Afraimovich et al. (1986); Pecora and Carroll (1998) that, in spite of their random-like behavior, the states of a collection of interacting chaotic systems that are identical can synchronize (i.e., be attracted toward a common chaotic evolution, ) provided that they are properly coupled. This phenomenon has been the basis for proposals for secure communication Cuomo and Oppenheim (1993); Argyris et al. (2008); Feki (2003), system identification Abarbanel et al. (2008); Creveling et al. (2008); Quinn et al. (2009); Sorrentino and Ott (2009a), data assimilation So et al. (1994); Duane et al. (2006), sensors Sorrentino and Ott (2009b), information encoding and transmission Hayes et al. (1993); Dronov et al. (2004), multiplexing Tsimring and Sushchik (1996), combatting channel distortion Sharma and Ott (2000), etc. In all of these applications it is typically assumed that one has accurate knowledge of the interaction between the systems, allowing one to choose the appropriate coupling protocol at each node (here we use the network terminology, referring to the chaotic systems as nodes of a connected network whose links correspond to the input that node receives from node ). In a recent paper Sorrentino and Ott (2008), an adaptive strategy was proposed for maintaining synchronization between identical coupled chaotic dynamical systems in the presence of a priori unknown, slowly time varying coupling strengths (e.g., as might arise from temporal drift of environmental parameters). This strategy was successfully tested on computer simulated networks of many coupled dynamical systems in which, at each time, every node receives only one aggregate signal representing the superposition of signals transmitted to it from the other network nodes. In addition, the strategy has also been successfully implemented in an experiment on coupled optoelectronic feedback loops Ravoori et al. (2009). Furthermore, a more generalized adaptive strategy, suitable for sensor applications, has also been proposed Sorrentino and Ott (2009b).

In past works, various other schemes for adaptive synchronization of chaos have also been proposed Mossayebi et al. (1991); Johnson Jr. and Thorp (1994); Parlitz and Kocarev (1996); Wu et al. (1996); Chua et al. (1996); di Bernardo (1996); Sobiski and Thorp (1998); Liao (1998); Zhou and Kurths (2006); De Lellis et al. (2008). So far, in all these studies, when the question of stability of the considered adaptive schemes has been studied, the question has been addressed using the Lyapunov function method (see e.g., Parlitz and Kocarev (1996); Wu et al. (1996); di Bernardo (1996); Liao (1998)), which provides a sufficient but not necessary condition for stability. While this technique has the advantage that it can sometime yield global stability conditions, it also has the disadvantages that its applicability is limited to special cases, and its implementation, when possible, requires nontrivial system specific analysis. In this paper, we address the stability of adaptive synchronization for the example of the scheme discussed in Ref. Sorrentino and Ott (2008). In particular, our analysis will extend the previously developed stability analysis of chaos synchronization by the master stability function technique Fujisaka and Yamada (1983); Pecora and Carroll (1998) to include adaptation. We will observe that the range in which the network eigenvalues are associated with stability, is dependent on the choice of the parameters of the adaptive strategy. The type of analysis we present, while for a specific illustrative adaptive scheme, can be readily applied to other adaptive schemes (e.g., those in Zhou and Kurths (2006); De Lellis et al. (2008)).

As compared to the Lyapunov technique, master stability techniques are much more generally applicable but they provide conditions for local, rather than global stability. We also note that, within that context, the master stability technique allows one to distinguish between stability of typical chaotic orbits and stability of atypical orbits within the synchronizing chaotic attractor [i.e., stability to ‘bubbling’ Ashwin et al. (1994a, b); Venkataramani et al. (1996); Hunt and Ott (1996); Restrepo et al. (2004); Ott (2002); see Secs. III and IV].

In Sec. II we review the adaptive synchronization strategy formulation of Ref. Sorrentino and Ott (2008), which applies to a network of chaotic systems with unknown temporal drifts of the couplings. In Sec. III, we present a master stability function approach to study linear stability of the synchronized solution in the presence of adaptation; we also consider a generalized formulation of our adaptive strategy and study its stability. Numerical simulations are finally presented in Sec. IV. Our work in Sec. IV highlights the important effect of bubbling in the dynamics.

## Ii Adaptive strategy formulation

As our example of the application of the master stability technique to an adaptive scheme, we consider the particular scheme presented in Ref. Sorrentino and Ott (2008). To provide background, in this section we present a brief exposition of a formulation similar to that in Ref. Sorrentino and Ott (2008), as motivated by the situation where the couplings are unknown and drift with time. We consider a situation where the dynamics at each of the network nodes is described by,

(1) |

where, is the -dimensional state of system ; determines the dynamics of an uncoupled () system (hereafter assumed chaotic), ; is a scalar output function, . We take to be a constant -vector, , with , and the scalar is a constant characterizing the strength of the coupling. The scalar signal each node receives from the other nodes in the network is,

(2) |

The quantity is an adjacency matrix whose value specifies the strength of the coupling from node to node . We note that if

(3) |

then Eq. (1) admits a synchronized solution,

(4) |

where satisfies

(5) |

which corresponds to the dynamics of an isolated system. We regard the as unknown at each node , while the only external information available at node is its received signal (2). The goal of the adaptive strategy is to adjust so as to maintain synchronism in the presence of slow, a priori unknown time variations of the quantities . That is, we wish to maintain approximate satisfaction of Eq. (3). For this purpose, as discussed in Ref. Sorrentino and Ott (2008), our scheme can be extended to the case where the output function is -dimensional, , where and is an dimensional matrix. For simplicity we consider . We assume that each node independently implements an adaptive strategy. At each system node , we define the exponentially weighted synchronization error , where

(6) |

and we evolve so as to minimize this error (a slightly more general approach is taken in Sorrentino and Ott (2008)). Hence we set equal to zero to obtain,

(7) |

By virtue of , we obtain the numerator and the denominator on the right hand side of Eq. (7) by solving the differential equations,

(8a) | |||

(8b) |

Since the dynamics of is imagined to occur on a timescale which is slow compared to the other dynamics in the network, we can approximate as constant . This essentially assumes that we are dealing with perturbations from synchronization whose growth rates (in the case of unstable synchronization) or damping rates (in the case of stable synchronization) have magnitudes that substantially exceed . Under this assumption, we note that Eqs. (1), (7), and (8) admit a synchronized solution, given by Eqs. (4), (5), and

(9a) | |||

(9b) |

To simplify the notation, in what follows, we take , , and ; e.g., we can now write,

(10) |

where . If the synchronization scheme is locally stable, we expect that the synchronized solution (4),(5), and (9) will be maintained under slow time evolution of the couplings .

## Iii Stability analysis

### iii.1 Linearization and master stability function

Our goal is to study the stability of the reference solution (4),(5), and (9). By linearizing Eqs. (1) and (8) about (5), and (9), we obtain,

(11a) | |||

(11b) |

where we have introduced the new variable .

Equations (11) constitute a system of coupled equations. In order to simplify the analysis, we seek to decouple this system into independent systems, each of dimension . For this purpose we seek a solution where is in the form , where is a time independent scalar that depends on and is a -vector that depends on time but not on . Substituting in Eqs. (11a),(11b), we obtain,

(12a) | |||

(12b) |

To make Eqs. (12) independent of , we consider and , where is a quantity independent of . Namely, the possible values of are the eigenvalues, , corresponding to linearly independent eigenvectors , where . This gives,

(13a) | |||

(13b) |

which is independent of , but depends on the eigenvalue . Considering the typical case where there are distinct eigenvalues of the matrix , we see that Eqs. (12) constitute decoupled linear ordinary differential equations for the synchronization perturbation variables and . All the rows of sum to . Therefore has at least one eigenvalue , corresponding to the eigenvector . Furthermore, since for all , we have by the Perron-Frobenius theorem that , and thus . For , Eq. (13a) becomes,

(14) |

This equation reflects the chaos of the reference synchronized state (Eq. (5)) and (because all the are equal) is associated with perturbations which are tangent to the synchronization manifold and are therefore irrelevant in determining synchronization stability. Stability of the synchronized state thus demands that Eqs. (12) yield exponential decay of and for all the eigenvalues , excluding this eigenvalue.

Then it becomes possible to introduce a master stability function Fujisaka and Yamada (1983); Pecora and Carroll (1998), , that associates the maximum Lyapunov exponent of system (13) with . In so doing, one decouples the effects of the network topology (reflected in the eigenvalues and hence the relevant values of ) from the choices of . In general an eigenvalue, and hence also , can be complex. For simplicity, in our discussion and numerical examples to follow, we assume that the eigenvalues are real (which is for instance the case when the adjacency matrix is symmetric). For any given value of stability demands that for all those values of corresponding to the eigenvalues .

Following Refs. Barahona and Pecora (2002); Nishikawa et al. (2003); Chavez et al. (2005); Hwang et al. (2005), we now introduce the following definition of synchronizability. Let us assume that the master stability function is negative in a bounded interval of values of , say . Then, in order for the network to synchronize, two conditions need to be satisfied, (i) , and (ii) , where () is the smallest (largest) network eigenvalue over all the eigenvalues . The network synchronizability is defined as the width of the range of values of , for which . Assuming that and are assigned (e.g., the network topology is given), then the network synchronizability increases with the ratio . In what follows, we will compare different adaptive strategies in terms of their effects on the synchronizability ratio .

In our analysis above, since we divide by , we have implicitly assumed that all the , i.e., that every node has an input. There is, however, a case of interest where this is not so, and this case requires separate consideration. In particular, say there is one and only one special node (which we refer to as the maestro or sender) that has no inputs, but sends its output to other nodes (which interact with each other), and we give this special node the label . Since node receives no inputs, we do not include adaption on this node, and we replace Eq.(1) for by . In addition, when investigating the stability of the synchronized state, it suffices to set (i.e., not to perturb the maestro). Following the steps of our previous stability analysis, we again obtain Eqs. (12) and (13), but with important differences. Namely, Eqs. (12) now apply for , the values of in Eqs. (13) are now the eigenvalues of the matrix for ; i.e., only the interactions between the nodes are included in this matrix. Note that is still given by , still including the input from the maestro node. Also since , all of the eigenvalues represent transverse perturbations and are therefore relevant to stability. (This is in contrast to the case without a maestro in which we had to exclude an eigenvalue, i.e., corresponding to . For a similar discussion for the case of the standard master stability problem with no adaptation, see Sorrentino et al. (2007).) The simplest case of this type (used in some of our subsequent numerical experiments) is the case , where there is one receiver node () and one sender/maestro node (). Since there is only one receiver node whose only input is received by the sender, reduces to the scalar and , yielding .

As stated above, in (4) is an orbit of the uncoupled system (4). In general, two types of orbits are of interest: (i) a typical chaotic orbit on the relevant chaotic attractor of (4), and (ii) the orbit that is ergodic on the maximally synchronization-unstable invariant subset embedded within the relevant chaotic attractor of (4). Here, by ‘relevant chaotic attractor’ we mean that, if the system (4) has more than one attractor, then we restrict attention to that attractor on which synchronized motion is of interest. Also, in (i), by the word ‘typical’, we mean orbits of (4) that ergodicly generate the measure that applies for Lebesgue almost every initial condition in the attractor’s basin of attraction. In this sense, the orbit in (ii) is not typical. In general the criterion for stability as assessed by (ii) is more restrictive than that assessed by (i). Conditions in which the synchronized dynamics is stable according to (i), but unstable according to (ii), are referred to as the ‘bubbling’ regime Ashwin et al. (1994a, b); Venkataramani et al. (1996); Hunt and Ott (1996); Ott (2002). In previous work on synchronization of chaos Ashwin et al. (1994a, b); Venkataramani et al. (1996); Hunt and Ott (1996); Restrepo et al. (2004); Ott (2002), it has been shown that, when the system is in the bubbling regime, small noise and/or small ‘mismatch’ between the coupled systems can lead to rare, intermittent, large deviations from synchronism, called ‘desynchronization bursts’ ^{1}^{1}1In addition, to desynchronism bursts, it is also possible that a large desynchronization orbit event can result in capture of the system orbit on a desynchronized attractor with a so-called riddled basin of attraction (e.g., see Refs. Ashwin et al. (1994a, b); Venkataramani et al. (1996); Ott and Sommerer (1994); Ott (2002)). This possibility, although a typical occurrence in such situations, did not manifest itself in the particular system used for our numerical studies in Sec. IV. Thus we, henceforth, restrict our discussion to the case that bubbling is associated with bursts.. By small system mismatch we mean that, for each node , the functions in (1) are actually different, , but that these differences are small (i.e., is small, where now denotes a reference uncoupled system dynamics; e.g., averaged over ). With reference to our adaptive synchronization problem (1), we shall see that, in addition to small noise and small mismatch in , bursting can also be induced by slow drift in the unknown couplings . From the practical, numerical perspective, the complete and rigorous application of the stability criterion (ii) is impossible, since there will typically be an infinite number of distinct invariant sets embedded in a chaotic attractor, and, to truly be sure of stability, each of these must be found and numerically tested. In practice, therefore, as done previously by others, we will evaluate stability for all the unstable periodic orbits embedded in the attractor up to some specified period. This will give a necessary condition for stability according to (ii), and furthermore, it has been argued and numerically verified in Ref. Hunt and Ott (1996) that stability, as assessed from a large collection of low period periodic orbits (and embedded unstable fixed points, if they exist in the relevant attractor), will extremely often yield the true delineation of the parameters of the bubbling regime, or, if not, an accurate approximation of it. Our numerical results of Sec. IV lend further support to this idea.

### iii.2 Generalized adaptive strategy

We now analyze a generalization of our adaptive strategy. Namely, we replace Eq. (8b) by,

(15) |

where is an arbitrary function of , normalized so that . The key point is that at synchronism , corresponding to ; and thus, since we take , the synchronized solution is unchanged. The stability analysis for this generalization is given in the Appendix I and results in the following master stability equations,

(16a) | |||

(16b) |

where , and denotes evaluated at . We then introduce a master stability function , that associates the maximum Lyapunov exponent of system (16) with and .

Thus we expect that, when our modified adaptive scheme is stable, it will again relax to the desired synchronous solution. The difference between the stability of the modified scheme (Eqs. (13)) and the stability of the original scheme (corresponding to Eqs. (16) with ), is that, by allowing the freedom to choose the value of , we can alter the stability properties of the synchronous state. We anticipate that, by properly adjusting , we may be able to tailor the stability range to better suit a given situation.

In the case of , Eq. (16b) reduces to,

(17) |

which has a Lyapunov exponent , where is the time average of , . For , Eq. (17) implies that decays to zero. Thus, if we choose a large enough value of , stability of the synchronized state is determined by (16a) with set equal to zero, and Eq. (16a) reduces to the master stability function for the determination of the stability of the system without adaptation Pecora and Carroll (1998). Therefore, in the case of , , the stable range of is independent of and is the same as that obtained for the case in which adaptation is not implemented ().

## Iv Numerical experiments

In our numerical experiments we consider the example of the following Rössler equation, for which, , ,

(18) |

with the parameters , and , and we use , and . In Fig. 1 the master stability functions calculated from Eq. (13) for the adaption scheme of Sec. II are plotted for three different values of , i.e., (dashed, dashed/dotted, and dotted curves, respectively). In addition, for comparison, we also plot the result of computations for the case in which no adaptation is introduced, corresponding to the reduced system (solid curves). The master stability function is shown in black (respectively, grey) for the cases that is a typical chaotic orbit in the attractor (respectively, the maximally unstable periodic orbit embedded in the attractor for periodic orbits of period up to four surface of section piercings; see Appendix II for a brief account of how the unstable periodic orbits were obtained). We say that synchronization is ‘high quality’ stable in the range of for which for all orbits (i.e., including the periodic orbits) is negative. As can be seen, by changing the parameter , the -range of stability can be dramatically modified. The bubbling range is given by the values of for which for a typical orbit but for the maximally unstable periodic orbit embedded in the attractor.

Figure 2 is a level curve plot of the values assumed by the master stability function evaluated for being a typical chaotic orbit. In the figure, the area of stability (corresponding to ) is delimited by the thick -level contour line. From the figure, we see that the width of the range of stability increases with . In Figs. 3(a,b) a comparison between the areas of stability is given, for the cases in which is a typical chaotic orbit in the attractor, and for the case that is the maximally unstable periodic orbit embedded in the attractor of period up to four. The thick solid (respectively, dashed) curves bound the area in which the master stability function is negative for corresponding to a typical chaotic orbit in the Rösller attractor (respectively, for corresponding to the maximally unstable periodic orbit embedded in the attractor of period up to four). The bubbling area falls between the dashed and the continuous contour lines.

Interestingly, we see that for , high-quality stability can never be achieved for any , while, in contrast, stability with respect to typical chaotic orbits (i.e., with bubbling) is achievable. Let denote the upper and lower values of at the borders of the stability regions with respect to a typical (t) chaotic orbit and with respect to unstable periodic orbits (p) in the synchronizing attractor. E.g., high-quality synchronism applies for and the bubbling regime corresponds to or . In terms of these quantities, useful measures for assessing the possibility of achieving stable synchronism for a given network topology are the ‘synchronizability’ ratios Barahona and Pecora (2002); Nishikawa et al. (2003); Chavez et al. (2005); Hwang et al. (2005),

(19) |

In what follows, where convenient, we drop the subscripts and with the understanding that the discussion may be taken to apply to stability based on either typical or periodic orbits. Noting that synchronism is stable for , and that , we consider the coupling network topology-dependent ratio where () denotes the maximum (minimum) eigenvalue of the adjacency matrix (not including the eigenvalue corresponding to the eigenvector ). Recall that . Since for stability, if

(20) |

then the system can be made stable by adjustment of the constant , but, if , then it is impossible to choose a value of for which for all the relevant eigenvalues , and stability is unachievable. Figure 3(c) shows plots of and versus for the same parameters as used in Figs. 3(a,b). Note that, for these computations, the values of without adaption (i.e., and ) always exceed the corresponding values with adaption. We have also found this to be true for the generalized adaptive scheme of Sec. III B (which includes the additional adaption parameter ). However, we do not know whether this is general, or is limited to our particular example (Eq. (18) with , and our choices of the parameters ,, and ).

To test our linear results in Fig. 3, we have also performed fully nonlinear numerical simulations for a simple network consisting of a sender system (labeled 1) connected to a receiver (labeled 2). In this case Eq. (1) becomes

(21a) | ||||

(21b) |

and is a scalar.

Each data point shown in Figs. 3(a,b) corresponds to a run, where the sender was given a random initial condition and random values for and were chosen in the plotted range. After waiting sufficient time to ensure that the sender state is essentially on the attractor, the -variable of the receiver state was initialized by a displacement of from the -variable of the sender state. A step-size of was used for a run time of time units, over which we recorded the normalized synchronization error,

(22) |

where indicates a time average and the subscript denotes evolution on the synchronous state (i.e., using dynamics from Eq. (4)). If, in that time span, never converged to and, at some point, exceeded , the run was considered to be unstable (corresponding to an in the figure). If converged to , a mismatch in the Rössler parameter was introduced to the receiver, and the run of duration time units was repeated with an initial separation of . If, at any time during the run, ever exceeded , the run was considered to be bubbling (corresponding to a green circle in the figure), otherwise the run was considered to be stable (corresponding to a red triangle in the figure). We see that the master stability computations of the high-quality stable, bubbling, and unstable regions (the solid and dashed lines) correspond well with these results. We also did a sampling of points up to period 5 and did not find that this altered our results. From Fig. 3(a), we observe the presence of a few green circles (i.e, bubbling) within the high-quality synchronization area, delimited by the dashed line. In reference to this observation, we note that (i) for the case in which a small parameter mismatch is present, the synchronization error is expected to vary smoothly with parameter variation, and there is no sharp transition from the stable to the bubbling regime; and (ii) our computations show that close to the dashed line, the master stability function associated with the most unstable invariant set embedded in the attractor is rather small. Facts (i) and (ii) explain our difficulty in using our nonlinear computations to clearly separate the bubbling from the stable regions about the dashed line in Fig. 3(a). An important point concerning Figs. 3(a,b) is that the area associated with bubbling in Fig. 3(a) is rather substantial. This observation would become particularly important in experimental realizations of adaptive synchronization, since small mismatches in the parameters and noise cannot be avoided in experiments.

Figure 4 shows a sample plot of the normalized synchronization error , versus . We implemented our adaptive strategy with values of and , corresponding to the bubbling regime (see Fig. 3), , and the receiver has mismatch in the parameter . The two insets are zooms showing phase-space projections in the plane , over two different time intervals. Inset (b) corresponds to a range of time between bursts () and shows that during this time the orbit is essentially that of a typical chaotic orbit. Inset (a) shows the orbit trajectory for a range of time during which a burst is growing. It is seen from inset (b) that during the time range of the growing burst the orbit closely follows a period 4 orbit embedded in the attractor. The burst is evidently caused by the instability of this period 4 orbit to perturbations that are transverse to the synchronization manifold.

We have also performed numerical master stability computations for our generalized adaptive strategy, presented in Sec. IIIb. This is shown in Fig. 5, where the and curves, corresponding respectively to the largest (smallest) values of for which (), are plotted versus for three different values of for typical chaotic orbits. For small (e.g., in the figure), the range of stability is almost independent of , while for larger values of the choice of can significantly affect the -range of stability. As expected, at , and are independent of .

Finally, we investigated whether, for our coupled systems with adaptation, bubbling can be caused by a slow drift in the coupling strength. For this purpose we now take the parameter in Eq. (21) to have a slow time drift,

(23) |

We implemented our adaptive strategy with values of and , corresponding to the bubbling regime (see Fig. 3). For most of the time there is good synchronization between the sender and the receiver, but we also observed the intermittent occurrence of short, intense desynchronization bursts. Figure 6 shows the synchronization error versus . Note that in the absence of parameter drift ( constant), the synchronization error would eventually become zero. This simulation shows that, similarly to the previously reported burst-inducing effect of small parameter mismatch or noise, drift also promotes the continuous intermittent occurrence of bursting.

## V Conclusion

This paper is concerned with the study of stability of adaptive synchronization of chaos in coupled complex networks (e.g., sensor networks). As an example addressing this issue, we consider a recently proposed adaptive scheme for maintaining synchronization in the presence of a priori unknown slow temporal drift in the couplings Sorrentino and Ott (2008). In contrast with previous approaches (e.g., Parlitz and Kocarev (1996); Wu et al. (1996); di Bernardo (1996); Liao (1998)), based on system specific use of the Lyapunov function technique, we present a master stability analysis which predicts the exact ranges of stability for the synchronized state. We observe that the stable range of synchronism can be sensitively dependent on the adaption parameters. Moreover, we are able to predict the onset of bubbling, which occurs when the synchronized state is stable for typical chaotic orbits but is unstable for certain unstable periodic orbits within the synchronized chaotic attractor. We define stability to be high quality when the synchronized state is stable with respect to all the orbits embedded in the attractor and numerically find the regions of ‘high quality stability’ for a given system of interest. We also found that, for our coupled systems with adaptation, bubbling can be caused by a slow drift in the coupling strength, in addition to small noise and small mismatch in . We emphasize that, since parameter mismatch, noise and drift are ubiquitous in experimental situations, and since (e.g., Fig. 3(a)) bubbling can occupy substantial regions of parameter space, consideration of bubbling can be expected to be essential for determining the practical feasibility of chaos synchronization applications.

We thank Anurag. V. Setty and Bhargava Ravoori for the enlightening discussions.

This work was supported by ONR grant N00014-07-1-0734.

## APPENDIX I: Stability of the generalized adaptive strategy

We note that the function in Eq. (15), when evaluated about (9), is equal to one. Then, by linearizing Eqs. (1), (8a), and (15) about (9), we obtain,

(24a) | |||

(24b) |

As in our derivation of Eqs. (13), we again set , where is a constant scalar that depends on and is a vector that depends on time but not on . Equations (24), then become

(25a) | |||

(25b) |

To make Eqs. (25) independent of , we again consider and take to be the eigenvalues of , resulting in Eqs. (16).

## APPENDIX II: Determination of unstable periodic orbits

To account for the phenomenon of bubbling, it is necessary to look not just at typical (that is, chaotic) orbits of the uncoupled oscillator, but the periodic orbits embedded in the chaotic attractor as well. As there are a (countably) infinite number of such orbits, it is impossible to account for them all. However, as shown by Hunt and Ott, the optimal periodic orbits of maximal transverse instability tend to be those of low period Hunt and Ott (1996). Thus, for our analysis, it was found to be sufficient to consider only those orbits with a period less than some appropriately chosen limit.

To find these low-period orbits for the Rössler attractor, we initialized an uncoupled oscillator with random initial conditions, waited for it to settle onto the attractor, then recorded its orbits for some suitable length of time at high temporal precision. We then noted each piercing of the surface of section in the positive- direction (). To a high degree of approximation, the coordinates of these points were found to lie a curve, thus suggesting that it is possible to reduce the three-dimensional flow to a one-dimensional map. We then plotted vs. ; that is, the coordinate of the piercing versus the coordinate of the . Each intersection of this curve with the line represents the coordinate of an initial condition for an orbit that starts on the surface of section and returns to its original position after piercing of the surface of section. With two coordinates (namely and ) known, all that remains is to find the value of such that lies on the attractor.

Of course, for , many of these intersections will be redundant, as every period orbit pierces the surface of section times, thus producing intersections on the curve. In addition, each curve will have intersections corresponding to orbits of any period that is a factor of . As an example, consider the curve vs. . The Rössler system used in this paper has three Period 4 orbits, one Period 2 orbit and one Period 1 orbit. Thus, the number of times vs. will intersect is .

As these orbits are inherently unstable, error accumulated through numerical integration can result in a trajectory leaving the periodic orbit after only a small number of periods. Thus, for the long term computation of Lyapunov exponents to obtain the master stability function, it is advisable to compute the trajectory for only a single period, then return the oscillator to its initial position , and repeat as often as needed.

## References

- Fujisaka and Yamada (1983) H. Fujisaka and T. Yamada, Prog. Theor. Phys. 69, 32 (1983).
- Afraimovich et al. (1986) V. S. Afraimovich, N. N. Verichev, and M. I. Rabinovich, Inv. VUZ Radiofiz. 29, 795 (1986).
- Pecora and Carroll (1998) L. Pecora and T. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
- Cuomo and Oppenheim (1993) K. M. Cuomo and A. V. Oppenheim, Phys. Rev. Lett. 71, 177 (1993).
- Argyris et al. (2008) A. Argyris, M. Hamacher, K. E. Chlouverakis, A. Bogris, and D. Syvridis, Phys. Rev. Lett. 100, 194101 (2008).
- Feki (2003) M. Feki, Chaos, Solitons and Fractals 18, 141 (2003).
- Abarbanel et al. (2008) H. D. I. Abarbanel, D. R. Creveling, and J. M. Jeanne, Phys. Rev. E 77, 016208 (2008).
- Creveling et al. (2008) D. R. Creveling, P. E. Gill, and H. D. I. Abarbanel, Phys. Lett. A 372, 2640 (2008).
- Quinn et al. (2009) J. C. Quinn, P. H. Bryant, D. R. Creveling, S. R. Klein, and H. D. I. Abarbanel, Phys. Rev. E 80, 016201 (2009).
- Sorrentino and Ott (2009a) F. Sorrentino and E. Ott, Chaos 19, 033108 (2009a).
- So et al. (1994) P. So, E. Ott, and W. P. Dayawansa, Phys. Rev. E 49, 2650 (1994).
- Duane et al. (2006) G. S. Duane, J. J. Tribbia, and J. B. Weiss, Nonlinear processes in Geophysics 13, 601 (2006).
- Sorrentino and Ott (2009b) F. Sorrentino and E. Ott, Phys. Rev. E 79, 016201 (2009b).
- Hayes et al. (1993) S. Hayes, C. Grebogi, and E. Ott, Phys. Rev. Lett. 70, 3031 (1993).
- Dronov et al. (2004) V. Dronov, M. Hendrey, T. M. Antonsen, and E. Ott, Chaos 14, 30 (2004).
- Tsimring and Sushchik (1996) L. S. Tsimring and M. M. Sushchik, Phys. Lett. A 213, 155 (1996).
- Sharma and Ott (2000) N. Sharma and E. Ott, Int. J. Bif. Chaos 10, 777 (2000).
- Sorrentino and Ott (2008) F. Sorrentino and E. Ott, Phys. Rev. Lett. 100, 114101 (2008).
- Ravoori et al. (2009) B. Ravoori, A. B. Cohen, A. V. Setty, F. Sorrentino, T. E. Murphy, E. Ott, and R. Roy, arXiv:0907.3894 (2009).
- Mossayebi et al. (1991) F. Mossayebi, H. K. Qammar, and T. T. Hartley, Phys. Lett. A 161, 255 (1991).
- Johnson Jr. and Thorp (1994) C. R. Johnson Jr. and J. S. Thorp, IEEE Sig. Proc. Lett. 1, 194 (1994).
- Parlitz and Kocarev (1996) U. Parlitz and L. Kocarev, Int. J. Bifurcation Chaos 6, 581 (1996).
- Wu et al. (1996) C. W. Wu, T. Yang, and L. O. Chua, Int. J. Bifurcation Chaos 6, 455 (1996).
- Chua et al. (1996) L. O. Chua, T. Yang, G.-Q. Zhong, and C. W. Wu, Int. J. Bifurcation Chaos 6, 189 (1996).
- di Bernardo (1996) M. di Bernardo, Int. J. Bifurcation Chaos 6, 557 (1996).
- Sobiski and Thorp (1998) D. J. Sobiski and J. S. Thorp, IEEE Trans. Circ. Syst. I 45, 194 (1998).
- Liao (1998) T.-L. Liao, Chaos Solitons Fractals 9, 1555 (1998).
- Zhou and Kurths (2006) C. Zhou and J. Kurths, Phys. Rev. Lett. 96, 164102 (2006).
- De Lellis et al. (2008) P. De Lellis, M. di Bernardo, and F. Garofalo, Chaos 18, 037110 (2008).
- Ashwin et al. (1994a) P. Ashwin, J. Buescu, and I. N. Stewart, Phys. Lett. A 193, 126 (1994a).
- Ashwin et al. (1994b) P. Ashwin, J. Buescu, and I. N. Stewart, Nonlinearity 9, 703 (1994b).
- Venkataramani et al. (1996) S. C. Venkataramani, B. R. Hunt, and E. Ott, Phys. Rev. E 54, 1346 (1996).
- Hunt and Ott (1996) B. R. Hunt and E. Ott, Phys. Rev. Lett. 76, 2254 (1996).
- Restrepo et al. (2004) J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 69, 066215 (2004).
- Ott (2002) E. Ott, Chaos in Dynamical Systems (Cambridge University Press, second edition, 2002).
- Barahona and Pecora (2002) M. Barahona and L. Pecora, Phys. Rev. Lett. 89, 054101 (2002).
- Nishikawa et al. (2003) T. Nishikawa, A. Motter, Y. Lai, and F. Hoppensteadt, Phys. Rev. Lett. 91, 014101 (2003).
- Chavez et al. (2005) M. Chavez, D. Huang, A. Amann, H. Hentschel, and S. Boccaletti, Phys. Rev. Lett. 94, 218701 (2005).
- Hwang et al. (2005) D. Hwang, M. Chavez, A. Amann, and S. Boccaletti, Phys. Rev. Lett. 94, 138701 (2005).
- Sorrentino et al. (2007) F. Sorrentino, M. di Bernardo, F. Garofalo, and G. Chen, Phys. Rev. E 75, 046103 (2007).
- Ott and Sommerer (1994) E. Ott and J. C. Sommerer, Phys. Lett. A 188, 39 (1994).