Complex behavior of chaotic synchronization under dual coupling channels

Most previous works on complete synchronization of chaotic oscillators focused on the one-channel interaction scheme where the oscillators are coupled through only one variable or a symmetric set of variables. Using the standard framework of master-stability function (MSF), we investigate the emergence of complex synchronization behaviors under all possible configurations of two-channel coupling, which include, for example, all possible cross coupling schemes among the dynamical variables. Utilizing the classic Rössler and Lorenz oscillators, we find a rich variety of synchronization phenomena not present in any previously extensively studied, single-channel coupling configurations. For example, in many cases two coupling channels can enhance or even generate synchronization where there is only weak or no synchronization under only one coupling channel, which has been verified in a coupled neuron system. There are also cases where the oscillators are originally synchronized under one coupling channel, but an additional synchronizable coupling channel can, however, destroy synchronization. Direct numerical simulations of actual synchronization dynamics verify the MSF-based predictions. Our extensive computation and heuristic analysis provide an atlas for synchronization of chaotic oscillators coupled through two channels, which can be used as a systematic reference to facilitate further research in this area.


Introduction
Synchronization has been widely observed in many natural, social and technological systems, and has attracted much attention during the last two decades [1][2][3][4][5][6][7][8][9][10][11][12]. In previous works, complete synchronization of identical chaotic oscillators has been a special focus due to the availability of analytic framework of the master-stability function (MSF) in determining the critical coupling strength at which synchronization arises and, more generally, the regions in the parameter space where stable synchronization may emerge [7,13,14]. In general, the MSF reduces the synchronization analysis to the interplay between two aspects: the dynamics of a single oscillator and the coupling topology characterized by the eigenvalues of the coupling matrix. With knowledge of the dynamics of the individual nodes and the nature of the pairwise coupling to infer the MSF, for any coupled system one can thus predict whether synchronization can occur based on the eigenvalues of the coupling matrix and the coupling strength, without the need of actual simulations. For example, a common class of MSFs have the property that, in the plot of the MSF versus some generalized coupling parameter, there exists a finite interval in which the MSF assumes negative values. For the normalized coupling strength either below or above this interval, the system cannot be physically synchronized. In this case, synchronization is determined by the largest and the smallest nontrivial eigenvalues of the coupling matrix, and their ratio is called eigenratio, which has been used extensively in the study of the synchronizability of complex networks of various types of topology [15][16][17][18][19][20]. A general observation is that a smaller eigenratio is beneficial to synchronization and thus is more desirable in designing synchronous networks [16]. The effect of different coupling schemes among the dynamical variables in a large number of chaotic oscillators was investigated [14], confirming the generality of the MSF property so discussed.
Most previous works in the field of complete synchronization of nonlinear and/or complex systems assumed one-channel coupling or some coupling scheme with certain symmetry. For example, for a network system of three-dimensional oscillators, coupling between any pair of interacting oscillators was usually assumed to exist between one pair of dynamical variables, one from each oscillator. In real-world systems, however, the dynamical interactions among the coupled oscillators can be more complicated [21]. In particular, the dynamics of each individual oscillator in the network might be affected simultaneously by several different dynamical variables from other oscillators. In this regard, two delayed signals [22] and multiple delayed signals [23] have been introduced to stabilize the fixed points of various chaotic dynamical systems. Multi-path propagation with multiple delays was proposed to improve the efficiency in network-controlled systems [24]. Conjugate coupling was introduced to regulate the oscillation death phenomenon [25,26] and, it was also found that multi-channel coupling can provide a much larger domain for oscillation death than the case of single-channel coupling [27]. In spite of these investigations, a systematic study of the effect of multi-channel coupling on synchronization has been lacking.
In this paper, we study synchronization of chaotic oscillators under two coupling channels. To be concrete but without loss of generality, we employ the classic Rössler and Lorenz chaotic oscillators and calculate their MSFs for all possible combinations of two coupling channels. We find three general types of synchronous dynamics with respect to the effect of more than one coupling channel: (1) enhancement of synchronization, (2) induction of synchronization, and (3) deterioration or even destruction of synchronization. In particular, for case (1), the MSF of the system is negative under one coupling channel, and the addition of another coupling channel with the same type can make the MSF even more negative. For case (2), the system cannot be synchronized under one coupling channel, but synchronization can emerge under two coupling channels. For case (3), the system can be synchronized under one coupling channel, but adding another coupling channel can counter-intuitively destroy synchronization. Our extensive analysis and computation provide a comprehensive picture of the complex synchronization behaviors that can arise under multiple coupling channels, and can serve as a systematic reference to facilitate study of synchronization in complex nonlinear dynamical systems.
In section 2, we briefly review the MSF formalism. In section 3, we provide a systematic analysis of the MSFs of the Rössler and Lorenz oscillators under two coupling channels and demonstrate an abnormal enhancement of synchronizable region for a coupled HR neuron system. Concluding remarks are provided in section 4.

The MSF formalism
The MSF formalism provides a general framework for determining the synchronizability of a system of coupled nonlinear oscillators [7]. Say the intrinsic dynamic of each isolated oscillator can be described by where x i is a d-dimensional vector of dynamical variables. The coupled-oscillator system with two coupling channels can be written as where N is the number of oscillators in the system, H x ( ) 1 and H x ( ) 2 are two distinct coupling functions, each characterizing one coupling channel, ɛ 1 and ɛ 2 are the corresponding coupling strength.
Here, for convenience, we write the coupling functions H 1 and H 2 separately to have independent control of the coupling strength. The coupling matrix G is determined by the undirected (weighted) connections among the oscillators and assumed to have real eigenvalues, and the sum over rows (∑ = G j N ij 1 ) is equal to zero, which guarantees that the synchronous state, , is a solution of the coupled system (2). Here we adopt the convention that the diagonal elements of G are positive while the non-diagonal elements are negative. The variational equations governing the time evolution of the set of infinitesimal deviations from the synchronous solution, δ where DF s ( ), DH s ( ) 1 and DH s ( ) 2 are the d × d Jacobian matrices of the corresponding vector functions evaluated at t s( ). The coupling matrix G can be diagonalized: are the set of eigenvalues satisfying μ μ μ = < ⩽ ⩽ 0 ... , and the columns of Q are the set of corresponding eigenvectors. The transform ξ δ = − Q x · 1 leads to the following block-diagonally decoupled form of equation (3): i N 2, , ) be the specific set of values of the normalized coupling parameters K 1 and K 2 , we see that each block of the above decoupled equation is structurally identical but differs only in K 1i and K 2i . This leads to the following generic form for all the decoupled blocks: The largest Lyapunov exponent of equation (4) is the MSFΨ K K ( , ) 1 2 [7], which is also the largest transverse Lyapunov exponent of equation (2).
For an oscillator network with two coupling channels, a necessary condition for synchronization is that all the normalized coupling points K K ( , ) In this case, a small deviation from the synchronization state will diminish exponentially so that the synchronous manifold is stable. Once the coupling configuration is determined, the MSF can be calculated and the negative region of Ψ K K ( , ) and then consider the matrix equation with initial condition , where I is the identity matrix of order d [28]. This matrix equation is solved together with equation (1) that yields the trajectory t s( ). Both equations are integrated using the fourth-order Runge-Kutta method. Let ν t ( ) i ( = … i d 1, , ) be the ordered eigenvalues of t O( ). The Lyapunov exponents are given by The largest λ i is the MSFΨ K K ( , ) 1 2 . In our computation, 3 × 10 3 time units of t s( ) are first integrated to allow the system to settle into the attractor, then 3 × 10 5 time units are used to calculate the Lyapunov exponents. Time step is chosen to be δ = × − t 1 10 3 .

MSFs for two-channel couplings
The MSFs for typical low-dimensional nonlinear oscillators with a single coupling channel was investigated systematically in a previous work [14]. For a pair of d-dimensional oscillators, there are d × d possible coupling configurations. For two coupling channels, the d × d configurations can lead to × C d d 2 distinct coupling schemes (note that the order of K 1 and K 2 are interchangeable). Taking d = 3 as an example, we have = C 36 9 2 twochannel coupling configurations. We calculate the MSFs for all these configurations for the Rössler and Lorenz oscillators, as shown in figure 1. Typically for different coupling schemes, the values of MSFs can be orders of magnitude different. As a result, in order to identify regions of negative values of MSF, instead of Ψ, we plot its normalized value Ψ Ψ | | sign( ) 1 6 , as in figure 1. Note that Ψ Ψ | | sign( ) 1 6 has the same sign as Ψ, so the negative region will be the same.

Rössler system
We use the classic Rössler oscillator [29]: with α = 0.2 and γ = 9. For any two-channel coupling scheme, the Jacobian matrix DF will have −K 1 and −K 2 added to two of its elements. For example, for the → → 3 3 & 3 1coupling channel, DF is modified to The MSFs calculated from all 36 two-channel coupling configurations of the Rössler system are shown in figure 1. Most cases are not surprising in the sense that, if for each coupling channel the MSF is negative (or positive) so that the system can be synchronized (or not synchronized), having an additional coupling channel of the same type will lead to a stronger synchronous (or unsynchronized) state, as exemplified by the → → 2 2and1 1or the → → 3 3 and 1 2 cases. If one coupling channel leads to a negative value of MSF while the other channel leads to a positive value of MSF, the effect of combining both channels will depend on the strength of the couplings. For example, if the coupling leading to positive MSF is too strong, the MSF will become However, there are counterintuitive cases. For instance, consider a particular K K ( , ) . By each coupling channel alone, the system cannot be synchronized, but when the two coupling channels are simultaneously present, we can haveΨ < K K ( , ) 0 1 2 , i.e., synchronization can emerge. To be specific, let us consider the → → 3 3 and 3 1 coupling configuration as an example. We plot the two-channel MSF instead of the normalized MSF for this case and also the MSFs for the two corresponding single-channel coupling configurations in figure 2. We see that, for the → 3 1 single-channel coupling case, the trend of Ψ decreases (with small oscillations) with the normalized coupling strength K and crosses zero around K = 50, as shown in figure 2(c). For the → 3 3 case, Ψ remains positive and generally increases with K, as shown in figure 2(b). Thus, for the two-channel coupling scheme, one would not expect Ψ to be negative for < K 50 31 , regardless of the value of K 33 . Nonetheless, as shown in figure 2(a), there is an interesting interplay between the → 3 3 and the → 3 1 coupling configurations, which results in a relatively large region (the right side of the white solid line) of Ψ < 0 even for < K 50 31 . More specifically, assuming that K 31 is fixed at 20 for = K 0 33 , we haveΨ > 0. As K 33 is increased, while in the absence of the → 3 1 coupling Ψ becomes more positive, we observe that, with = K 20 31 , Ψ decreases and crosses zero around = K 4 33 , and continues to decrease until = K 10 33 , where Ψ reaches minimum. Increasing K 33 further leads to an increase in Ψ and it becomes positive about = K 16 33 . This behavior can be seen more clearly by plotting the cross section of Ψ at = K 20 31 , as shown in figure 3(a) as the lower curve. For comparison, the MSF versus K 33 for = K 0 31 is also shown in figure 3 (a).
Having examined all the two-channel coupling combinations, we find that the above counterintuitive example is not special, but holds for many more cases such as 0.95 0 21 11 , and → → . These phenomena, some being quite counterintuitive, are a consequence of the simultaneous presence of two coupling channels. In the following we carry out a detailed analysis to understand these phenomena heuristically. Since the MSF can only inform us about the local stability of the synchronous state, the global stability must be verified by numerics. Therefore, we carry out direct numerical simulation of the synchronization dynamics to verify the counterintuitive behavior as exemplified in figure 3(a).
Note that, since the role of the coupling matrix in affecting the stability of the existed synchronization state is reflected only by the eigenvalues, typically it is sufficient to consider a system of two coupled oscillators, for which the coupling matrix is simply and the eigenvalues are 0 and 2. We thus have = ɛ K 2 1 1 and = ɛ K 2 2 2 , and the two oscillators will synchronize if Ψ ɛ ɛ < (2 , 2 ) 0 1 2 . However, due to the complexity of the synchronization boundary in the two-dimensional plane K K ( , ) 1 2 , the synchronization condition for a real coupled dynamical system with > N 2 identical oscillators can be highly nontrivial. In particular, the system will synchronize only if all the normalized coupling points K K ( , ) In our simulations, unless otherwise specified, we consider a random-network system of 100 identical oscillators, where the coupling configuration is determined by the probability p that any pair of two nodes are connected. To be concrete, we set p = 0. and thus W = 0. Otherwise, W has nonzero values. In our simulation, we first run 4 × 10 3 time units to eliminate transient behavior and then use the subsequent 10 3 time units to calculate W. An ensemble of 100 random initial conditions is used, which yields the average value W and the standard deviation. We use the same setting for all subsequent simulations of synchronization dynamics in this paper.
For the two-channel coupling scheme → → 3 1 and 3 3, we fix ɛ = 31.0945 31 so that μ ɛ = 20 31 2 . The corresponding cross section Ψ versus K 33 for = K 20 31 is shown in figure 3(a). From figure 2(a), we see that, by varying ɛ 33 , all the coupling points K K ( , )  large critical point, it is the first triangle point (corresponding to μ 2 ) that moves out of the negative region first, thus the critical coupling point is determined solely by μ 2 . The simulation results are shown in figure 3(b), we see that the synchronizable region of the system matches well with the region in which the MSF assumes negative values. This provides direct evidence for the 'counterintuitive' synchronous behavior under the two-channel coupling configuration. Note that, for μ ɛ = 25 31 2 , from figure 2(a) we see that when ɛ 33 is decreased, the first point moves out the negative region can be another point. In this case, the critical coupling follows a different criterion.
Next, we examine the MSF from the two-channel coupling scheme → → 2 1and1 1. The contour plot of the MSF is plotted in figure 4(a). For the → 2 1single-channel coupling configuration, Ψ is always positive and increases with the normalized coupling strength, as shown in figure 4(a) when = K 0 11 . Thus for this singlechannel configuration synchronization is ruled out. With respect to the → 1 1 coupling, it has a finite negative region of Ψ < 0, as shown in figure 4(a) when = K 0 21 . Therefore, one may think that a finite K 21 value may not help the system to achieve synchronization. However, we find that, as the → 2 1coupling channel is introduced, for a finite value of K 21 the negative region for the → 1 1 configuration is enlarged. As a result, for some value of K 11 (i.e., 6.0) for which the system cannot be synchronized in the single-channel configuration, synchronization can now be realized with the introduction of the → 2 1coupling channel. From figure 4(a) and also more cases shown in figure 1, the boundary separating the synchronous and unsynchronized regions appears quite regular. Here we provide a heuristic argument that in some limiting cases, the boundary can be approximately determined, as exemplified in figure 4(b).
For the → → 2 1and1 1coupling scheme, we have  If  DF s ( )can be approximated as a constant matrix, the Lyapunov exponents will be the real part of the eigenvalues of  DF. Then a zero largest Lyapunov exponent means that the determinant of  DF is zero: 11 21 or ). 21 11 The time average value of the last term is 0.02, leading to α = − K K 0.98 21 11 . Hence, the case of constant  DF s ( ) matrix can be used as a good approximation of the boundary. In fact, the fitting curve for the lower boundary of figure 4(a) is 21 11 , which agrees with the above formula reasonably well since α = 0.2. An alternative assumption is to use the average value of the dynamical variables, = = − x y z s ( , , ) (0.1650, 0.8252, 0.8252) in  DF, and calculate the largest eigenvalue for a given point K K ( , ) 21 11 to approximate the largest Lyapunov exponent. The results are shown in figure 4(b). Comparing with figure 4(a), we observe some similar patterns, despite some difference as large as 20%. Again, this originates from the independence of the dynamical variables of the leading terms relating K 21 and K 11 in the expansion of 0. To demonstrate the above effect via direct simulation, we consider three cases: = K 0, 0.4, 0.8 21 , and plot Ψ versus K 11 , as shown in figure 4(c). The enlargement of the negative region can be seen as K 21 is increased from zero. The corresponding simulations using two coupled oscillators with coupling matrix (7) are shown in figure 4(d), which agree with the MSF results well, directly verifying the counterintuitive role of the → 2 1 channel in inducing synchronization. The MSF of the corresponding two-channel coupling configuration is shown in figure 5(a), where the white solid curve is the contour line of Ψ = K K ( , ) 0 12 21 , which indicates the boundary of the synchronous region. Taking = K 15 21 and = K 50 12 as an example, we see that the MSF for each single-channel configuration is negative, but it becomes positive when both channels are simultaneously present, leading to a loss of synchronization. To gain further insights, we plot Ψ Ψ + K K ( ) ( ) 12 21 in figure 5(b) with the rectangular box specifying the region in which either the single-channel or double-channel coupling configurations can lead to synchronization. We see that the negative region deviates significantly from that due to a naive combination of the negative regions from the isolated single-channel configurations. To be more specific, we fix = K 11 21 and vary K 12 , which corresponds to the cross section as indicated by the dashed line in figure 5(a) and also shown in figure 5(c). We also include the case of = K 0 21 . . Now by increasing K 12 , we see thatΨ K ( ) 12 increases and becomes positive around = K 40 12 , as shown in figure 5(c). Therefore, for fixed ɛ 21 such that μ ɛ = 11 N 21 , the synchronization condition is μ ɛ < 40 N 12 . This indicates the effect of combining the two single-channel coupling configurations in suppressing synchronization. We also carry out direct simulations using the system of 100 coupled Lorenz oscillators, as shown in figure 5(d). There is a good agreement between the direct simulation results that the prediction from the MSF behaviors.

Lorenz system
There are also examples where, a region that cannot be synchronized for either single-channel coupling configuration, when both are present, becomes synchronizable. For example, for the → 3 3 coupling channel, the MSF has three cross points with the K axis: 1.368, 9. . As the other coupling strength ɛ 33 from 0 is increased to reach the value of 13.7724, i.e., ɛ ɛ = 1.1073 33 13 , the last normalized coupling point (the down triangle), which corresponds to μ 2 , enters the negative region. Consequently, all the normalized coupling points in between the down triangle and the upper triangle enter the negative region. The system is synchronizable. As ɛ 33 is increased further, e.g., to ɛ = 21.2550 33 so that ɛ ɛ = 1.7089 33 13 , the upper triangle (corresponding to μ N ) moves out the negative region, and synchronization is lost. Figure 6(b) shows the simulation results, which agree with the above analysis well. Note that, since the system is evolved into a fixed point, we only run ten time units to reduce the transient effects and 10 time units to calculate W with an ensemble average over 200 random initial conditions. Here the simulation is carried on a network with 100-oscillators. It would be interesting to see how does synchronization evolve with network size or network structure (scale-free, small-world etc). Since the MSF method separates the nodal dynamics and the coupling network structure, the effects of network structure are revealed in its eigenvalues, especially μ 2 and μ N if the synchronizable region has a simple concave shape in the normalized coupling parameter space. Since both μ 2 and μ N depend on network structure or network size [31], and because of the complex interplay between network structure and the MSFs of dynamical systems, there could be certain cases that some network structures are better synchronized via single-channel coupling for a given dynamical system, and some other structures by multi-channel coupling. However, the straightforward extension to say that certain networks are better synchronized via single-channel or multi-channel coupling regardless of the dynamical systems is in general not true. A detailed examination of all the eigenvalues (or all the points on the normalized coupling parameter plane) will be necessary to predict synchronization.
Since the MSF Ψ, the largest transverse Lyapunov exponent averaged over the global trajectory, is determined by the largest local transverse Lyapunov exponent at each trajectory point, to gain insights, we plot the largest local transverse Lyapunov exponents along a typical trajectory from the chaotic Lorenz attractor and compare the one-channel coupling scheme → 2 3for = K 18 23 with the two-channel configuration → → 2 3and3 3with = K K ( , ) (18, 20) 23 33 . The results are shown in figure 7. We see that, when only the single-channel coupling → 2 3is present, most of the local transverse Lyapunov exponents are positive, as shown in figure 7(a), leading to a positive global transverse Lyapunov exponent. When the → 3 3 coupling channel is added, as shown in figure 7(b), for the same trajectory the local transverse Lyapunov exponent becomes mostly negative for the trajectory, providing insights into the behavior of Ψ.

Abnormal enhancement of synchronization of coupled HR neurons
Multi-dimensional dynamical systems are often invoked as models of real phenomena, where typically adding dimensions increases the 'realisticity' of the model. For instance, many biological functions, such as gene regulation or neural interactions, can be modeled via multi-dimensional dynamical systems, which can then be coupled via various schemes. To be specific, it has been well recognized that the coupling between neurons are often involved with different effects, e.g., the electric gap junction and the diffusion of the extracellular ions [32][33][34][35]. It has been found that the electrical and chemical synaptic connections of different types of inhibitory neurons are specific, and may allow each inhibitory network to function independently [36]. Odor representations in the olfactory bulb were found to be stabilized by interneurons that were densely coupled to the output neurons by electrical and GABAergic synapses [37]. In particular, Hindmarsh and Rose proposed a mathematical model of the neuron spiking, which provided a qualitative description of the membrane potential  1.831 23 33 . and transport of ions in a single neuron bursting [38]. Synchronization of HR neurons has thereafter attracted much attention [39][40][41], and experimental observation of synchronization between neurons has been obtained via electrical coupling [42]. Most of the phenomenal studies of synchronization of dynamical neuron models assume one channel coupling. However, as we demonstrate below, when more channels are employed in the coupling between neurons, their synchronization behavior can be intriguing and beyond the expectation from the result of single channel couplings. The dynamics of HR neuron is given by [ , where I = 3.2 is the external current input, r = 0.006, s = 4. x represents the membrane potential, y is the spiking variable representing for the transport of sodium and potassium ions that is made through fast ion channels, and z is for the transport of other ions through slow channels. The Jacobian matrix is The neurons, when coupled by the slow ion channel z, cannot get synchronized for any strength of the coupling parameter. When the coupling is from the fast ion channel y to the membrane potential x, there is only a small synchronization window in the normalized coupling parameter K, i.e.,[0.286, 1.233] [14]. However, when both coupling channels are on, as shown in figure 8, the synchronizable region is greatly enhanced.

Conclusion and discussion
To summarize, we have systematically analyzed the MSFs of the coupled chaotic Rössler and Lorenz oscillators for all possible two-channel coupling configurations, which can serve as the library and provide guidelines for future research concerning synchronization behavior under multiple-channel coupling configurations. While for most configurations, the effects of the two channel couplings can lead to expected behaviors of synchronization in accordance to each of the isolated cases, there are situations where the interplay between the two coupling channels can lead to counterintuitive phenomena. For example, two coupling channels, each when isolated leading to no synchronization, can generate global synchronization when both are present simultaneously, and vice versa. We have demonstrated various aspects of the counterintuitive synchronous behavior for a large number of two-channel coupling configurations, and provided an understanding of such behavior based on examining the time averaged values of the dynamical variables for situations where the leading terms of K 1 and K 2 are independent of the variables. In general, synchronization behaviors predicted by the MSF for any isolated, single-channel coupling configuration, may not be expected when any combination of two such channels are present. Due to the inherent nonlinearity of the system, deeper insight (e.g., analytic insight) into these phenomena is still lacking. Our result can serve as an atlas to greatly facilitate further research on chaotic synchronization dynamics under multiple-channel coupling configurations. Indeed, while there is a large body of literature on chaotic synchronization, vast majority of the previous works focused on the setting of single-channel coupling. There has been little understanding of nonlinear synchronization dynamics for multiple-channel coupling, and our present work has filled this gap. More specifically, our extremely detailed and systematic probing of the parameter space structure provides, for any given system under multichannel coupling, a quick assessment of the likelihood of synchronization with minimal amount of information, namely, the eigenvalues of the coupling matrix. In addition to numerical computations, we have provided a physical understanding of the parameterspace structure by approximating the Jacobian matrix to calculate the largest Lyapunov exponents and by analyzing the behaviors of the local transverse Lyapunov exponents. Furthermore, our result for HR neurons indicates that, multi-channel coupling, which is believed to be more realistic especially in coupled neuron systems, can enhance significantly the synchronizable region in the normalized coupling parameter space.
Our work illustrates that, in the study of synchronization dynamics in systems of coupled nonlinear oscillators, the extension from single-channel to two-channel or multiple-channel coupling configurations can be highly nontrivial. In fact, new and significantly richer phenomena of synchronization can emerge as a result of providing additional coupling channels to the system, a demonstration that nonlinear dynamical systems never stop to present us with surprises, justifying further research even in the relatively well studied area of synchronization. In addition, synchronization is the simplest form of collective behavior. Which other kinds of new collective behaviors emerge in the coupled oscillator systems, particularly in the parameter regions away from the simple synchronization? It would be interesting to see how does the spectrum of collective behaviors depend on the coupling scheme (even more than how it depends on the interaction strengths/parameters). This can be a promising avenue for the coupled dynamical systems and deserves future investigation.