Limits to Detection of Generalized Synchronization in Delay-coupled Chaotic Oscillators

We study how reliably generalized synchronization can be detected and characterized from time-series analysis. To that end, we analyze synchronization in a generalized sense of delay-coupled chaotic oscillators in unidirectional ring configurations. The generalized synchronization condition can be verified via the auxiliary system approach; however, in practice, this might not always be possible. Therefore, in this study, widely used indicators to directly quantify generalized and phase synchronization from noise-free time series of two oscillators are employed complementarily to the auxiliary system approach. In our analysis, none of the indices provide the consistent results of the auxiliary system approach. Our findings indicate that it is a major challenge to directly detect synchronization in a generalized sense between two oscillators that are connected via a chain of other oscillators, even if the oscillators are identical. This has major consequences for the interpretation of the dynamics of coupled systems and applications thereof.


I. INTRODUCTION
Over the past few decades scientists from different disciplines have been trying to identify and quantify synchronization in coupled dynamical systems and a large number of measures have been proposed [1,2]. These studies are based on the understanding and knowledge of the complex dynamics in coupled systems. Most of these analyses focused on the study of identical synchronization, in which coupled oscillators follow the same trajectory [3,4]. Even if their trajectories do not present a one-to-one correspondence, there still remains the possibility for the systems of being synchronized in a generalized sense. Generalized synchronization (GS) occurs when one of the oscillators completely determines the dynamical evolution of the other oscillator. Therefore, identical synchronization can be considered as a particular case of GS. Another type of synchronization is phase synchronization (PS), in which phases of the oscillators are correlated, while their amplitudes might be uncorrelated.
The synchronous behavior in coupled dynamical systems has been extensively analyzed in many types of networks, for instance, in biological systems [5], neural networks [6], and lasers [7]. The way the elements interact makes a strong impact on their collective dynamical behaviors. A key factor for the synchronization behavior is the network connectivity in which the dynamical elements are interacting with each other. The influence of the underlying network on the synchronization properties in coupled dynamical systems has been intensively analyzed [8][9][10][11][12][13][14]. In many of these studies, correlation-based measures were often used as an indicator for synchronization. However, when the elements synchronize in a generalized sense, their synchronization can be overlooked if one analyzes it only with the correlation-based measures. One typical example occurs in a drive-response configuration network [14][15][16][17][18] with the auxiliary system approach proposed in Ref. [15]. In these studies, the oscillators at the end of the original driveresponse configuration and the auxiliary system converge to an identical synchronized state after some transients, although the oscillators are set to different initial conditions and exhibit sensitive dependency on the initial condition when they were uncoupled. Identical synchronization between the response oscillator and the auxiliary system arises under the condition that their initial conditions are set to the same basin of attraction [15,19]. The drive-response system is, therefore, considered to be synchronized in a generalized sense. In recent studies [17,18], correlation and mutual information measures between the drive and response systems were found to decrease if an increasing number of elements were placed between the drive and response elements, in a scheme as shown in Fig. 1 (or the more general ones depicted in Fig. 1 of Ref. [18]). However, the auxiliary system approach indicated that synchronization between the elements n/2 and n /2 remains. As in the simple drive-response case, it can be concluded that X and Y (or Z) exhibit synchronization in a generalized sense.
In these previous studies, synchronization was identified using the auxiliary system. By using the auxiliary system approach, we have also tested in this paper that the time series generated by the coupled systems depicted in Fig. 1 are synchronized in a generalized sense. Since synchronization in a generalized sense exists between X and Y, we expected that some of the recently proposed indices can detect this synchronization. Therefore, the aim of this study is to evaluate whether some indices have the ability to identify synchronization in a generalized sense between elements 0 and n/2 in Fig. 1 and to understand the limits of these indices.

II. COUPLED MACKEY-GLASS SYSTEM
The system we consider in this study is a ring configuration consisting of n identical nonlinear oscillators that are delay coupled, as depicted in Fig. 1. We chose the Mackey-Glass oscillator as the nonlinear elements. The normalized equation Mackey-Glass oscillators are chosen as network nodes and arranged in a ring. The oscillators in the ring are delay coupled, with a total delay in the loop of τ d . We focus on elements X, Y, and Z. The part surrounded by gray dashed line is an optional chain of coupled oscillators considered in Refs. [17,18] and equivalent to oscillators 1 and n/2 on the ring configuration.
of the j th oscillator in the system is given by with and The total delay τ d in the loop was chosen as 30. Throughout this paper, GS and PS between the oscillators X and Y in Fig. 1, which correspond respectively to oscillators 0 (X) and n/2 (Y), were quantified by using time series consisting of 20 000 data points (sampling rate t = 0.1). All data for GS and PS evaluations in this paper were obtained from noise-free numerical simulations in order to explore the fundamental limitations of the different methods. Previous studies of the ring configuration showed that these two oscillators do not exhibit identical synchronization within the chaotic regime [17], while in-phase synchronization and rotating waves can be observed eventually in the case of periodic dynamics [20,21]. The parameters in Eqs. (1)-(3) are the same as in Ref. [18] and set to a = 2.1, b = 1/3, and c = 9.9. The characteristic time constant is 0.47. In Fig. 2 we show amplitude time traces of elements X and Y for an increasing number of elements n. The oscillators exhibit chaotic fluctuations, independently of the number of oscillators in the system, due to the delayed coupling. Because of these fluctuations the Mackey-Glass oscillators cover a wide frequency range, as shown in Fig. 3. Distinct peaks can be seen in the spectra, e.g., in the case of n = 2. Such peaks can be distinguished when the number of elements n is less than 8. For a larger number of elements the spectra become flat. The cross correlation (CC) functions between the oscillators X and Y are shown in Fig. 4. Since the total delay in the system is 30 in the ring configuration, the delay between the oscillators X and Y is 15. For this reason, the peaks of the CC functions appear at time shifts of T s = ±15. In agreement with Refs. [17,18], the amplitudes of the largest peaks gradually decrease as the number of oscillators in the system increases until the peaks practically disappear. The sign of the largest CC coefficient changes depending on the number of oscillators in the system, indicating that it determines the dominance of in-phase or antiphase synchronization between X and Y.

III. GENERALIZED SYNCHRONIZATION
Our aim at this point is to check whether different indicators commonly used in the literature are able to detect GS when the linear methods fail to do so. In this section we introduce indices based on correlation integrals, namely cross redundancies (CRCI) and joint probability of recurrences (JPR) [22,23], and the GS from the reconstructed state space of X and Y, namely the generalized synchronization index (GSI) and the convergent cross mapping (CCM) [24,25].

Cross redundancies
Information theory can be used to unveil relationships between dynamical systems [22,26]. One of the most prominent measures is mutual information. The mutual information for a multivariate time series is also called redundancies [27]. For estimating local probability distribution functions, a fixed volume approach is used in redundancies based on correlation integrals [27]. The use of the correlation integral to estimate the correlation dimension of strange attractors was first introduced in Ref. [28]. Let us consider a time series {x(k)} (k = 1,2, . . . ,L). When applying Takens embedding, the kth point on the multidimensional state can be written as x(k) = (x(k),x(k − τ e ),x(k − 2τ e )), . . . ,x(k − (m − 1)τ e ), where m and τ e are embedding dimension and embedding delay, respectively [29]. For this multidimensional state space  (4) where N = L − (m − 1)τ e is the number of embedding vectors, (·) is the Heaviside step function, and ε is a threshold.

X, the correlation integral is defined by
To quantify the dependence of x m on x 1 conditioned on x 2 , . . . ,x m−1 we compute the conditional redundancy as Based on this equation, the conditional redundancies were extended to multiple variables [27]. Assuming that H (X) ≈ − log C[X], Prichard and Theiler [27] proposed the timelagged mutual information measure between two time series {x(k)} and {y(k)} (k = 1,2, . . . ,L) as by replacing x 1 , x 2 , and x 3 of Eq. (5) in a case of m = 3 with x(k), τ s y(k), and y(k), respectively. τ s y(k) denotes the difference of y(k) and y(k + τ s ).

Joint probability of recurrences
The cross-redundancy correlation index introduced in the previous section is a measure in which the correlation integrals are applied to the conditional mutual information. Romano et al. [23] proposed another application of the correlation integrals for GS identification. The main idea in their indicator is to quantify the coincidence probability of recurrences between multidimensional states X and Y . For this, they modified the threshold in the definition of the correlation integral [Eq. (4)] and proposed to compute the coincidence probability of recurrences as follows. For the multidimensional state X, a constant value N N is introduced. This constant value satisfies N l=1 (ε k x − x(k) − x(l) ) = N N for ∀k where the threshold ε k x is chosen such that the number of nearest neighbors for each column becomes N N (analogously for Y ). For the multidimensional states X and Y , using the number of shifted points τ s , the degree of coincidences is quantified by a coefficient S XY (τ s ), which is defined as where R R = N N /N is a free parameter, which lies between 0 and 1. The joint probability of recurrences is calculated from the normalized maximum value of S XY (τ s ) as This indicator is robust for a wide range of parameters [23].

Generalized synchronization index
In the previous two GS indicators, we used the correlation integral to compute their values. Now we introduce a GS indicator that does not use the correlation integral but the distance information of the geometry of multidimensional states. For each vector x(k), let us assume that r k [i] (i = 1,2, . . . ,R) stores the time indices of R points that are geometrically nearest to x(k) in the multidimensional states. The neighborhood is formed by the basis of Euclidean distances and the R closest neighbor points of x(k) that have smaller Euclidean distances in the multidimensional state. If we focus on the R neighboring points between x(k) and x(r k [i]) for the R points, the average Euclidean distance is given by [24] As it happens for the state space X, variable s k [i] stores the time indices of R nearest-neighbor points of y(k); therefore, the average Euclidean distance with the R mutual neighbors for the bivariate series is defined by If two elements generating the time series, {x(k)} and for independent elements, indicating that mutual neighbors are more spread.
From the definitions, D R k (X) refers to the self time indices (r k [i]), whereas D R k (X|Y ) refers to those of another multidimensional state vector (s k [i]). Therefore, the distances calculated inside of the sum of Eq. (9) must be smaller than those of Eq. (10). For this reason, the value D R k (X|Y ) is always larger than D R k (X). If the observed time series are noisy and/or the interactions between the two elements are noisy, the mutual neighbors spread further in the multidimensional state space. In such a case, D R k (X|Y ) further increases with respect to D R k (X) and it becomes difficult to detect the synchronization between the two elements. In order to emphasize the contribution of Y to X, Bhattacharya et al. [24] diminished the self-influence, namely the original neighbors defined in Eq. (9). To diminish its influence, they proposed to use more distant points in the multidimensional state. Let us redefine r k [i] by assuming that it stores the time indices of 2R (i = 1, . . . ,2R) nearest points in the multidimensional state space X.
Combining Eqs. (10) and (11), the strength of influence of Y on X is defined as This index is a natural extension described by Eqs. (9) and (10), which is called the similarity index [24]. In an analogous way, the dependence of X on Y, E(X|Y ), can be also estimated as follows: where D k (Y ) and D R k (Y |X) are defined as follows: where s k [i] also stores 2R nearest-neighbor points in the multidimensional state space Y .

Convergent cross mapping
The GSI directly quantifies GS between two oscillators with the geometrical information of the multidimensional state spaces. As another approach to detect GS, Sugihara et al. [25] proposed a method to predict time series with similar multidimensional state information as the GSI. The state spaces X and Y , which are embedded with dimension m and embedding delay τ e , are prepared to predict the original time series {x(k)} and {y(k)}. Here we describe the predicted time series as {x(k)} and {ŷ(k)}. When redefining the variables r k [i] and s k [i] that store the m + 1 nearest-neighbor indices of the kth point on the multidimensional state spaces, x(k) is predicted from x(s k [i]) (i = 1,2, . . . ,m + 1) by a locally weighted average aŝ where w k [i] is a weight based on the distance between x(k) and x(s k [i]). The weights are determined by the following equation: where In a similar way, y(k) can be predicted aŝ where Synchronization is estimated by computing the crosscorrelation coefficient between the original time series and the predicted one.

IV. PHASE SYNCHRONIZATION
Complementary to the GS indicators, in this section, we define quantifiers to detect PS. In addition, we present a brief introduction to the Hilbert transform, which is often used to extract phase information from time series of real variables.

A. Hilbert transform
Before being able to compute some phase synchronization indicators, it is necessary to estimate phases. Among the three basic approaches to estimate phases of time series, Hilberttransform-based phase estimation [22], wavelet-transformbased phase estimation [30], and event-timing-based phase estimation [26], we adopt the method based on the Hilbert transform, because it is mathematically well defined for any time series and does not require stationarity. Although, in a strict sense, it does only provide meaningful phase information for phase-coherent attractors, we apply it in the following without explicitly testing this property. For a time series {x(k)}, the analytic signal is where x H (k) is the Hilbert transform of x(k) and P.V. stands for the Cauchy-Principal value. The analytic signal can be decomposed further as where a(k) and φ(k) are the instantaneous amplitude and phase, respectively. Using the phase series of φ x (k) and φ y (k) computed from Eq. (25), the phase difference between two signals is simply defined as The instantaneous phases φ(k) computed from Eq. (25) are real and not constrained to the interval zero and 2π . In fact, several bound intervals can be considered, but here we define and restrict (k) to [0,2π ) using the simple modulo 2π operation.

B. Mean phase coherence
One approach to PS is based on phase locking. Phase locking requires the rigorous condition that relative phases or phase differences have a constant value. However, this condition is weakened when estimating the existence of PS between noisy or chaotic oscillators. In this case, the phase-locking condition is fulfilled if phase differences remain bounded. The difference must be smaller than a certain value that is bounded to 2π . From the relative phases, as estimated above, the circular variance of the angular distribution is obtained by transforming the relative phase angles on the unit circle in the complex plane. With this circular variance , Mormann et al. [31] proposed mean phase coherence (MPC) as a measure for the phase synchronization, using the index based on conditional probability proposed in Ref. [32]. The MPC is given by where L is the data length [31]. The index MPC ranges from zero, which corresponds to a nonsynchronous state, to 1, which indicates that the two oscillators are perfectly synchronized in phase. 062924-5

C. Phase lag index
In the MPC, instantaneous relative phases are projected onto the circular unit and the length of the average resultant vector is computed to quantify PS. Another approach to quantify PS is to define the asymmetry of the relative phase distribution. When two oscillators do not achieve PS, this distribution is expected to be flat. In turn, it is considered that the deviation from a flat distribution can be attributed to PS between two oscillators. The phase lag index (PLI) measures this asymmetry [33]. The asymmetry of the phase difference indicates the likelihood of the phase differences in the interval (−π,0) and the interval (0,π ) and is quantified by The definition in Eq. (27) demands for the phase differences to be in the interval of −π < (k) π . However, due to our definition of the phase difference via the Hilbert transform, it lies within the range of 0 (k) < 2π . Therefore, Eq. (27) is modified as

D. Conditional mutual information of phase
For the identification of the relation between elements in a system, measures derived from information theory can be effective [22,26].
Let us consider two variables X and Y with marginal probability distributions p(x) and p(y). From the probability distribution p(x), the Shannon entropy of X is defined as With the joint probability distribution function p(x,y), the joint entropy of X and Y is computed by Using these entropies, the mutual information between X and Y is defined as The conditional mutual information (CMI), which is the mutual information between X and Y conditioned by a third variable Z, can be written as Paluš and Stefanovska [34] introduced the time-shifted time series of the second component in Eq. (32), so the CMI is written, in this case, as I (φ x (k); φ y (k + τ s )|φ y (k)). Paluš and Vejmelka [26] showed that it is possible to identify the direction of the connections in the coupled Rössler system with I (φ x (k); φ y (k + τ s )|φ y (k)) and I (φ y (k); φ x (k + τ s )|φ x (k)). However, they also revealed that the use of τ s φ x,y (k), instead of φ x,y (k + τ s ), gives a better identification of the directionality. For this reason we also use, in this study, the definition of τ s φ x,y (k) as the CMI of phase.

E. Phase synchrony index based on entropy
In the CMI introduced above, PS is quantified via information theory. Another application of information theory is to directly use relative phases or phase differences. In the next PS indicator, relative phases or phase differences are applied to entropy. The instantaneous phases of individual time series are computed by Eq. (25), and the relative phases are defined as their differences as explained above. These relative phases have been projected onto the unit circle. The restriction of (k) being within the interval of [0,2π ) is to eliminate the effect of noise-induced phase slips when the signals are noisy. The distribution of (k) is estimated with the optimal bin number of M = e 0.626 + 0.4 ln(L − 1) [35]. In order to quantify the strength of phase synchronization, the following index is computed [32]: where H max = ln(M) is the maximum entropy that is obtained when the probabilities of all the bins are equivalent. This index measures the deviation of the relative phase distribution from the uniform distribution, the entropy of which is precisely H max = ln(M).

F. Correlation of probability of recurrences
We have already introduced four indicators for PS. So far, all the indicators use instantaneous phases estimated by the Hilbert transform. Here we will introduce a PS indicator that does not use instantaneous phases computed from the Hilbert transform. Romano et al. [23] proposed the use of recurrences of chaotic systems to detect PS indirectly since the diagonal lines of a recurrence plot show that some determinism exists in a system. Additionally, the characteristic time scale of the system is reflected in the vertical distances between these diagonal lines. The distances become constant for periodic systems but several distances can be found for chaotic systems, reflecting that chaotic systems present several time scales [2].
If two oscillators are in PS, the patterns of their recurrence plots coincidently appear because their phases are adapted with each other, even though amplitudes of their oscillators are uncorrelated. Namely, if a high probability is observed in the recurrence plot of one oscillator at τ s steps, the recurrence plot of the other oscillator also exhibits a high probability at τ s steps. The ratio of τ s recurrence is directly estimated from the recurrence plot as where τ s is time shift. This probability is called the generalized auto correlation. Correlation of probability of recurrences (CPR) is the CC between P X (τ s ) and P Y (τ s ), which quantifies their PS and is given by where P X,Y (τ s ) indicates that the mean value has been subtracted; σ X,Y are the standard deviation of P X,Y (τ s ).

V. RESULTS
In this section, we evaluate the performance of the indicators described in the previous sections by applying them to a system of n delay-coupled Mackey-Glass oscillators unidirectionally connected in a ring configuration. We particularly discuss the influence of an increasing number of elements in the considered system for each indicator.

A. Evaluation of generalized synchronization indices
First, we evaluate GS indices from the cross redundancies based on correlation integrals (CRCI). Figure 5 shows the CRCI as a function of the time shift T s . I (x(k); τ s y(k)|y(k)) for n = 2 exhibits prominent peaks separated by the total delay  τ d = 30. This function reflects the symmetry of the system and also the second peaks can be seen in it. For n = 6, the largest peaks are still distinguishable but the second ones have disappeared under the baseline. The amplitude of the largest peaks diminishes with an increasing number of oscillators in the system and disappear when n is larger than 10 (upper panels for n = 8 and 18 in Fig. 5). All properties are analogous to those of the CC function (Fig. 4) except for a change of the sign of the peaks.
To illustrate the influence of an increasing number of oscillators n, we plot the maximum value of the CRCI vs n in Fig. 6. I (x(k); τ s y(k)|y(k)) exponentially decreases for an increasing number of oscillators. The value of I (x(k); τ s y(k)|y(k)) drastically decays from n = 2 to 4. A clear difference between these cases can be identified depending on whether two oscillators interact directly or indirectly, suggesting that mediating oscillators make GS detection much harder. A further increase in the number of oscillators leads to a decrease of I (x(k); τ s y(k)|y(k)) resulting in a value very close to zero for n = 10. In fact, the dependence of I (x(k); τ s y(k)|y(k)) with n saturates for n > 8. The fact that the maximum amplitude of the peaks decreases for increasing n is a property shared with the CC (see Ref. [18]); however, I (x(k); τ s y(k)|y(k)) monotonically decreases since the CC globally decreases but oscillating for an increasing number of oscillators. The dependence of I (y(k); τ s x(k)|x(k)) is a mirror image of the result of I (x(k); τ s y(k)|y(k)).
The second indicator is the joint probability of recurrences (JPR). This measure is the normalized maximum of S XY (τ s ), as detailed in Eqs. (7) and (8). For comparison, we plot in Fig. 7 the latter index as a function of the time shift and in Fig. 8 the peak value of JPR as a function of the number of intermediate oscillators, n. Clearly, the profile of S XY (τ s ) mimics that of the CRCI (cf. Fig. 5 vs Fig. 7). In particular, the peaks are located at precisely the same time shifts for both indices. An apparent difference between them, however, is that the peaks for S XY (τ s ) are wider than for the CRCI. This result indicates that S XY (τ s ) [and therefore also JPR XY (τ s )] is less sensitive than CRCI to the time delay of the coupling between the systems.
Following the same scenario with the CRCI, we analyze how an increasing number of oscillators affects the JPR (Fig. 8). Basically the JPR has the same properties as the CRCI and the results of Fig. 8 entirely agree with those of Ref. [36] in which τ d = 300. From this fact, the performance of the JPR is independent of the total delay in the ring configuration.
The next tested indicator for GS is the generalized synchronization index (GSI). The dependence of the GSI on the time shift T s is shown in Fig. 9. Compared to the previous indicators, the GSI has a sharper peak. Due to the influence of Y on X (E(X|Y )), the peaks for positive time shift are higher than those for negative shifts. Taking into account that causal consequences appear only in the future, it is natural that peaks emerge for positive time shifts, while peaks for negative shifts are smaller. In addition, for n = 2, the third peak is visible only when the time shift T s is positive. This was not observed for the other indicators. In contrast to E(X|Y ), E(Y |X) exhibits a horizontally mirrored function (lower panels for n = 2,6, and 8 in Fig. 9), a fact that is caused by the opposite direction of the influence. These results reveal the unidirectional ring configuration, which can therefore be deduced from the dynamical properties.
As before, we also carry out the analysis as the number of elements in the ring increases and plot the maximum value of the GSI in Fig. 10, with similar results as compared with previous indicators.
The last considered method is the cross convergence mapping (CCM). The CCM is not a measure to quantify the synchronous degree of GS but a method to predict time series. As explained above, x(k) on the multidimensional state space is predicted by using the indices of m + 1 nearest-neighbor points of y(k) that are used as k-conditioned predictors. Then the correlation coefficient between the original time series ({x(k)}) and the predicted one ({x(k)}) is the indirect measure of GS degree between X and Y. Here we denote the CC coefficient between the original and the predicted time series as CC x and CC y , respectively. Figure 11 shows CC x and CC y as a function of the time shift T s . Similarly to the GSI, the asymmetry function is observed in the CC x function when n = 2 (upper panel). Namely, the peaks for the positive shifts are higher than those for the negative shifts when we predict {x(k)}. When {y(k)} is predicted, the mirror function of CC x is obtained (lower panel). For larger n (Fig. 11 n = 6,8,18), the peak heights decrease and disappear for n = 18. These properties indicate that the basic performance of CCM to detect GS is similar to other indicators. One difference is that the peaks of the CCM are not sharp and keep a high value around the coupling delay, suggesting that the CCM is more sensitive to the small changes originated from the influence of other elements. In agreement with the previous GS indicators, the dependence with an increasing number of elements of the maximum value of the CC coefficient between the original and the predicted time series exhibits a decreasing function with n (see Fig. 12). The difference is that apparently the CC coefficient is still higher than the baseline (about 0.8) when n = 10. From this point of view the CCM has the highest ability to detect GS for an increasing number of elements.

B. Evaluation of phase synchronization indices
In the previous section, we showed that almost all indicators for detecting GS fail when a ring configuration system includes more than about eight oscillators. The oscillators X and Y are, even in cases of large n, still synchronized according to the auxiliary system approach [14,17,18]. Hence, for further analysis, we introduce PS indicators here. Although we use PS indicators, our interest is this partial aspect of GS rather than claiming PS. All results for the PS indicators as a function of the number of oscillators in the ring are presented in Fig. 13.
The phase lag index (PLI) is almost zero already for the case of two mutually coupled oscillators, indicating that they do not exhibit any coherent phase relationship at phases different from 0 and π . Increasing the number of oscillators, the index remains very small, indicating that the oscillator X is asynchronous in phase with the oscillator Y.
The mean phase coherence (MPC) takes higher values compared with the PLI. When increasing the number of oscillators the MPC dependence becomes flat as in the case of the PLI. However, this index identifies some PS as the number of oscillators in the system increases. This is in disagreement with the result of the PLI and requires further analysis. To investigate this, we conducted the following process. Since we have the data of the oscillators X and Y in the cases of n = 2 to 24, we computed the MPC between X for a given n and Y for another n. These time series came from independent realizations and, thus, cannot be synchronized. We computed their average value as a baseline of non-PS. For example, to know the significance of PS in the case of n = 2, we compute the average of the MPC between X of n = 2 and Y of n = k (k = 4,6, . . . , 24). Hereafter, we call this method a data-exchange analysis. The oscillator X of n = 2 should be independent of Y of n = k because these oscillators belong to different systems. Nevertheless, our exchanged data analysis shows that the MPC yields almost the same value for all the cases. Taking into account this result, we conclude that the highest value of the MPC does not indicate PS.
The phase synchrony index ρ exhibits a similar trend as the MPC although ρ is lower than the MPC, even though both indices quantify the PS and are in the range between zero and 1. However, ρ has also a nonzero value that differs from the PLI. This inconsistency among the PS indicators occurs when the distribution of the relative phases or the phase differences is symmetric but there is a significant peak around zero or π (Fig. 14). Although their definitions differ, both the phase synchrony index (ρ) and the MPC assume a nonflat distribution of relative phases or phase differences when two elements are in PS. The PLI was also developed to quantify PS as the other PS indicators but this measure does not quantify the variance of the relative phase distribution but the asymmetry of the distribution. Therefore, even though there exists a peak in the phase distribution, the PLI will be small because the distribution is symmetric when the peak is at zero (2π ) or π . As expected, for all n, we observe the relative phase distributions whose peak is located at zero, as shown in Fig. 14. This distribution originates from the Mackey-Glass dynamics and it is not related to generalized synchronization. Next we estimate the conditional mutual information (CMI) of the phase. Although this index has direction sensitivity and we estimate both directions of the index, the CMI is always zero. This indicates that the index does not detect any PS.
As the last index to quantify PS, we calculate the correlation of probability of recurrence (CPR). This is the only index that does not use the phases estimated by the Hilbert transform. This index yields higher values when the system has less than eight oscillators but the values gradually decrease (Fig. 13).  For more than eight oscillators, this index is almost constant and yields a value about 0.6. That is, the CPR maintains at a higher value, which suggests that the oscillators are still phase synchronized, but for n 8 this value remains constant. The asymptotic value of the CPR is almost as high as in the case of the MPC, so we apply the data-exchange analysis to the CPR. As a result, the value of the CPR for the exchanged data is superimposed on the original one except for 2 n 6, indicating that the CPR can detect PS properties between X and Y for a small number of oscillators in the ring. This is a trend observed in the results of the GS indicators. As mentioned before, we assume that this difference from the other quantifiers might originate from the different definition of the phase.
To be sure that the results of the Hilbert-based PS indices were not due to the time series presenting a broad band spectra, we also tested the PS indicators with the Hilbert transform on narrow bands, band-pass-filtered signals (results not shown). Even though the filtered signals were used, we did not find any significant change in the result shown in Fig. 13.

VI. DISCUSSION AND CONCLUSION
In this study, we computed commonly applied GS and PS quantifiers in a system of coupled Mackey-Glass oscillators in unidirectional ring configurations. Our study is relevant because GS and PS can occur in various delay-coupled chaotic systems, as well as transitions among them [37]. The indices we adopted here are widely used to measure synchronization between nonlinear systems, particularly in the field of neuroscience in order to estimate functional or effective brain connectivity [3]. If brain sites interact via many intermediate neural oscillators as considered in this paper, given time series of a certain length, and with a certain amount of noise, most of the tested indices cannot identify the interactions between the sites, and the connectivity patterns derived from the indices are not reliable.
The measures adopted in this study, to assess the degree of GS between two time series, showed almost equivalent performance and similar characteristics. Particular peaks were identified at time shifts corresponding to the total delay of the coupling between the two oscillators that are most distant to each other. An interesting finding in our analyses of GS was that only the GSI and the CCM exhibit clearly different amplitudes with respect to positive and negative shifts. This asymmetry arises from the existence of time-lagged directed interactions. When the influence of Y on X is quantified, the amplitude of the peaks at positive shifts is significantly higher than that at negative shifts. The peaks at positive time shifts correspond in this case to a directed causal interaction. Due to the ring configuration, Y is also influenced by X with a time delay. Since the peaks at negative time shifts correspond to this recurrent influence of X onto itself via Y, these peaks have a smaller amplitude. Actually, the peaks at negative time shifts can be interpreted as the fingerprint of a closed feedback loop for X via Y. This might be a practical way to identify closed feedback loops since the smaller peaks would not appear if X and Y had a configuration like a drive-response system. Such connection identification is impossible for frequently used measures like the CC function or the basic mutual information function because of their symmetry. However, further analyses are needed to understand how reliable GSI and CCM are to detect a given coupling configuration.
In contrast to the GSI and the CCM, the two other quantifiers, the CRCI and JPR, exhibit symmetric functions of the time shift similarly to the CC functions. The difference of these quantifiers from the GSI and the CCM stems from what information of the constructed multidimensional state is used. The GSI and the CCM directly use the information of the multidimensional states, while for the recurrence plot (RP) analysis only a two-dimensional projection is employed. The RP is useful for visualizing the recurrences of the trajectory on the multidimensional state but loses the distance information between any pairs of points and the distance information is replaced into the simplified information: near or far. Thus, the symmetry of the function can result in the loss of distance information of the multidimensional state when computing the CRCI and the JPR. However, the distance between points of the multidimensional state plays an important role to distinguish the actual influence from mere correlation.
All the measures detected GS between the two most distant oscillators in the system with a small total number of oscillators, in particular fewer than 10 oscillators for our ring configuration. Up to this number the indicators gradually decayed with increasing number of elements and then saturated to a constant value that depends on each measure. For the quantification of GS, all the indices have similar performance and dependence on the number of elements in the ring.
In our PS analysis, we introduced five indices, four of which were estimated with phases obtained from the Hilbert transform. These four indices exhibited different degree of PS even though we analyzed the same data. In particular, the MPC and the phase synchrony based on entropy tended to better detect PS than the PLI because of the symmetry of the relative phase or the phase difference distribution, and the existence of a significant peak around zero. The MPC and the CPR exhibit relatively high values but behave differently for an increasing number of oscillators. The CPR decayed like the GS indices, whereas the MPC was constant for an increasing number of oscillators. This difference might come from the alternative definition of phase or the limitations of the Hilbert transform. In fact, all phase-related indices using the Hilbert transform exhibited no dependence on the number of oscillators. Indeed, we identified that the MPC and CPR detected pseudo-PS. The origin of this pseudo-PS would be the characteristic time scale of the Mackey-Glass oscillator or the limitations of the Hilbert transform. In addition, taking into account the baseline of PS in the CPR, the CPR also indicates weak PS when two oscillators are directly connected. Therefore, the indicators imply that the two oscillators exhibit weak or no PS.
In this manuscript, we applied previously proposed indices to GS and PS identifications to directly identify synchronization without the auxiliary system. In Refs. [17,18], independently of the number of oscillators, Y and Z in the system, as depicted in Fig. 1, exhibited identical synchronization. Therefore, it was concluded that X synchronized with Y in a generalized sense. In our approach to directly identify GS from time series of finite length, we found, however, that the more oscillators mediated between X and Y, the lower the GS indices became. Beyond a certain number of oscillators, the generalized synchronization was not detectable anymore.
Nevertheless, considering that the auxiliary system approach can still identify synchronization in a generalized sense, the information of X must be received by Y. The auxiliary system is a copy of a chain from X to Y, hence, the process of the information transfer from X to Y, namely, an oscillatormediating connection, is completely reconstructed. In contrast, when synchronization in a generalized sense is quantified with GS indicators, the oscillator-mediating process is unknown. The corresponding relationship might be complicated and embedded within a high-dimensional state space. Up to eight oscillators in a ring configuration network, the GS indicators can identify synchronization in a generalized sense, even though the connectivity between X and Y is unknown. For a large number of coupled oscillators, no quantifier captured the state of GS that is identified by the auxiliary system approach.
It remains an open question, whether for many mediating oscillators a measure can be found that it is capable to identify GS from the time series or whether this is fundamentally impossible given a time series of finite length and precision. This would also have consequence for the possible application of generalized synchronization for encryption schemes.