## 3 Methods
### Synchronization in Phase Space  
A system could be represented in the phase space, i.e., a topological representation of its behavior under different initial conditions. This method assumes that each signal represents a projection of a higher-dimensional dynamical system evolving in time, whose trajectories are embedded into a manifold, i.e., a region of its phase space. The phase space of a system could be reconstructed by means of the Takens's Theorem [38]. For a generic single temporal observation of a system \(u(t)\), the trajectory \(\mathbf{x}_{i}\) is expressed as:  
\[\mathbf{x}_{i}=(u_{i},u_{i+\tau},\ldots,u_{i+(m-1)\tau}), \tag{1}\]  
where \(m\) is the embedding dimension and \(\tau\) is the time delay. Both parameters must be properly selected to avoid redundancy in the phase space. The dimension \(m\) of the reconstructed phase space should be large enough to preserve the properties of the dynamical system: \(m\geq 2D+1\), where \(D\) is the correlation dimension of the original phase space. The correct time delay \(\tau\) should be chosen by determining when the samples of the time series are independent enough to be useful as coordinates of the time delayed vectors. For the estimation of the embedded parameters \(m\) and \(\tau\) several techniques have been proposed. As an example, the first local minimum of average mutual information algorithm [39] can be used to select the proper time delay. The minimum embedding dimension is usually estimated through the false nearest neighbors (FNN) algorithm [40]. After the reconstruction of the dynamic trajectories of the two systems in the phase space, it is possible to quantify their interacting behavior by projecting the phase space into the bidimensional cross-recurrence plots (CRPs) [41]:  
\[CR_{i,j}(\epsilon)=\Theta(\epsilon-||\mathbf{x}_{1,i}-\mathbf{x}_{2,j}||)\ \ i,j=1,\ldots,N, \tag{2}\]  
where \(\Theta\) is the Heaviside function, \(\epsilon\) is a threshold for closeness, \(N\) is the number of considered states for each system and \(||\cdot||\) the maximum norm function. A CRP represents an effective method to reduce the dimensionality of the phase space and compare the trajectories of the interacting systems as it is a matrix whose entries include information on the degree of closeness of each state of the first system with each state of the second system. A CRP exhibits characteristic patterns that show local time relationships of the segments of the trajectories of the two interacting systems. As an example, diagonal lines occur when the evolution of the states is similar at different times and their lengths are related to the periods during which the two systems move in similar ways remaining close to each other [42]. The main diagonal of a CRP represents a line of synchronization (LOS), since it implies the identity of the states of the two systems in the same time intervals. The LOS structure can be analyzed to extract information about the synchronization of the two time series [43]. In particular, the presence of LOS suggests that the two time series are fully synchronized, while discontinuities appear when the two signals are not fully synchronized. We defined the synchronization time metric (SYNC) to quantify the mean period during which the two systems are synchronized in order to reflect the dynamical synchronization behavior of the series throughout the observation period [22; 25]. SYNC is proportional to the ratio of the sum of the lengths of the subsegments \(l_{j}\) along the LOS to the total number of samples \(N\):  
\[SYNC=\frac{1}{N}\frac{\sum_{j=1}^{N_{d}}l_{j}}{N_{d}}, \tag{3}\]  
The coupling behavior in phase space can be illustrated by referring to a pair of simple sine waves.