Dominant Attractor in Coupled Non-Identical Chaotic Systems

The dynamical interplay of coupled non-identical chaotic oscillators gives rise to diverse scenarios. The incoherent dynamics of these oscillators lead to the structural impairment of attractors in phase space. This paper investigates the couplings of Lorenz–Rössler, Lorenz–HR, and Rössler–HR to identify the dominant attractor. By dominant attractor, we mean the attractor that is less changed by coupling. For comparison and similarity detection, a cost function based on the return map of the coupled systems is used. The possible effects of frequency and amplitude differences between the systems on the results are also examined. Finally, the inherent chaotic characteristic of systems is compared by computing the largest Lyapunov exponent. The results suggest that in each coupling case, the attractor with the greater largest Lyapunov exponent is dominant.


Introduction
Chaos theory is an inviting subject of study in nonlinear dynamics. The concepts in line with chaotic dynamics, such as strange attractors, fractal dimensions, Lyapunov exponents, and synchronization, have been extensively addressed in every research field [1,2]. For example, in neuroscience, the presence of deterministic chaos has been proved both experimentally and theoretically [3,4]. Networks consisting of chaotic oscillators have also received considerable attention over the past decades [5,6].
In the context of networks, coupling and its configuration are inseparably studied. Adjustment of coupling strength may result in a correlated in-time behavior, universally known as synchronization [7]. Such a collective state of coupled systems indicates dynamical interactions, having a critical role in many natural networks [8]. For example, it is vital in nervous system tasks such as the human memory process or enhancement of live verbal communication [9][10][11]. Regarding chaotic systems, two trajectories starting from initial conditions tend to diverge exponentially over time [12,13]. As a result, it was first realized that chaotic systems intrinsically defy synchronization until, in 1990, Pecora and Carroll [14] proved quite the contrary. They showed that synchronization is achievable by drive-response coupling of two identical chaotic systems [12].
Chaotic systems exhibit various forms of synchronization such as phase [15], antiphase [16], lag [17], explosive [18], cluster [19,20], generalized [21], and complete [22]. Complete synchronization, the perfect hooking of chaotic trajectories, is only achievable in identical systems [13,23]. In such systems, attractors structures begin to assimilate until they match completely [24]. However, when oscillators are non-identical, even the classic definition of the synchronization manifold is not analytically valid [23]. Consequently, few studies have discussed this collective behavior in non-identical chaotic systems [25][26][27]. Our results may apply to any conditions where the emergence of weak connections between different oscillators could impair their primary function. As each oscillator's inherent dynamics is considered, any factor impairing it would be functionally important. For instance, the brain's connectivity is a prominent feature in the brain based on the neuroanatomy. Circuits consisting of various types of neurons are observed throughout the nervous system [33]. An example of connectivity among a diverse set of neurons is corticosubcortical basal ganglia circuits that carry out the vital task of movement modulation. In simple terms, if we take the basal ganglia as a controller based on its function (facilitating voluntary and inhibiting unwanted movements [34]), this circuit is a closed loop in which cortical neurons provide the input of the controller, and the desired response is inserted back into the frontal lobe through the thalamus [35]. Any abnormal communication between neural elements may lead to the malfunction of the whole circuit. The term abnormal means stronger or weaker than its functioning state. In the case of Parkinson's disease, the strengths of former connections are altered, which leads to tremors at rest or a paucity of volunteer movement [36]. Assuming that under pathological conditions, connections in primary pathways are weakened, neighboring neurons exert disturbing inputs to the abovementioned components. Taking the intrinsic dynamics of a neural element in a functioning circuit as primary and its disturbing dynamics as secondary, we can point out some dualities between this scenario and our study. The operative tasks of neighboring neurons can be examined through the coupling of non-identical systems, and low coupling strengths simulate their weak connections. The functional dynamics of oscillators are deemed to be their uncoupled attractors.
In Section 2, as materials of the study, we briefly review the mathematical description, time series, attractor properties of systems, and also the coupling scheme. Then, we introduce the time and amplitude scaling and the comparison methods. In Section 3, we evaluate networks of two non-identical systems and discuss the costs in given coupling strength ranges. The return maps of coupled systems are provided to visualize the relative differences in cost values. By the end of this section, we present the LLEs of systems and their values after scaling. We conclude the results in Section 4.

Materials and Methods
In this section, the used chaotic systems are described mathematically. To investigate the effect of the amplitude and frequency of oscillations, amplitude and frequency scaling are defined further. The end of this section introduces the applied method for evaluating the similarity of attractors of coupled systems with their original attractors.

Coupled Systems
To investigate the effects of different factors on the dominant attractor of two coupled non-identical oscillators, we consider three paradigmatic chaotic systems: the Lorenz system, the Rössler system, and the Hindmarsh-Rose system. The dynamical equations of these systems are as follows: Lorenz system: where the parameters are σ = 10, ρ = 28, and β 1 = 2.
where α = 0.2, β 2 = 0.2, and γ = 9. Hindmarsh-Rose system: where I = 3.2, r = 0.006, and s = 4. Note that the systems at the defined parameters are monostable. Some of the previous studies [37][38][39] have attempted to investigate the synchronization of coupled chaotic systems and find the best coupling variable for optimal synchronization. As, in this study, the aim is not to reach synchronization, for simplicity, we consider that the systems are coupled linearly through x variables. Considering the describing equations of two systems as F = (f 1 , f 2 , f 3 ) and G = (g 1 , g 2 , g 3 , the coupled systems can be described bẏ where d denotes the coupling strength. After coupling two systems, the coupling strength is increased, and the attractors of the two coupled systems are compared with their original ones. The time series of x states of the Lorenz, Rössler, and Hindmarsh-Rose (HR) systems are represented in Figure 1a-c, respectively. It can be observed that the frequency of oscillations is different in these systems. The Lorenz oscillations are faster than the other two, and the HR system is slowest. The time scale of the chaotic oscillations is equal: f L = 0.94 for the Lorenz system, f R = 0.17 for the Rössler system, and f HR = 0.03 for the HR system. frequencies. Thus, the parameter ω is defined as the ratio of frequencies between two systems, i.e., the ratio of the frequency of the faster attractor to the slower one. Then, the frequency scaling is performed by taking dt→dt/ω in the equations of the slower system in Equation (4) (the faster system remains unchanged). Consequently, considering the frequency-adjusting parameter, the frequency of the slower system equates with the faster one. The attractors of the systems in 3D phase space are depicted in Figure 2. The Lorenz, Rössler, and Hindmarsh-Rose attractors are demonstrated in parts a, b, and c, respectively. As can be seen, the size of the attractors is different. As systems are coupled on x states, the impact of such a dissimilarity only in the x state draws our attention. For example, the difference in Lorenz (or Rössler) attractor size in comparison with Hindmarsh-Rose is noticeable. Therefore, when x states of Lorenz and HR systems are connected, the effect of Lorenz on HR is much higher than vice versa. To examine the impact of the x state amplitudes, a scale on amplitude is also performed.
For scaling the domain of the x state in each system, the absolute global maximum of the x state (|maxX A |) and the absolute global minimum (|minX A |) are considered. Then, a scaling parameter a is defined: By converting x → ax in all equations of the system, the x domain becomes a subset of [−1, 1].

Comparison of the Attractors
A suitable similarity indicator is required to detect the error between the original attractor and the one resulting from adding the coupling term. In the literature, some This frequency difference can significantly affect the dynamics of systems after coupling. To examine this effect, we scale the frequencies of systems to have almost the same frequencies. Thus, the parameter ω is defined as the ratio of frequencies between two systems, i.e., the ratio of the frequency of the faster attractor to the slower one. Then, the frequency scaling is performed by taking dt → dt/ω in the equations of the slower system in Equation (4) (the faster system remains unchanged). Consequently, considering the frequency-adjusting parameter, the frequency of the slower system equates with the faster one.
The attractors of the systems in 3D phase space are depicted in Figure 2. The Lorenz, Rössler, and Hindmarsh-Rose attractors are demonstrated in parts a, b, and c, respectively. As can be seen, the size of the attractors is different. As systems are coupled on x states, the impact of such a dissimilarity only in the x state draws our attention. For example, the difference in Lorenz (or Rössler) attractor size in comparison with Hindmarsh-Rose is noticeable. Therefore, when x states of Lorenz and HR systems are connected, the effect of Lorenz on HR is much higher than vice versa. To examine the impact of the x state amplitudes, a scale on amplitude is also performed.
frequencies. Thus, the parameter ω is defined as the ratio of frequencies between two systems, i.e., the ratio of the frequency of the faster attractor to the slower one. Then, the frequency scaling is performed by taking dt→dt/ω in the equations of the slower system in Equation (4) (the faster system remains unchanged). Consequently, considering the frequency-adjusting parameter, the frequency of the slower system equates with the faster one. The attractors of the systems in 3D phase space are depicted in Figure 2. The Lorenz, Rössler, and Hindmarsh-Rose attractors are demonstrated in parts a, b, and c, respectively. As can be seen, the size of the attractors is different. As systems are coupled on x states, the impact of such a dissimilarity only in the x state draws our attention. For example, the difference in Lorenz (or Rössler) attractor size in comparison with Hindmarsh-Rose is noticeable. Therefore, when x states of Lorenz and HR systems are connected, the effect of Lorenz on HR is much higher than vice versa. To examine the impact of the x state amplitudes, a scale on amplitude is also performed.
For scaling the domain of the x state in each system, the absolute global maximum of the x state (|maxX A |) and the absolute global minimum (|minX A |) are considered. Then, a scaling parameter a is defined: By converting x → ax in all equations of the system, the x domain becomes a subset of [−1, 1].

Comparison of the Attractors
A suitable similarity indicator is required to detect the error between the original attractor and the one resulting from adding the coupling term. In the literature, some For scaling the domain of the x state in each system, the absolute global maximum of the x state ( |maxX A |) and the absolute global minimum ( |minX A |) are considered. Then, a scaling parameter a is defined: By converting x → ax in all equations of the system, the x domain becomes a subset of [−1, 1].

Comparison of the Attractors
A suitable similarity indicator is required to detect the error between the original attractor and the one resulting from adding the coupling term. In the literature, some methods have been proposed for quantifying the difference between attractors [40,41]. Rybin et al. [41] proposed a method to estimate the differences between two return maps based on the 2D histogram. Jafari et al. proposed an objective function for parameter estimation in which the similarity of the attractor's topologies was employed [40]. Although we do not seek a feasible solution to an optimization problem in this paper, the same approach can provide us with the proper means of similarity estimation, which is our case study. Applying the method proposed by Rybin et al. [41] may also provide improved results, which is suggested to be considered for future works.
To apply the method, return maps of systems before and after coupling implementation are constructed. A return map (or Poincaré map) is an ordered pattern of a chaotic flow obtained by recording local maxima in a given run time. Local maxima samples constitute a time series P(n). The return map is a plot of P(n + 1) as a function of P(n). Figure 3 shows the ordered patterns of return maps obtained from original cases. In Figure 3a, all three return maps of Lorenz states are depicted. Any given initial condition in the basin of attraction of the same attractor results in a similar pattern. Return maps of Rössler and Hindmarsh-Rose are illustrated in Figure 3b,c, respectively.

Results
As previously mentioned, the network of two non-identical systems is considered, and the coupling strength varies. The introduced cost function is calculated for each system, and its similarity with the original attractor is found. The fourth-order Runge-Kutta method is used to provide the numerical solutions of differential equations. The required computational time of each coupling case varies according to the speed at which they evolve. These intervals, concerning the integrational steps of 0.005, are long enough to construct return maps with sufficient local maxima samples. They also vary when systems are scaled in frequency. Consequently, the used computational time for each case is mentioned in the corresponding section. Moreover, 40% of data are omitted to avoid any transient. Initial conditions are taken from the original trajectories.
As the coupling strength of the network increases, the interaction of oscillators changes the inherent dynamics of the systems. Therefore, as well as the quantitative changes in each oscillator's fluctuations, qualitative alternations may occur in the system dynamics, such as bifurcations or unboundedness. Under such conditions, the value and trend of the cost function are no longer valid to determine the overcoming attractor (as Then, we take the cost function of each state to be the average of all calculated minimum distances over the whole dataset: Note that in each step, we have three datasets for each state x, y, and z. Thus, the final value of cost for each system is as follows: Cost =Cost x × Cost y × Cost z Of course, a fair number of local maxima from original attractors is required, so we must have enough data in hand to plot return maps. This provides us with maximally covered patterns, and the possible error arising from the paucity of data is minimized. Here, an average of 5000 samples is taken to plot original patterns. The number of local maxima in coupled systems need not be the same as the original data.

Results
As previously mentioned, the network of two non-identical systems is considered, and the coupling strength varies. The introduced cost function is calculated for each system, and its similarity with the original attractor is found. The fourth-order Runge-Kutta method is used to provide the numerical solutions of differential equations. The required computational time of each coupling case varies according to the speed at which they evolve. These intervals, concerning the integrational steps of 0.005, are long enough to construct return maps with sufficient local maxima samples. They also vary when systems are scaled in frequency. Consequently, the used computational time for each case is mentioned in the corresponding section. Moreover, 40% of data are omitted to avoid any transient. Initial conditions are taken from the original trajectories.
As the coupling strength of the network increases, the interaction of oscillators changes the inherent dynamics of the systems. Therefore, as well as the quantitative changes in each oscillator's fluctuations, qualitative alternations may occur in the system dynamics, such as bifurcations or unboundedness. Under such conditions, the value and trend of the cost function are no longer valid to determine the overcoming attractor (as the former attractor no longer exists). Given the aim of the study, such alternations in the quality of system dynamics are undesirable. Thus, parameter d is taken within a range in which these conditions are satisfied.
The following results are presented for three different coupling configurations of Lorenz, Rössler, and HR systems. For each coupling configuration, three cases are considered: For each case, the cost functions are computed. The frequency and amplitude scaling results are shown in blue, red, and green colors.

Lorenz-Rössler Coupling
The network of coupled Lorenz and Rössler systems is described in Equation (8).
For these systems, the coupling coefficient is raised to d = 0.1 as the analyses have shown that no qualitative change occurs in this range. As shown in Figure 1, the Rössler fluctuations are less frequent than the Lorenz system. Therefore, the frequency scaling is applied to make Rössler as fast as Lorenz. Thus, Equation (8) is transformed into Equation (9): where ω R = 5.53. To apply the amplitude scaling, two control parameters are required, which are the absolute value of extrema (maxima or minima) of the x state in each system. Taking the x amplitudes of Lorenz (C L ) and Rössler (C R ) as control parameters, Equation (8) turns into Equation (10): where C L = 17.96 and C R = 17.49. Both frequency and amplitude scaling are performed by combining Equations (9) and (10). To avoid repetition, the resulting equation is not presented. The equations are solved for 7000 time units in general form and for 5000 time units after frequency scaling. Figure 4 represents the cost functions concerning the coupling strength. Lorenz and Rössler's costs are depicted in Figure 4a,b, respectively. The blue, red, and green cost functions refer to the case of frequency scaling, amplitude scaling, and both, respectively. The coupling strength is raised to d = 0.1 to observe the evident growth in cost trends. As it is shown in Figure 4, for d = 0 where no coupling is between oscillators, the cost is also equal to zero as each attractor is precisely the same as the original one. By increasing d, a gradual rise appears in the value of the costs. It can be observed that for all three cases, the Rössler cost range is 10,000 times higher than the Lorenz range. To be specific, in the absence of amplitude scaling, the Lorenz cost range is on the order of 10 −5 and the Rössler cost is within the range of 10 −1 . As amplitudes are scaled, the former ranges convert to 10 −6 and 10 −2 , respectively. Thus, regardless of the difference in frequency and attractor size, the cost ranges of Lorenz are much smaller than the cost ranges in Rössler. This implies that Rössler's dynamics do not succeed in affecting Lorenz's dynamics as much as Lorenz impairs the Rössler attractor. Therefore, it can be concluded that the Lorenz attractor overcomes the Rössler attractor in the given range of coupling strength.
For visualization purposes, return maps of Lorenz and Rössler systems when the coupling strength is set to d = 0.05 are shown in Figure 5a,b, respectively. The left, middle, and right columns refer to the frequency scaling, amplitude scaling, and both. Return maps of uncoupled systems are plotted with light green in the background to compare original and coupled patterns. As shown in Figure 5a, the Lorenz attractor almost remains unchanged in all cases. On the other hand, Rössler samples are less ordered and the coupling results in a different and scattered return map compared with its original pattern (Figure 4b).  For visualization purposes, return maps of Lorenz and Rössler systems when the coupling strength is set to d = 0.05 are shown in Figures 5a,b, respectively. The left, middle, and right columns refer to the frequency scaling, amplitude scaling, and both. Return maps of uncoupled systems are plotted with light green in the background to compare original and coupled patterns. As shown in Figure 5a, the Lorenz attractor almost remains unchanged in all cases. On the other hand, Rössler samples are less ordered and the coupling results in a different and scattered return map compared with its original pattern (Figure 4b). Light green plots represent uncoupled systems return maps. The blue, red, and dark green dots correspond to frequency scaling, x amplitude scaling, and both, respectively.

Lorenz-Hindmarsh-Rose Coupling
The network of coupled Lorenz and Hindmarsh-Rose systems is described in Equation (11). Light green plots represent uncoupled systems return maps. The blue, red, and dark green dots correspond to frequency scaling, x amplitude scaling, and both, respectively.

Lorenz-Hindmarsh-Rose Coupling
The network of coupled Lorenz and Hindmarsh-Rose systems is described in Equation (11).
In this case, the coupling strength varies within the range [0,0.1] in which the attractors are only affected quantitatively. As previously mentioned in Section 2.1, Hindmarsh-Rose evolves much slower than Lorenz. Thus, in this coupling case, the frequency of the HR system is scaled. Under this transformation, Equation (11) changes to Equation (12): where ω H = 31.25.
To study the role of size of attractors in the results, x state amplitudes are scaled by two control parameters C L and C H that are Lorenz and HR scaling parameters, respectively. Thus, Equation (11) is transformed into Equation (13): where C L = 17.96 and C H =1.8. To achieve the equations of the Lorenz-HR network after both frequency and amplitude scaling, it suffices to combine Equations (12) and (13). To avoid repetition, the resulting equation is not represented. Due to the slow dynamics of HR, simulations are performed up to 50,000 time units. When the frequency of HR is scaled, a computational time of 7000 suffices.
The cost functions concerning coupling coefficient d are provided in Figure 6. Lorenz and HR costs are represented in Figure 6a,b, respectively. Similar to the previous case, the coupling strength is elevated until d = 0.1. We observe that the value of costs when no coupling is applied (d = 0) is zero, and an incremental trend is seen as coupling strength increases. Furthermore, in all cases, the Lorenz cost range is much smaller than the HR cost range, implying that despite the type of applied scaling, the overcoming attractor is Lorenz. The HR cost ranges tend to be 10,000 higher than those of Lorenz. When both scalings are applied, the HR cost slightly improves and lowers from 0.1575 and 0.22 down to 0.09. On the other hand, the Lorenz cost is also enhanced, and the ratio of HR and Lorenz costs remains the same. To investigate the abovementioned difference in HR costs, return maps of Lorenz and HR systems, when the coupling strength is d = 0.01, are provided in Figure 7. For comparison, original patterns are also shown by light green in the background. Visual evidence in Figure 7a confirms the low value of Lorenz cost as the coupled Lorenz has its former return map patterns in all cases. HR return maps are illustrated in Figure 7b, indicating that coupled HR return map patterns, notably in x and y states, are thoroughly impaired by the emergence of scattered points. In general, distributed points of coupled HR seem to construct a rough order. However, this order differs when no frequency scaling is performed . To be more accurate, in the red graph, a different line of points that is nowhere near the original pattern appears, and fewer points exist near the original pattern of the HR x state. This reasoning seems to justify the previously mentioned cost difference. The trend of HR cost in the case of x amplitude scaling (red plot) differs considerably from the HR cost in two other cases. This dissimilarity is well recognizable at d = 0.01 as cost values in blue and green graphs (when frequency scaling is applied, whether x is scaled or not) are approximately equal to zero. In contrast, the cost does not equate to zero in the red plot (where HR frequency is not scaled). Although the dominant attractor is the Lorenz attractor, it can be deemed that slow HR is more prone to be impaired by Lorenz.
To investigate the abovementioned difference in HR costs, return maps of Lorenz and HR systems, when the coupling strength is d = 0.01, are provided in Figure 7. For comparison, original patterns are also shown by light green in the background. Visual evidence in Figure 7a confirms the low value of Lorenz cost as the coupled Lorenz has its former return map patterns in all cases. HR return maps are illustrated in Figure 7b, indicating that coupled HR return map patterns, notably in x and y states, are thoroughly impaired by the emergence of scattered points. In general, distributed points of coupled HR seem to construct a rough order. However, this order differs when no frequency scaling is performed. To be more accurate, in the red graph, a different line of points that is nowhere near the original pattern appears, and fewer points exist near the original pattern of the HR x state. This reasoning seems to justify the previously mentioned cost difference. Unlike the results obtained from Lorenz-Rössler coupling, Lorenz-HR coupling indicates that, when non-identical oscillators vary notably in the time scale, matching frequencies might somewhat help the slower system (in this case, HR) to preserve its Unlike the results obtained from Lorenz-Rössler coupling, Lorenz-HR coupling indicates that, when non-identical oscillators vary notably in the time scale, matching frequencies might somewhat help the slower system (in this case, HR) to preserve its attractor. Note that the frequency-adjusting parameter of HR in Lorenz-HR coupling is ω H = 31.25, whereas in Lorenz-Rössler coupling, it is ω R = 5.42.

Rössler-Hindmarsh-Rose Coupling
The equations of our last case of study, the coupling of Rössler and Hindmarsh-Rose, are as follows: In this case, the slower oscillator is HR. Thus, frequency scaling is applied to HR equations, and the frequency-adjusting parameter is taken as ω HR to avoid repetition. With this transformation, Equation (14) changes to Equation (15): To scale x amplitudes, scaling parameters C R and C H turn Equation (14) to Equation (16): where x scaling parameters are set to C R = 17.49 and C H = 1.8.
Moreover, to evaluate the roles of time and amplitude scale simultaneously, Equations (15) and (16) are combined. As both coupled systems are slow, a computational time of more than 60,000 is needed. Scaling HR frequency decreases the simulation time to 50,000.
The cost results of this coupling configuration are depicted in Figure 8. In this case, it suffices to raise the coupling strength to d = 0.01, as both attractors become vulnerable as the coupling becomes stronger. Figure 8a,b stand for Rössler and HR costs, respectively. Again, the trends of costs follow the expected pattern, and they grow by increasing the coupling strength. According to the blue plot, Rössler and HR cost orders are within the same range when the HR frequency is scaled to be as fast as Rössler. Although the HR cost is 0.0013 and Rössler owns a cost of 0.0018, this slight difference can be ignored when return maps are analyzed further (both are on the order of 10 −4 ). However, the Rössler attractor tends to overcome HR when x amplitudes are also scaled (as its cost becomes 10 to 100 times lower than HR cost). Compared to Lorenz-Rössler and Lorenz-HR couplings, cost orders are much closer in this case. To shed light on the dominant attractor, return maps concerning d = 0.005 are provided in Figure 9. Comparison of Figures 9a,b, which are the corresponding return maps of Rössler and HR, respectively, implies that the Rössler attractor remains unaltered. Even though Figure 9b shows that adjacent points of the coupled HR return map exist near the original plot (light green in the background), which is the root of its low cost, its pattern changes. Hence, in this case, the Rössler attractor is the dominant attractor. To shed light on the dominant attractor, return maps concerning d = 0.005 are provided in Figure 9. Comparison of Figure 9a,b, which are the corresponding return maps of Rössler and HR, respectively, implies that the Rössler attractor remains unaltered. Even though Figure 9b shows that adjacent points of the coupled HR return map exist near the original plot (light green in the background), which is the root of its low cost, its pattern changes. Hence, in this case, the Rössler attractor is the dominant attractor.

Lyapunov Exponents of Systems
To compare the inherent characteristics of our studied systems, the largest Lyapunov exponents (LLE) are computed in all cases and presented in Table 1. As can be seen, the

Lyapunov Exponents of Systems
To compare the inherent characteristics of our studied systems, the largest Lyapunov exponents (LLE) are computed in all cases and presented in Table 1. As can be seen, the most significant value of LLE belongs to the Lorenz system. An important observation is that when the size of an attractor is scaled, it has no impact on LE values. Thus, LLE values of amplitude-scaled systems are not presented, to avoid repetition. However, time scaling affects the LEs of a system directly. For instance, consider the Lorenz-Rössler coupling configuration. When the control parameter scales the Rössler frequency (ω R ), its LLE is also multiplied by ω R (turning 0.08 to 0.433). As the results in previous sections suggest, in Lorenz-Rössler and Lorenz-HR cases, the Lorenz attractor is dominant. In both cases, even by scaling Rössler and HR timing, the Lorenz LLE is twice the LLEs of Rössler and HR. As the cost data in Figure 8 suggest, when the HR frequency is scaled in Rössler-HR coupling, Rössler barely overcomes HR. There is speculation that the reason for such a very close competition lies in the fact that the LLEs of the two attractors are so close (Rössler LLE = 0.08 and time-scaled HR LLE = 0.075). Consequently, one can generalize that in all cases, the attractor with a greater value of LLE is the dominant one. Kaplan-York dimensions (D KY ) of attractors are also calculated using Lyapunov exponents and are shown in Table 1. Any type of scaling, whether amplitude or frequency, fails to alter the D KY of attractors.

Conclusions
This study investigated the dominant attractor in Lorenz-Rössler, Lorenz-HR, and Rössler-HR coupled systems. The similarity of each system's former and latter attractor was assessed by calculating a cost function based on Euclidean distances of return map samples. The impact of attractor size and system oscillation frequency was also considered. Relative comparisons were carried out between cost ranges in a given period of coupling strength. In the Lorenz-Rössler case, the Lorenz cost range was 10,000 times lower than that of Rössler regardless of scaling in frequency and amplitude. The ratio of costs in Lorenz-HR coupling resulted in the same number of 10,000. In the coupling case of Rössler-HR, costs tended to be within the same range, and the most significant disparities in ratios were within the range of 10 to 100. Altogether, we surmised that an attractor with a lower cost range was the one that succeeds in maintaining its previous dynamics. An illustration of return maps confirmed our speculation based on cost results. It can be concluded that the rerun map of an attractor with a lower cost range tends to preserve its previous patterns. It was obtained that in each case, the attractor with a higher value of LLE dominated the other. Lorenz, with the greatest LLE value of 0.819, held the lowest cost values and the least impairments in the return map pattern. Our observations also implied that in Hindmarsh-Rose coupling cases, scattered points appeared and impaired its return map pattern. The effect of frequency difference was highlighted by the fact that when the HR frequency was scaled, some scattered points were omitted, and its cost was enhanced. The LLE value of HR also supports this outcome as a slow system. When frequency scaling was applied to HR, its LLE was multiplied by the scaling parameter. It raised from 0.013 to 0.075, while Rössler held an LLE value of 0.08. In this case, the closeness in LLE was in line with the close range of Rössler and HR costs. Thus, the raised value of LLE helped HR preserve its attractor. Finally, our results suggest that an attractor size has no role in determining the overcoming attractor, and it is supported by the fact that it does not change LLE values either.