Synchronization in Dynamically Coupled Fractional-Order Chaotic Systems: Studying the Effects of Fractional Derivatives

0is study presents the effectiveness of dynamic coupling as a synchronization strategy for fractional chaotic systems. Using an auxiliary system as a link between the oscillators, we investigate the onset of synchronization in the coupled systems and we analytically determine the regions where both systems achieve complete synchronization. In the analysis, the integration order is considered as a key parameter affecting the onset of full synchronization, considering the stability conditions for fractional systems. 0e local stability of the synchronous solution is studied using the linearized error dynamics. Moreover, some statistical metrics such as the average synchronization error and Pearson’s correlation are used to numerically identify the synchronous behavior. Two particular examples are considered, namely, the fractional-order Rössler and Chua systems. By using bifurcation diagrams, it is also shown that the integration order has a strong influence not only on the onset of full synchronization but also on the individual dynamic behavior of the uncoupled systems.


Introduction
Synchronization is an emergent physical phenomenon caused by the interaction of two or more dynamic entities that pervade the natural world [1][2][3][4]. In the case of oscillating units, the synchronization phenomenon can be defined as the adjustment of temporal evolution to a common rhythm.
For the case of integer-order systems, there exists a vast and mature literature where we can find different interconnection schemes for synchronizing dynamic systems, like, for example, master-slave synchronization scheme, adaptive synchronization, and synchronization based on state observers, to name a few [5][6][7][8][9][10]. Although each of these strategies is effective, there are limitations in their applications, e.g., there are cases where these schemes have marginal ranges for which the synchronous response is achieved or have poor robustness to maintain a stable synchronous state under the influence of external disturbances. is is one of the reasons why dynamic interconnections have emerged as an alternative to the classical static schemes. In this case, the interaction between agents is indirectly achieved through a suitably designed dynamic coupling [11][12][13]. is type of synchronization strategy has shown better performance than static couplings. In particular, dynamic coupling increases the intervals of coupling strength values for which it is possible to achieve synchronized behavior, and it may also be possible to synchronize systems that cannot be synchronized with static coupling [11].
On the other hand, the use of fractional calculus has been extensively studied in nonlinear systems (see, e.g., [14][15][16][17][18]) and also, there exist notable contributions related to the study of synchronization in fractional-order systems (see, e.g., [16,[19][20][21][22][23]). For example, there is work based on applying sliding modes to fractional-order models to achieve synchronization [24][25][26][27][28][29]. e modeling and analytical study of fractional-order systems is also a fruitful field, e.g., the use of the Razumikhin approximation for fractional-order systems with delay [30,31], the extrapolation of Lyapunov theory to fractional systems [32,33], and the existence and uniqueness of equilibrium points of the Mittag-Leffler criteria [34,35]. However, the use of dynamic couplings in the context of fractional-order systems seems to be unexplored so far.
Consequently, in this study we present a synchronization scheme of fractional order based on dynamic coupling. In particular, a master-slave interconnection is considered, in which the systems have an indirect interaction via a fractional-order linear system. In the analysis, the Rössler equation [36] and the Chua double-scroll oscillator [37] are considered. Among the questions to be addressed is whether a dynamic coupling designed for integer-order systems is also effective in its fractional-order version? If so, how does the derivative order influence the occurrence of synchronization in the coupled systems? e local stability of the synchronous solution in the coupled systems is investigated by analyzing the error dynamics, and furthermore, the onset of synchronization is also numerically investigated by computing some statistical metrics like Pearson's correlation between time series. Additionally, using bifurcation diagrams, we have investigated the dynamic behavior of the uncoupled systems. e obtained results show that the integration order has a strong influence on the stability of the synchronous solution, and interestingly, it also produces a period-doubling cascade route to chaos in the uncoupled systems.
e rest of the study is organized as follows: Section 2 describes the basics of fractional calculus and gives a brief introduction to fractional-order systems. en, Section 3 describes the proposed synchronization scheme and the local stability of the synchronous solution in the coupled systems is discussed. Subsequently, in Section 4, the performance of the dynamic coupling is investigated using the Rössler and Chua systems as application examples. Finally, Sections 6 and 7 are dedicated to the discussion and conclusions, respectively.

Preliminaries
is section presents a brief overview of some basic concepts about fractional-order systems. In particular, the Caputo derivative, the general representation of a fractional-order system, and the stability of linear time-invariant fractionalorder systems are revisited.

Fractional Caputo Derivative.
In the literature, there are various definitions of fractional-order derivatives, the most common being the Riemann-Liouville and Caputo operators [38,39]. e fractional Caputo derivative for a timeinvariant system described by the vector field f(x) is defined as with n � q for the integration order 0 < q < 1, being Γ the gamma function defined as follows: 2.2. Fractional-Order Dynamic System. A commensurate fractional-order time-invariant system can be described, in general, as follows: subject to initial conditions where n 1 , n 2 , . . . , n k are rational numbers, such that n k > n k−1 > · · · > n 1 > 0, n j − n j−1 ≤ 1 for all j � 2, 3, . . . , k, and 0 < n 1 ≤ 1. e least common multiple of the denominator of n 1 , n 2 , . . . , n k is defined by M and set q � 1/M and N � Mn k . en, equation (3) can be expressed as follows [38]:

Stability of Time-Invariant Fractional-Order System.
A linear time-invariant fractional-order system is described by where x ∈ R n is the state vector, A ∈ R n×n is a constant matrix, and 0 < q < 1 is the fractional commensurate derivative order. en, the stability of the system described by equation (6) is determined as follows [38]: (i) System (6) is stable, if and only if |argλ j | ≥ qπ/2, ∀j � 1, 2, . . . , n (ii) System (6) is asymptotically stable, if and only if |argλ j | > qπ/2, ∀j � 1, 2, . . . , n (iii) System (6) is unstable, if and only if |argλ j | < qπ/2, for at least one λ j , j � 1, 2, . . . , n From these conditions, it is clear that the local stability of fractional-order systems depends on the integration order q, so that the stability of an equilibrium point can be modified by the fractional order, and therefore, the stability region at the complex plane is as in Figure 1 [40].

Proposed Synchronization Scheme Based on
Dynamic Coupling e dynamic interconnection considered in this study has been presented in [11] for integer-order systems and is adapted here to the fractional-order case. e scheme, where the interaction between the systems is indirect via a dynamical system ( Figure 2) is described by the following set of equations: where x m , x s ∈ R n represents the state vectors of both the master and slave systems and h � (h 1 , h 2 ) T for h i ∈ R, i � 1, 2 is the state variables of the dynamic coupling. It is assumed that the vector field F is smooth enough, which can be either linear or nonlinear, and the coupling force between the systems is denoted by k.
On the other hand, the design of a dynamic coupling involves two coupling matrices, denoted B 1 ∈ R n×2 and B 2 ∈ R 2×n . ese matrices are generated under the premise that only one of the elements of each of these matrices is equal to 1, and the other entries are zero, which means that the coupling is applied only in one state variable of the slave system and that the coupling considers only one measured variable.
Finally, the matrix G from equation (7) is given by where c 1 , c 2 , and α c are design parameters of the dynamic coupling. e construction of the coupling system, for the integer-order case, is inspired by the so-called Huygen's coupling [11,41], which in its simplest form can be interpreted as a damped oscillator.
Since in this study the focus is on synchronization, it is necessary to give the following definition.
e coupled systems, equation (7), are said to be asymptotically synchronized if 3.1. Local Stability Analysis. In order to investigate the stability of the synchronous solution defined in equation (9), we proceed as follows. First, it is assumed that function F in equation (7) can be written as the sum of linear and nonlinear components, i.e., it is assumed that where P ∈ R n×n is a constant matrix and E(x i ) ∈ R n is a vector containing nonlinear terms. Next, the synchronization error is defined as e p :� (x m − x s , h) T . Note that in the definition of the error, we have included the state h of the dynamic coupling. e reason is because the parameters in the dynamic coupling should be chosen such that, when the systems synchronize, the coupling vanishes. en, by replacing equations (10) into (7), and computing the corresponding synchronization error dynamics, we obtain where where O � (0, 0) T . Furthermore, note that the term g p (t, e p ) is a vanishing perturbation [42] because g p (t, 0) � 0. en, the stability properties of system equation (11) are fully determined by the eigenvalues of the matrix A. In particular, following the results presented in Section 2.3, we have that the synchronization error dynamics (equation (11)) is locally asymptotically stable if  us, if it is possible to find values of k, c 1 , c 2 , and α c such that the above condition is satisfied, then the coupled systems described by equation (7) will achieve complete synchronization, according to Definition 1.

Statistical Metrics Used for Detecting Complete
Synchronization. In this study, the onset of synchronization in the coupled systems (equation (7)) is also numerically studied by computing the following synchronization index.
where n is the dimension of the systems to be synchronized, |e i | is the absolute value of the synchronization error between the i-th state variable of the master and slave systems, i.e., and C i is the Pearson correlation [43], computed from the i-th state variable of the master and slave systems described by where σ N mi N si is the covariance between the data obtained from the time series of the state variables of the master and slave systems, and σ N mi σ N si is the standard deviation obtained from the i-th state variable of the master (slave) oscillator. Finally, when S � n, the systems are synchronized.
In the next section, the onset of synchronization with dynamical coupling is studied for two particular fractionalorder chaotic systems, namely the Rössler and the Chua oscillators.

Application Example 1: Rössler System
e fractional-order version of the well-known Rössler system [44] is given by where x i , i � 1, 2, 3 denotes the state variables and a, b, c are constants.
It has been reported that every system has a limit of integration order for which it is possible to use a fractionalorder derivative without stabilizing its dynamics [16,17,45].
is could be interpreted to mean that the dynamics described in equation (17) are restricted to have at least one eigenvalue in the unstable region, and this being true only if |arg(λ)| < qπ/2 for at least one of its eigenvalues (λ). With a � 0.2, b � 0.2, c � 5.7, the eigenvalues obtained by the Jacobian matrix evaluated at the equilibrium point E 1 � (0.007, −0.0351, 0.0351) are λ � [0.0970 ± 0.9952i; − 5.6870]. Since |min(arg(λ))| � 1.4736, the critical order of integration is defined as q c � |min(arg(λ))|(2/π) (Section 2), and by substitution it is obtained as q c � 0.9381. is result is confirmed by the bifurcation diagram as shown in Figure 3, where the local maxima in x 1 are plotted as a function of the variation of the integration order. Note that the derivative order is the same for all state variables of the system. e inset shows the attractors and the corresponding integration order considered for the developed analysis.
It is worth noting that the bifurcation diagram shown in Figure 3 was numerically calculated using the Adams-Bashforth-Moulton (ABM) method [46] and following the guides for plotting a proper bifurcation diagram [47].

Dynamically Coupled Fractional-Order Rössler Systems.
Now, we consider a pair of identical Rössler systems [44] described by where x m,s denotes the state variables of the master and slave systems, respectively, h 1,2 denotes the states of the dynamic coupling, k indicates the coupling force between the oscillators, a, b, c are constants of the Rössler model, and q j , j � m, s, d denotes the integration order of the master, slave, and dynamic coupling, respectively. Here, we consider the case in which all orders of integration are equal, i.e., q m � q s � q d . en, the dynamic behavior of system equations (18)-(20) is numerically studied as a function of the coupling strength k and the integration order q. For this purpose, equations (18)- (20) are numerically integrated with the following parameter values cf. [11]: On the other hand, the coupling strength is varied in the interval 0 ≤ k ≤ 20 with a test size Δ k � 0.2 and the integration order is varied in the interval 0.96 ≤ q ≤ 1 at Δ q � 2e − 4 .
e obtained results are shown in Figure 4(a), where the colors indicate the value of the synchronization index S (see equation (14)). Synchronous behavior is indicated by the dark blue color (S � 3), while the remaining colors indicate unsynchronized dynamics.
Moreover, it is evident from Figure 4(a) that there are regions where, for a fixed coupling strength, the limit behavior is completely determined by the integration order. For example, for a fixed coupling strength of k � 10, the synchronization index S abruptly changes when the 4 Complexity integration order q is within the interval 0.988 < q < 0.995, as can be seen in Figure 4(b), but if the coupling strength is increased up to k � 15, the systems always achieve complete synchronization.
To validate the synchronization regions obtained by the time series analysis, we conduct a stability analysis following the results presented in Section 3.1. First, let the following synchronization errors be defined as follows: e j � x jm − x js , for j � 1, 2, 3, and e 4 � h 1 , e 5 � h 2 . en, the resulting error dynamics are given by Note that the error dynamics in equation (21) can be written in the form of equation (11) Note that the term g p (t, e p ) in equation (21) is indeed a vanishing perturbation since when the systems synchronize it follows that x 1m � x 1s , x 2m � x 2s , and x 3m � x 3s and therefore g p (t, 0) � 0. Consequently, the local stability of the synchronous solution in the coupled systems described by equations (18)-(20) can be determined from the condition in equation (13). In particular, we compute condition equation (13) as a function of the integration order q and the coupling strength k. e obtained results are shown in Figure 5(a), where the blue region corresponds to values of k and q for which condition (13) is not satisfied and thus the synchronous solution is unstable, whereas on the white region, condition (13) is satisfied and then the synchronous solution is expected to be stable. For the sake of comparison, Figure 5(b) shows the overlap of Figure 4 with Figure 5(a). It can be seen that there is a good agreement between the numerical and the analytical results.

Remark 1.
In the previous analysis, we have considered the case where the integration order of the systems and the dynamic coupling are the same. However, we also have conducted a numerical study in which the integration orders are different. In particular, we have numerically integrated equations (18)-(20) using the parameter values considered before, as a function of the integration orders of the master and slave systems, while the integration order of the dynamic coupling remains fixed. e integration orders q m , q s are varied in the interval 0.98 ≤ q j < 1, j � m, s, and considering the integration order q d � 0.985. e obtained results are shown in Figure 6(a) where the blue areas correspond to synchronization (S � 3). From the obtained results, it is clear to see that the integration order of the master and slave systems should be almost the same to observe a synchronized behavior and that larger differences are tolerated as long as the integration order of both systems approaches to one.
On the other hand, Figure 6(b) shows the obtained results for the case that only the integration order of the dynamic coupling is varied, while the oscillators are assumed to have integer order, i.e., q m � q s � 1. In this case, equations (18)- (20) are numerically integrated by varying the integration order q d of the dynamic coupling in the interval 0.8 < q d < 1 and the corresponding synchronization index, given by equation (14), is calculated. e obtained results are shown in Figure 6(b) for two different coupling forces, where the choice of these k values corresponds to those reported in [11] and those obtained in Figure 4(b).

Application Example 2: Chua System
If the Chua system described in [37] is modeled with derivatives of fractional order, then the system described in equation (23) is analyzed.
where D q is the fractional-order derivative by Caputo's definition, x i , i � 1, 2, 3 denotes the system state variables, σ, β are constants of the Chua circuit, and ϕ(x 1 ) is a nonlinear function defined in equation (24), with constant values a, b.
In the same way as for the Rössler fractional-order model, the system described in (23) is analyzed to identify the minimum fractional order that can be modeled without stabilizing the dynamics, namely, q c � 0.9541, since is result is confirmed by the bifurcation diagram as shown in Figure 7, where the local maxima in x 1 are plotted as a function of the integration order variation. e inset shows the attractors and the corresponding integration order.

Dynamically Coupled Fractional-Order Chua Systems.
In the same way as for the Rössler system described above, a pair of Chua oscillators [48] coupled by dynamical connections and defined by the system of equations (25)-(27) is considered.
master system where x m,s denotes the state variables of the master and slave systems, respectively, h 1,2 is the states of the dynamic coupling, ϕ(x 1m,s ) is the nonlinear function defined in equation (28), and k is the coupling force between the systems. 6 Complexity As reported in [11], the following values are used in this study for all analyzes developed; σ � 10, β � 14.87, a � −1. To identify synchronization regions in the coupled Chua systems of fractional order, an analysis of the coupling force as a function of the integration order is developed. e obtained results are shown in Figure 8. Note that the color map represents the value of the metric S and that the darkest shade of blue represents S � 3, which means that the systems have reached full synchronization.
Analogous to the stability analysis performed for the coupled Rössler system, the dynamic error model of the coupled pair of Chua oscillators described in equations Due to the nature of the nonlinearity of the Chua circuit, it is not possible to perform the same analysis as in the Rössler system; instead, it is necessary to use the Jacobian of the error model, equation (30), evaluated in one of the equilibrium points of the system. Since the Chua circuit has symmetric equilibrium points located at E 1 � (−1.841, 0.0004474, 2.179) for the previously defined values, the choice of one of these points does not affect the analysis.
where p(x) � 1/2(b sign(|x| − 1) + b) − 1/2(a sign(|x| − 1) − a)). After defining the system shown in equation (30), it is possible to perform the stability system analysis, where the stability of the dynamic model of the Chua coupling error is sought by modifying the integration order and the coupling force. e analytical result is shown in Figure 9(a) and then compared with the map obtained from the time series  (14) as a function of integration order, for a fixed k � 13 and k � 20. A value of S � 3 indicates complete synchronization. 8 Complexity analysis in Figure 8(b). As with the Rössler case, the analytical result is able to describe the boundary at which the system is unstable.

Remark 2.
Similarly, to the case of the Rössler systems, we have also investigated the onset of synchronization in the fractional-order Chua oscillators as a function of the derivative orders in the oscillators while keeping the coupling system at fixed q d , resulting in the map as shown in Figure 10(a). e results from the analysis of the behavior of the Chua systems under the dynamic coupling integration order variation are shown in Figure 10(b), where the integration orders of the oscillators remain fixed at q m,s � 1.

Discussion
From the results presented in this study, it is possible to confirm the research question formulated in the introduction of this study, according to which the use of dynamic couplings in chaotic systems of fractional order is able to induce complete synchronization, as it has been reported in their counterparts of integer order. Likewise, it is noteworthy to mention that the transition from unsynchronized behavior to synchronization is not abrupt, since there exists a region where the coupled systems may exhibit some sort of intermittency phenomenon. ese areas are indicated by the blurred areas in Figures 4 and 8.
If the oscillators are modeled with fractional derivatives but a fixed integration order is maintained in the coupling system, the desired synchronous behavior is achieved only under the condition that both oscillators have the same integration order. Small variations in the integration order in some of the models cause the systems to lose their synchrony, as shown in Figures 6(a) and 10(a). In contrast, the synchronization seems to have some robustness against variations in the integration order of the dynamic coupling, provided that the oscillators have the same integration order, as shown in the numerical results presented in Figures 6(b) and 10(b).
It is also noteworthy that in the (q, k) plane there exist regions where the systems are easier to synchronize. is is explained by the bifurcation diagrams of the isolated oscillators, given in Figures 3 and 7, where the modification in the integration order causes important qualitative changes in the dynamics of the system, where both models are able to present chaotic or periodic behavior, for a set of parameters where the integer-order dynamics is always chaotic, only due to the modification of the derivative order. is is not only an indication that the synchronization between the systems requires lower coupling forces for periodic and quasiperiodic behaviors but also a clear indication that the change in the integration order can be associated with a modification of the vector field, which can also be achieved in the integerorder system by modifying the system parameters [17].
e local stability analysis is in good agreement with the numerical analysis, as shown in Figures 5 and 9. It should be noted, however, that the stability conditions are only necessary conditions. is result is similar when talking about the stability of the equilibrium points in a chaotic system, where obtaining unstable saddle points with index 2 favors the occurrence of chaotic behavior but does not guarantee the occurrence of a strange attractor in the system [49].
It is worth mentioning that in the cases of analysis where the equations to solve do not have the same derivative order, the algorithm proposed by Petráš is implemented [50]; otherwise, the Adams-Bashforth-Moulton (ABM) method is used [46], which is a generalization of the classical ABM integrator that is well known in the resolution of first-order switching system problems [40,51].
Notice that the results presented here have been obtained under the assumption of identical oscillators. It is still necessary to determine the robustness of the dynamic coupling against parameter mismatches or external disturbances in the oscillators.

Conclusions
We have analyzed the onset of synchronization in fractionalorder chaotic systems interacting via a linear time-invariant dynamic coupling, which also is described by fractional derivatives. e obtained results have shown the ability of the dynamic coupling to induce synchronization in the systems, and the strong influence of the integration order on the onset of synchronized behavior has been demonstrated.
Among the observed limitations is that dynamic coupling is sensitive to variations in the integration order of the systems. erefore, the master and slave systems should have the same integration order. Furthermore, it has been shown that the linearization approach used here to study the local stability of the synchronous solution only provides necessary conditions. Further investigation is needed to derive stronger stability conditions. Perhaps the use of transverse Lyapunov exponents can solve this problem. It remains as future work to extend these results to the bidirectional case and also to the case of networks. Also, it would be interesting to investigate whether any emergent behavior or other types of synchronous behaviors can occur in the coupled systems depending on the integration order of both the systems and the coupling.
Finally, we would like to point out that the results presented here apply to the Caputo definition of fractional derivative. Moreover, these results were obtained using two different numerical approaches, namely, the ABM and the Petráš integrator method. As future work, we plan to compare these results with different fractional operators, such as the Riemann-Liouville operator or the Atangana-Baleanu operator, to name a few.

Data Availability
e data used to support the findings of the study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.