Recurrence networks—a novel paradigm for nonlinear time series analysis

This paper presents a new approach for analysing the structural properties of time series from complex systems. Starting from the concept of recurrences in phase space, the recurrence matrix of a time series is interpreted as the adjacency matrix of an associated complex network, which links different points in time if the considered states are closely neighboured in phase space. In comparison with similar network-based techniques the new approach has important conceptual advantages, and can be considered as a unifying framework for transforming time series into complex networks that also includes other existing methods as special cases. It has been demonstrated here that there are fundamental relationships between many topological properties of recurrence networks and different nontrivial statistical properties of the phase space density of the underlying dynamical system. Hence, this novel interpretation of the recurrence matrix yields new quantitative characteristics (such as average path length, clustering coefficient, or centrality measures of the recurrence network) related to the dynamical complexity of a time series, most of which are not yet provided by other existing methods of nonlinear time series analysis.


Introduction
Since the early stages of quantitative nonlinear sciences, numerous conceptual approaches have been introduced for studying the characteristic features of dynamical systems based on observational time series [1]- [4]. Popular methods that are increasingly used in a variety of applications (see, for example, [5]) include (among others) Lyapunov exponents, fractal dimensions, symbolic discretization and measures of complexity such as entropies and quantities derived from them. All these techniques have in common that they quantify certain dynamically invariant phase space properties of the considered system based on temporally discretized realizations of individual trajectories.
As a particular concept, the basic ideas of which originated in the pioneering work of Poincaré in the late 19th century [6], the quantification of recurrence properties in phase space has recently attracted considerable interest [7]. One particular reason for this is that these recurrences can be easily visualized (and subsequently quantified in a natural way) by means of the so-called recurrence plots obtained from a single trajectory of the dynamical system under study [8,9]. When observing this trajectory as a scalar time series x(t) (t = 1, . . . , N ), one may use a suitable m-dimensional time delay embedding of x(t) with delay τ [10], 3 x (m) (t) = (x(t), x(t + τ ), . . . , x(t + (m − 1)τ )), for obtaining a recurrence plot as a graphical representation of the binary recurrence matrix where (·) is the Heaviside function, · denotes a suitable norm in the considered phase space and is a threshold distance that should be reasonably small (in particular, much smaller than the attractor diameter [9,11,12]). To simplify our notation, we have used the abbreviation x i = x (m) (t = t i ) (with t i being the point in time associated with the ith observation recorded in the time series 7 ) wherever appropriate. Experimental time series often yield a recurrence plot displaying complex structures, which are particularly visible in sets of recurrence points (i.e. R i, j ( ) = 1) forming diagonal or vertical 'line' structures. A variety of statistical characteristics of the length distributions of these lines (such as maximum, mean, or Shannon entropy) can be used for defining additional quantitative measures that characterize different aspects of dynamic complexity of the studied system in more detail. This conceptual framework is known as recurrence quantification analysis (RQA) [13]- [15] and is nowadays frequently applied to a variety of real-world applications of time series analysis in diverse fields of research [16]. However, most of these RQA measures are sensitive to the choice of embedding parameters, which are found to sometimes induce spurious correlations in a recurrence plot [17].
Recent studies have revealed that the fundamental invariant properties of a dynamical system (i.e. its correlation dimension D 2 and correlation entropy K 2 ) are conserved in the recurrence matrix [18,19], the estimation of which is independent of the particular embedding parameters. What is more, recurrence plots preserve all the topologically relevant phase space information of the system, such that we can completely reconstruct a time series from its recurrence matrix (modulo some rescaling of its probability distribution function) [20,21].
A further appealing paradigm for analysing the structural features of complex systems is based on their representation as complex networks of passive or active (i.e. mutually interacting) subsystems. For this purpose, classical graph theory has been systematically extended by a large variety of different statistical descriptors of the topological features of such networks on local, intermediate and global scales [22]- [24]. These measures have been successfully applied for studying real-world networks in various scientific disciplines, including the structural properties of infrastructures [25], biological [26], ecological [27] and climate networks [28,29], to give some prominent examples. The corresponding results have triggered substantial progress in our understanding of the interplay between the structure and dynamics of such complex networks [30]- [32].
The great success of network theory in various fields of research has recently motivated first attempts to generalize this concept for a direct application to time series [33]- [43]. By means of complex network analysis, important complementary features of dynamical systems (i.e. properties that are not captured by existing methods of time series analysis) can be resolved, which are based on spatial dependences between individual observations instead of temporal correlations. For example, a corresponding idea has been successfully applied for describing causal signatures in seismic activity by means of networks of recurrent events [44]- [46]. However, concerning applications to the field of time series analysis, a substantial number of the recently suggested techniques have certain conceptual limitations, which make them suitable 4 only for dealing with distinct types of problem. As an alternative that may provide a unifying conceptual and practical framework for nonlinear time series analysis using complex networks, we reconsider the concept of recurrences in phase space for defining complex network structures directly based on time series. For this purpose, it is straightforward to interpret the recurrence matrix R( ) as the adjacency matrix A( ) of an unweighted and undirected complex network, which we suggest to call the recurrence network associated with a given time series. To be more specific, the associated adjacency matrix is given by where δ i, j is the Kronecker delta introduced here in order to avoid artificial self-loops, corresponding to a specific choice of the Theiler window in RQA [9]. A related conceptual idea has recently been independently suggested by different authors [35,36,39,40,42], but not yet systematically studied. In this work, however, we aim to give a rigorous and detailed interpretation of a variety of quantitative characteristics of recurrence networks, yielding novel concepts for statistically evaluating distinct phase space structures captured in recurrence plots. In this sense, our presented work continues and complements the developments in the field of RQA within the last two decades [9,43]. It shall be noted that a generalization to weighted networks (as partially studied in [33,34]) is straightforward if the recurrence matrix is replaced by the associated distance matrix between pairs of states. In any case, recurrence networks based on the mutual phase space distances of observational points on a single trajectory are spatial networks, i.e. fully embedded in an m-dimensional space, which has important implications for their specific topological features. We will raise this point in detail within this paper.
The consideration of recurrence plots as graphical representations of complex networks allows a reinterpretation of many network-theoretic measures in terms of characteristic phase space properties of a dynamical system. According to the ergodicity hypothesis, we suppose that we can gain full information about these properties by either ensembles of trajectories, or sufficiently long observations of a single trajectory. Following this line of ideas, we approximate the (usually unknown) invariant density p(x) (which is related to the associated invariant measure µ by dµ = p(x)dx) of the studied system by some empirical estimatep ( ) (x) obtained from a time series, where defines the level of coarse-graining of phase space involved in this procedure. Transforming the time series into a recurrence network then allows us to quantitatively characterize the higher-order statistical properties (i.e. properties that are based on joint features of different states) of the invariant density p(x) by means of network-theoretic measures.
According to the above argumentation, quantitative descriptors of the topological features of recurrence networks can be considered as novel and complementary measures. Specifically, in contrast to the estimation of both dynamical invariants (in particular, K 2 , D 2 and related nonlinear characteristics), as well as traditional RQA measures that are exclusively based on line structures in recurrence plots (i.e. the presence of temporal correlations between different parts of a trajectory), the network-theoretic approach to quantitatively characterizing recurrence plots is distinctively different as it relates to spatial dependences between state vectors in some appropriate phase space. As we will demonstrate in this paper, our technique exhibits additional deep insights into some phase space properties of dynamical systems such as the local attractor fragmentation or the presence of dynamically invariant objects, which are directly related to their complex dynamics, but are not provided by other existing methods. For this purpose, we will take seriously the duality of adjacency matrices of complex networks, on the one hand, and recurrence matrices of dynamical systems, on the other hand, which would in turn also allow transferring concepts from dynamical systems theory (given that the corresponding recurrence plot-based estimates are invariant under temporal reordering) to complex networks. In this work, however, we will concentrate on a detailed discussion of how phase space properties can be further quantified in terms of network theory. The remainder of this paper is organized as follows. Section 2 presents a review of existing approaches for extracting complex networks from time series, including a comprehensive discussion of their potentials and possible conceptual problems (with a special emphasis on how to interpret the resulting networks' topological properties). The concept of recurrence networks as a natural alternative is further discussed in section 3. In particular, it is demonstrated that many network-theoretic measures yield sophisticated quantitative characteristics corresponding to certain phase space properties of a dynamical system that have not yet been explicitly studied in terms of other dynamical invariants or measures of complexity based on RQA. In order to support our theoretical considerations, section 4 provides some examples of how network-theoretic measures reveal phase space properties of various low-dimensional dynamical systems. Finally, we summarize our main results and outline future directions of further methodological developments based on our proposed technique.

Approaches for transforming time series into complex networksa comparative review
In this section, a brief review and classification of existing approaches for studying the properties of time series by means of complex network methods is presented in chronological order (see table 1). In particular, the strengths and possible limitations of the existing approaches will be discussed.

Partitioning of phase space: transition networks
The concept of symbolic dynamics [47] allows characterizing the properties of a dynamical system based on a partition {S 1 , . . . , S K } of its phase space into K mutually disjoint sets, yielding a transformation of every possible trajectory into an (in principle, infinite) sequence of abstract symbols. Formally, the proper application of concepts of symbolic time series analysis (such as mutual information or other entropic quantities) requires the existence of a generating 6 partition that corresponds to a unique assignment of symbolic sequences (i.e. sequences of class identifiers) to every trajectory of the system. Note that this prerequisite is usually violated in real-world applications due to the presence of noise; however, even in the ideal noise-free case, generating partitions either do not exist or can hardly be estimated (see [48] and references therein). Nevertheless, applications of symbolic time series analysis have recently attracted considerable interest in numerous applications [49]- [52].
Partitioning of the phase space of a dynamical system can also be used for transforming a time series into a complex network representation. In the simplest possible case, we identify the different sets S i with the vertices of a network and consider the n-step transition probabilities for characterizing the edges of a fully connected weighted and directed graph [53]. By setting a suitable threshold p min to these transition probabilities, we obtain the adjacency matrix of an associated unweighted, directed network, which we suggest referring to as the n-step transition network, by setting We note that related graph-theoretic methods have been recently used for the set-oriented approximation of almost invariant sets in dynamical systems [54,55]. Recent applications of transition networks in the field of time series analysis include the characterization of transitions in cellular automaton models of vehicular traffic by means of the associated degree distributions [56]- [59] (which has been recently generalized to the consideration of weighted networks [60]), as well as the analysis of stock exchange time series [61,62] in terms of identifying vertices with the highest relevance for information transfer using betweenness centrality and inverse participation ratio. From the time series analysis point of view, the main disadvantage of the transition network approach is that it induces a significant loss of information on small amplitude variations. In particular, two observations with even very similar values are not considered to belong to the same class if they are just separated by a class boundary. This might influence the quantitative features of a corresponding network, as it is not exclusively determined by the widths of the individual classes, but also by their specific definition. Corresponding results are known for other methods of symbolic time series analysis such as estimates of mutual information [63,64] or block entropies [52,65]. In this respect, the recurrence network approach introduced in this work is more objective as it only depends on a single parameter . Note, however, that coarsegraining might be a valid approach in the case of noisy real-world time series, where extraction of dynamically relevant information hidden by the noise can be supported by grouping the data.

Cycle networks for pseudo-periodic time series
In 2006, Zhang and Small [33]- [35] suggested to study the topological features of pseudoperiodic time series (representing, for example, the dynamics of chaotic oscillators like the Rössler system) by means of complex networks. For this purpose, individual cycles (defined by minima or maxima of the studied time series) have been considered as vertices of an undirected network, and the connectivity of pairs of vertices has been established by considering a generalization of the correlation coefficient to cycles of possibly different length or, alternatively, their phase space distance.

7
A potential point of criticism concerning this method is that the definition of a cycle is not necessarily straightforward in complex oscillatory systems. In [33]- [35], the authors mainly considered nonlinear oscillators in their phase-coherent regimes; however, it is not clear how a cycle could be defined for non-phase-coherent oscillations, for example, in the Funnel regime of the Rössler system. The same problem arises for systems with multiple timescales, which are hence hard to treat this way. Furthermore, it is not intuitively clear how to interpret correlations of cycles, since the values of the corresponding measures are not exclusively determined by the proximity of the corresponding parts of the trajectory in phase space, but depend also on the specific choice of sampling. This could yield rather different estimates of the cycle correlation coefficient (or, alternatively, the phase space distance) between two cycles, even if the two parts of the trajectory are very close to each other.

Correlation networks of embedded state vectors
A generalization of the cycle networks method used by Zhang and Small that can also be applied to time series without obvious oscillatory components has been suggested by Yang and Yang [37] using a simple embedding of an arbitrary time series. In their formalism, individual state vectors x (m) i in the m-dimensional phase space of the embedded variables are considered as vertices, from which the Pearson correlation coefficient can be easily computed as withx being the zero-mean components of the respective embedding vectors (for scalar time series, the generalization to vector-valued time series is straightforward). If r i, j is larger than a given threshold, the vertices i and j are considered to be connected, resulting in an undirected network representation. An equivalent approach called fluid-dynamic complex networks (FDCN) has recently been proposed by Gao and Jin [42,66] and successfully applied for characterizing the nonlinear dynamics of conductance fluctuating signals in a gas-liquid two-phase flow. One potential conceptual problem of this particular technique is that the consideration of correlation coefficients between two phase space vectors usually requires a sufficiently large embedding dimension m for a proper estimation with low uncertainty (more specifically, the standard error of the correlation coefficient is approximately proportional to 1/ √ m − 1). Hence, local information about the short-term dynamics captured in a time series might get lost when following this approach. Even more, since embedding is known to induce spurious correlations to a system under study, the results of the correlation method of network construction may suffer from related effects.
With respect to the interpretation of the resulting network patterns, we note that for vertices corresponding to mutually overlapping time series segments, the consideration of correlation coefficients, as applied in the papers cited above, corresponds to studying the local autocorrelation function of the signal. Hence, the presence of edges between these vertices is exclusively determined by linear correlations within the signals. In principle, we might think of replacing the correlation coefficient by other measures of interrelationships such as the mutual 8 information, which are also sensitive to general statistical dependences [28,29,67]; however, the appropriate estimation of such nonlinear quantities would require a considerably larger amount of data, i.e. a very large embedding dimension m.
Finally, when studying time series with pronounced cycles (like trajectories of the Rössler system), different cycles could be included in one embedding vector (depending on the sampling rate), which casts additional doubts with respect to the direct interpretability of the resulting network properties.

Visibility graphs
An alternative to the latter two threshold-based concepts has been suggested by Lacasa et al [38] in terms of the so-called visibility graph. In this formalism, individual observations are considered as vertices, and edges are introduced whenever a partial convexity constraint is fulfilled, i.e. x(t a ) and holds. Since the visibility condition is symmetric with respect to a and b, visibility graphs are undirected. Visibility graphs have been used to study the behaviour of certain fractal as well as multifractal stochastic processes [68,69], energy dissipation in three-dimensional turbulence [70] and the nonlinear properties of financial [71,72] and environmental time series [73]. Recently, a slightly modified approach of horizontal visibility graphs has been proposed and applied to certain random time series [74]. Although visibility graphs are easily established and allow an alternative estimation of the Hurst exponent of fractal time series, an interpretation of the convexity constraint in terms of other phase space properties of the considered system has not yet been provided. Moreover, in its present formulation, the application of this approach is restricted to univariate time series.

Complex networks based on neighbourhood relations in phase space
As it has already been mentioned, the transformation of time series into complex networks by means of neighbourhood relationships has already been discussed by different authors. In particular, there are two possible approaches that can directly be related to slightly different definitions of recurrence plots [9].
On the one hand, a neighbourhood can be defined by a fixed number k of nearest neighbours of a single observation, i.e. a constant 'mass' of the considered environments [35,36]. We refer to this method as a k-nearest neighbour network in phase space. This setting implies that the degrees k v of all vertices in the network are kept fixed at the same value k. Hence, information about the local geometry of the phase space, which is mainly determined by the invariant density p(x), cannot be directly obtained by most traditional complex network measures (see section 3) 8 . Note that the adjacency matrix of a k-nearest neighbour network is in general not symmetric, i.e. the fact that a vertex j is among the k-nearest neighbours of a vertex i does not imply that i is also among the k-nearest neighbours of j. Hence, k-nearest neighbour networks are directed networks.

9
On the other hand, we can define the neighbourhood of a single point in phase space by a fixed phase space distance, i.e. considering a constant 'volume' [39,42,43,75]. This approach has the advantage that the degree centrality k v gives direct information about the local phase space density (see section 3). The consideration of neighbourhood relationships within a fixed phase space volume corresponds to the standard definition of a recurrence plot as mentioned in the introduction section. The resulting networks will therefore be referred to as recurrence networks in the following. Unlike the k-nearest neighbour networks, these recurrence networks are undirected graphs by definition. We note that Gao and Jin [42] termed an equivalent approach as fluid-structure complex networks (FSCN) and used it for analysing gas-liquid two-phase flow and the Lorenz system as a toy model in terms of link density. In addition, they related their observations to the presence of unstable periodic orbits (UPOs) in a dynamical system. We will come back to this point in section 4.3.
We note that since all complex network approaches based on the proximity of different parts of the trajectory (i.e. cycle, correlation, k-nearest neighbour and recurrence networks) do not preserve information about the temporal order of the respective state vectors (i.e. are invariant with respect to random permutations of these vectors), the resulting characteristics of network topology represent dynamically invariant properties related to the state density in the corresponding (abstract) phase space. As a consequence, considering a recurrence network based directly on a univariate time series (or a one-dimensional dynamical system), we cannot distinguish between stochastic and deterministic systems with the same phase space density (for a simple example, see appendix A.2). However, in such cases, embedding or consideration of phase space dimensions larger than one provides a feasible solution for this identification problem. Furthermore, as we will show in the following sections, the statistical characterization of recurrence networks provides useful and complementary insights into the phase space structures that cannot be obtained by other existing approaches of nonlinear time series analysis.

Recurrence networks as a unifying framework for complex network-based time series analysis
Following the considerations in the previous subsections, many existing methods for transforming time series into complex networks suffer from specific conceptual limitations. In particular, the concepts used for defining both vertices and edges of the networks, which differ across the various techniques, are in some cases rather artificial from a dynamical systems point of view (table 1).
In contrast to the other recently suggested approaches, the identification of a recurrence matrix with the adjacency matrix of a complex network is a straightforward and natural idea that conserves many local properties of the invariant density of the studied dynamical system captured in single observational time series. In particular, individual values of the respective observable can be directly considered as vertices of the recurrence network (similar to the visibility graph concept), while the existence of an edge serves as an indicator of a recurrence, i.e. pairs of states whose values do not differ by more than a small value in terms of a suitable norm in phase space.
It should be pointed out that not only recurrence networks but also other existing methods of complex network-based time series analysis are based on the concept of recurrences in phase space. Apart from the k-nearest neighbour networks originating from the idea of a fixed local recurrence rate (i.e. a fixed mass of the neighbourhoods), the idea of considering a threshold Table 2. Relationships between recurrence network entities and corresponding geometrical objects and their properties in phase space.

Recurrence network
Phase space Recurrence of states Path Overlapping sequence of -balls value to the proximity of two vertices can also be found in other previously suggested methods.
In particular, correlation networks [37] (see section 2.3) can be regarded as a special case of recurrence networks where the usual metric distance has been replaced by the correlation distance [76] where r i, j is the correlation coefficient between two embedded state vectors defined in equation (5). Similar considerations apply to the case of cycle networks. Note, however, that the advantage of using the concept of recurrences of individual states in phase space defined by metric distances instead of correlations or cycle properties associated with state vectors of dynamical systems is that it allows for creating networks without any embedding or consideration of groups of states. On the one hand, this independence from a particular embedding is beneficial when dynamical invariants of the studied system are of interest. On the other hand, the statistical properties of the resulting recurrence networks reflect exclusively the invariant density of states in phase space (in terms of certain higher-order statistics), because time-ordering information is lost in this framework. Following the above considerations, it can be argued that the concept of recurrence networks yields a unifying framework for transferring time series into complex networks in a dynamically meaningful way. In particular, this approach can be applied (i) to both univariate and multivariate time series (phase space trajectories), (ii) with and without pronounced oscillatory components and (iii) with as well as without embedding 9 . Moreover, similar to traditional RQA, studying network properties for sliding windows in time also allows for coping with non-stationary time series [43]. Consequently, unlike for most of the existing techniques, there are no fundamental restrictions with respect to its practical applicability to arbitrary time series.

Quantitative characteristics of recurrence networks
While the definitions of edges and vertices in our approach have already been given above (table 1), we now provide a geometrical interpretation of a third important network entity, the path, within the framework of recurrence networks (table 2). A path between two vertices i to j in a simple graph without multiple edges can be written as an ordered sequence of the vertices it contains, i.e. (i, k 1 , . . . , k l i, j −1 , j), where the associated number of edges l i, j measures the length of the path. In phase space, a path in the recurrence network is hence defined  10 Due to the natural interpretation of vertices, edges and paths, the topological characteristics of a recurrence network closely capture the fundamental phase space properties of the dynamical system that has generated the considered time series. In the following, we will present a detailed analysis of the corresponding analogies for different network properties that are defined on a local (i.e. considering only the direct neighbourhood of a vertex), intermediate (i.e. considering the individual neighbourhoods of the neighbours of a vertex) and global (i.e. considering all vertices) scale (table 3) 11 . It has to be emphasized that these quantities can be considered as (partly novel and complementary) measures within the framework of RQA.

Local network properties
3.1.1. Degree centrality (local recurrence rate). As a first measure that allows us to quantify the importance of a vertex in a complex network, the degree centrality [77] of a vertex v, k v , is defined as the number of neighbours, i.e. the number of vertices i = v that are directly connected with v: Note that, in general, the sum is taken over all i = v. However, according to our definition (2), we skip the corresponding condition in the following. Normalizing this measure by the 12 maximum number of possible connections, N − 1, we obtain the local connectivity which, from the recurrence plot point of view, corresponds to the local recurrence rate R R v of the state v. Thus, degree centrality and local connectivity yield an estimator for the local phase space density, since for a vertex v located at position (when using the maximum norm) and, hence, In complex network studies, one is often interested in the frequency distribution of degree centralities, P(k), in particular, the presence of an algebraic scaling behaviour, which is a characteristic of scale-free networks [22]. In the case of recurrence networks, however, we note that the emergence of a scale-free property in the degree distribution requires distinct structures of the phase space density, which are probably only present in rather specific cases. Moreover, although several authors have recently focused their attention on this characteristic obtained from different types of complex networks derived from time series [33], [37]- [39], [42], [56]- [59], [68]- [70], we would like to underline that for a complete characterization of the phase space properties of a dynamical system, not only degree centralities but also other higher-order statistical measures have to be studied.

Edge density (global recurrence rate).
In some situations, it is useful not to consider the full distribution of degree centralities in a network, but to focus on the mean degree of all vertices as a simple characteristic quantity of this distribution, where is the total number of edges in the recurrence network. The mean degree centrality k is directly proportional to the edge density ρ of the network or, alternatively, its recurrence plot equivalent, the global recurrence rate R R, Note that the recurrence rate coincides with the definition of the correlation integral C 2 ( ), which is commonly used to estimate the correlation dimension D 2 , for example, using the Grassberger-Procaccia algorithm [78]. The connection between the edge density and the correlation dimension can be understood by the fact that the local recurrence rate R R v of a vertex v corresponds to the measure of an m-dimensional ball B (x v ) of radius centred at the point x v in the m-dimensional phase space in the limit that time goes to infinity (N → ∞). When considering the Euclidean norm as a distance measure in phase space, these balls are defined as hyperspheres, for the maximum norm as hypercubes, etc. Then, the pointwise (information) dimension of the probability measure µ at [4]. Owing to the heterogeneity of the phase space visited by the trajectory (i.e. the non-uniform phase space density that results in different degree centralities k v in different parts of this space), the proper estimation of D p is a nontrivial task and often requires expensive computational power and a high data quality and quantity. Thus, we expect a better statistics for D 2 , since it more heavily weights regions of the phase space that have a higher probability measure µ. Although the correlation integral has been well established in the literature for estimating the correlation dimension, we point out the improvement in estimating D 2 based on the diagonal lines in R i, j ( ), which yields an algorithm that is independent of the embedding parameters [19]. Consequently, the recurrence network representation A i, j of a time series fully conserves the geometric properties of the phase space of the underlying dynamical system.

Intermediate scale network properties
3.2.1. Local clustering coefficient. The clustering coefficient of a vertex v, C v , characterizes the density of connections in the direct neighbourhood of this vertex in terms of the density of connections between all vertices that are incident with v. In many networks, such loop structures formed by three vertices occur more often than one would expect for a completely random network. Hence, high clustering coefficients reveal a specific type of structure in a network, which is related to the cliquishness of a vertex [24].
In this work, we consider the definition of clustering coefficient proposed by Watts and Strogatz [24], where N v is the total number of closed triangles including vertex v, which is bound by the maximum possible value of k v (k v − 1)/2. For vertices of degree k v = 0 or 1 (isolated or treelike points, respectively), the clustering coefficient is defined as C v = 0, as such vertices cannot participate in triangles by definition. Equation (16) can be rewritten in terms of conditional probabilities as using Bayes' theorem with 14 and a similar expression for P(A i, j = 1, A v,i = 1, A v, j = 1). As for a recurrence network, the value of A i, j depends only on the phase space distance and the choice of ; the latter relationship can be used to derive analytical results at least for one-dimensional systems based on their invariant density. Corresponding details can be found in appendix A, including the corresponding treatment of the Bernoulli and logistic maps as specific examples.

Global clustering coefficient.
As for the degree centrality, we consider the average value of the clustering coefficient taken over all vertices of a network, the so-called global clustering coefficient as a global characteristic parameter of the topology of a network. We expect that the value of C is-for a given dynamical system with a phase space density p(x)-in the asymptotic limit N → ∞ exclusively determined by the choice of , which defines the scale of resolution. A more detailed discussion of the corresponding effects and their implications for certain model systems will be given in section 4.

Mean nearest neighbour degree.
The mean nearest neighbour degree k nn v of vertex v gives the average degree in the neighbourhood of v, The degree centrality k v is a measure of the density of states in the immediate neighbourhood of state v, whereas k nn v can be interpreted to indicate the mean density of states in the next neighbourhood (next topological shell of neighbours) of state v. Hence, both measures taken together contain information about the local density anomaly in the vicinity of v, which we propose to measure by the local degree anomaly Vertices with a positive degree anomaly ( k v > 0) hence indicate local maxima of phase space density, while those with k v < 0 correspond to local density minima. Hence, the local degree anomaly may be considered as a proxy for the local heterogeneity of the phase space density. In a similar way, the average absolute value of the local degree anomaly, | k v | v , serves as a measure of the overall spatial heterogeneity of the phase space density profile.

Assortativity.
A network is called assortative if vertices tend to connect preferentially to vertices of a similar degree k. On the other hand, it is called disassortative if vertices of high degree prefer to connect to vertices of low degree, and vice versa. Hence, assortativity can be quantified by the Pearson correlation coefficient of the vertex degrees on both ends of all edges [24,79], If the density of states in phase space hardly varies within an -ball, the degrees on either ends of an edge will tend to be similar and hence the assortativity coefficient R will be positive. This means that the more continuous and slowly changing the density of states is, the closer R will be to its maximum value one. Within the framework of recurrence networks, R can hence be interpreted as a measure of the continuity of the density of states or put differently of the fragmentation of the attractor. Note that this aspect has not yet been specifically addressed by other nonlinear measures, in particular, within the RQA framework.

Matching index (twinness).
The overlap of the neighbourhood spaces of two vertices i, j is measured by the matching index where µ i, j = 0 if there are no common neighbours, and µ i, j = 1 if the neighbourhoods coincide [24]. Using the notion of -balls around points in phase space, we can alternatively write .
Owing to the spatial constraints of the recurrence network, the neighbourhood spaces of i, j can only overlap if i.e. µ i, j = 0 for all pairs of vertices (i, j) with d i, j > 2 . Moreover, the matching index µ i, j decreases on average with increasing spatial distance d i, j between the two considered states. Note that since A i, j = 0 already for d i, j > , unconnected points with a matching index µ i, j > 0 may exist. The matching index of pairs of vertices in a recurrence network is closely related to the concept of twins [80], which has recently been successfully applied for constructing surrogate data (twin surrogates) in the context of statistical hypothesis testing for the presence of complex synchronization [81,82]. Twins are defined as two states of a complex system that share the same neighbourhood in phase space, i.e. the two vertices of the recurrence network representing these states have a matching index µ i, j = 1. Hence, the matching index can be used for identifying candidates for twins. Note that pairs of vertices i and j in a recurrence network with µ i, j 1 can still be considered as potential twins, since µ i, j = 1 is in some cases approached by only slight changes of the threshold . Consequently, we suggest interpreting the matching index as a measure of the twinness of i and j. Furthermore, it should be noted that adjacent pairs of edges (i, j) (A i, j = 1) with a low matching index µ i, j 0 connect two distinct regions of the attractor and may therefore be indicative of geometrical bottlenecks in the dynamics (cf our discussion of the betweenness centrality in section 3.3.5).

Global network properties
3.3.1. Shortest path length. As we consider recurrence networks as undirected and unweighted, we assume all edges to be of unit length in terms of graph (geodesic) distance. Consequently, the distance between any two vertices of the network is defined as the length of the shortest path between them. Note that time information is lost after transforming the trajectory into a network presentation. Therefore, the terminology of the shortest path length l i, j in the recurrence network reflects the minimum number of edges that have to be passed on a graph from a vertex i to a vertex j. In the same spirit, l i, j is related to the distance of states i and j in phase space.
In order to better understand the meaning of shortest path lengths, let us study their calculation for two toy model series: first, we consider a periodic trajectory x(t) = sin(6π · 0.1t), y(t) = cos(6π · 0.1t), with t = 0, 1, . . . , 10, i.e. there are N = 10 points in the phase space ( figure 1(a)). The corresponding recurrence plot is shown for = π 5 ( figure 1(b)). As it has already been mentioned above, the recurrence matrix R i, j and the adjacency matrix A i, j of the associated recurrence network are basically equivalent. Adopting a common visualization of connectivity patterns from the literature on complex networks, we illustrate the recurrences  1(c)). In this representation, the shortest path length (in the network sense) between two vertices i and j corresponds to the smallest number of 'jumps' in phase space via pairs of neighbours (i.e. recurrences) in phase space. Obviously, the number of such jumps is determined by the prescribed value of and the spatial distance between i and j. For instance, the shortest path from vertex i = 1 to j = 10 is l 1,10 = 3 as indicated by the matrix of mutual shortest path lengths ( figure 1(d)). Note that this list is symmetric by definition, i.e. l i, j = l j,i . The same heuristic analysis can also be performed for a general non-periodic trajectory in phase space as shown in figures 2(a)-(d).
We wish to underline that the terminology of shortest path lengths in networks does not have a direct relevance for the dynamical evolution of the observed system. In contrast, l i, j measures distances among a discrete set of points on the attractor in units of the neighbourhood size . For example, in the periodic case displayed in figure 1, it takes nine iterations (time points) from vertex 1 to vertex 10 in the time domain, while the shortest path to cover the phase space distance has only a length of l 1,10 = 3. Hence, shortest paths do not allow to infer the temporal evolution of the system. Even more, for the path concept in a recurrence network, no information about the temporal order of the individual observations is considered (for example, the shortest path between vertices 1 and 2 in figure 1 is given by the sequence (1, 8, 5, 2), which is not ordered in time).
One should note that if the phase space is strongly fragmented (for instance, in the period-3 window of the logistic map, which has been discussed elsewhere [43], the phase space consists of three discrete points), the resulting recurrence networks can be composed of different disconnected clusters. Furthermore, there might be more than one shortest path connecting two nodes. For example, in the aperiodic example in figure 2, the shortest path from node 1 to node 7 (l 1,7 = 3 as shown in figure 2(d)) can be obtained by three different choices: (1, 2, 6, 7), (1, 5, 3, 7) and (1, 5, 6, 7).

Average path length.
The average path length L is defined as the mean value of the shortest path lengths l i, j taken over all pairs of vertices (i, j), Here, for disconnected pairs of vertices, the shortest path length is set to zero by definition. Note that in most practical applications, this has no major impact on the corresponding statistics. The average phase space separation of states d i, j i, j serves as an -lower bound to L, since d i, j l i, j (27) due to the triangular inequality, and hence Interpreted geometrically, this inequality holds because L approximates the average distance of states along geodesics on the recurrence network graph (which can be considered as the geometric backbone of the attractor) in multiples of , while d i, j i, j gives the mean distance of states in R m as measured by the norm · . (27) for the average path length, the diameter

Network diameter. By a similar argument as used in equation
of the recurrence network (i.e. the maximum path length) serves as an -upper bound to the estimated diameter of the attractor in phase space:

Closeness centrality.
The inverse average shortest path length of vertex v to all others in the recurrence network is measured by the closeness centrality [77] c If there is no path connecting two vertices i and j, the maximum shortest path length in the graph, N − 1, is used in the sum by definition. In a recurrence network, c v can be geometrically interpreted as measuring the closeness of v to all other states with respect to the average length (in units of ) of geodesic connections on the recurrence network graph. In other words, c v is large if most of the other vertices are reachable in a small number of -jumps from state to state. From equations (27) and (32), we can see that the inverse closeness c −1 v is bounded from below by the average phase space distance of vertex (state) v to all other vertices (states) d v,i i in units of (geometrical closeness), as measured by the norm · , Put differently, geometrical closeness provides an upper bound for topological closeness,

Betweenness centrality.
Let σ i, j be the number of shortest paths between two vertices i and j, and σ i, j (v) the number of such paths that include a specific vertex v of the network. Then, the betweenness centrality of v is defined as [83] b In addition to degree and closeness centralities, betweenness centrality yields another possibility to identify especially relevant vertices. Note that unlike the degree centrality k v , it is defined locally but depends on global adjacency information. For general complex networks, the meaning of b v can be understood in terms of the importance of individual vertices for the transport of information or matter, assuming that both typically travel through the network on shortest paths. In general, there are σ i, j different shortest paths connecting two vertices i and j. We then regard a vertex v to be an important mediator for the transport on the network, if it is traversed by a large number of all existing shortest paths. According to these conceptual ideas, in equation (35), the contribution of shortest paths is weighted by their respective multiplicity σ i, j , the physical rationale for this normalization being that the total volume of information flow between two vertices, when summed over all shortest paths connecting them, should be the same for all pairs in the network.
For a recurrence network, the notion of information transfer is not useful anymore. However, we can still argue in a geometric way that high betweenness states are typical for regions of sparse phase space density that separate different high-density clusters (referring to the information flow analogy mentioned above, we consider the corresponding vertices as geometric bottlenecks). Thus, the occurrence of high betweenness values can be a sign of highly fractionated attractors (on the scale resolved by the considered threshold ). A more detailed discussion of the corresponding implications for some simple model systems will be given in section 4.

Edge betweenness.
While betweenness centrality refers to vertex properties of a network, an equivalent measure can also be defined based on the number of shortest paths on the network that include a specific edge (i, j). We refer to the corresponding property as the edge betweenness b i, j . Note that although there is a conceptual difference between vertex-related and edge-related betweenness, both quantities are indicators for regions of low phase space density that separate regions with higher density (or, to phrase it differently, of regions of high attractor fractionation) and thus have practically the same dynamical meaning.

Examples
In the following, we will show the potentials of the network-theoretic measures discussed in the previous section for recurrence networks obtained from three paradigmatic chaotic model systems.

Model systems
Basic results for one-dimensional maps have already been described for the logistic map (see [43]) based on numerical calculations and are supplemented by some further computations in the appendix. At this point, we prefer to discuss in some more detail the properties of systems that are defined in somewhat higher dimensions. In particular, we consider the Hénon map as an example for a chaotic two-dimensional map, and the Rössler system as well as the Lorenz system as two examples for three-dimensional chaotic oscillators. In all the following considerations, no additional embedding will be used. Note, however, that for the continuous systems, temporal correlations between subsequent observations have been excluded by removing all sojourn points from the recurrence matrix [9,84]. Figures 3 and 4 show examples of typical trajectories of these three model systems. In addition, the shortest paths between the first and last points of the individual realization are indicated, underlining the deep conceptual differences between the concepts of the trajectory (in phase space) and path (in a recurrence network, see section 3).

-dependence of network properties
In order to evaluate the robustness of the topological properties of recurrence networks, their dependence on the free parameter of the method, , has to be explicitly considered [12]. In the following, we will briefly summarize the main findings for the three model systems.  figure 5 and verify the existence of an inverse relationship of a corresponding lower bound postulated in equation (28).
For the global clustering coefficient, the dependence on is more complicated and depends on the specific properties of the considered system (figure 5). In particular, while for too small , problems may occur, since the recurrence network will generally decompose into different disconnected clusters for a length N of the considered time series, for intermediate threshold values, an approximately linear increase of C with seems to be a common feature of all three examples. Following the discussion of the behaviour of one-dimensional maps in appendix A, we argue that this increase is most likely related to the effect of the attractor boundaries.
Finally, concerning the assortativity coefficient R, we observe that for small , the recurrence networks are highly assortative (e.g. R is close to 1). This behaviour can be related to the fact that in the case of small neighbourhoods, these phase space regions are usually characterized by only weak variations of the phase space density, so that neighbouring vertices have a tendency to obey a similar degree. As becomes larger, larger regions of the phase space are covered, where the density may vary much stronger, which implies that the degrees of neighbouring vertices become less similar. Note, however, that since in this case, the mutual overlap of the different neighbourhoods becomes successively larger, there is still a significantly positive correlation between the degrees of neighbouring vertices. We further observe that the decrease of R with is interrupted by intermediate increases, which are probably related to some preferred spatial scale of the separation of certain dynamically invariant objects such as UPOs. We will come back to this point in section 4.3.3. (Note that the square gives an idealized representation of the considered neighbourhood.) The red line displays the shortest path between the initial condition and the final value on this trajectory (l = 5). (c, d) One example trajectory of the Lorenz system with T = 5 time steps, and the resulting shortest path between the initial and the final state ( = 1.5, l = 9).

Vertex properties.
Similar to the global network properties discussed above, the properties of the individual vertices of a recurrence network show a corresponding dependence on the recurrence threshold as well. This observation is rather trivial for the degree centrality, where not only the global mean value (proportional to the recurrence rate R R) but also the local value at an individual vertex must be a monotonously increasing function of [9]. We note that for comparing the properties of recurrence networks between different dynamical systems (or for different values of the control parameters of the same system) that generally have attractors with different phase space diameters, it is indeed favourable to fix the recurrence rate R R instead of the metric distance . For a detailed discussion on this question see, e.g., [85].
Concerning other vertex properties such as b v or C v , we note that the spatial distribution of the individual values becomes successively smoother as increases, since the small-scale features of the phase space density captured by these measures cannot be resolved with large recurrence thresholds anymore. However, for sufficiently small , we find that the main  Figure 5. Dependence of (a-c) the average path length L, (d-f), the global clustering coefficient C and (g-i) the assortativity coefficient R for the Hénon map, Rössler and Lorenz system (from left to right) using the maximum norm. Dashed lines in the plots of L( ) indicate the approximate presence of the theoretically expected 1/ dependence of the average path length. Note that although a normalized threshold might yield a better comparability of the results for the different systems, we prefer using the absolute values here since the typical normalization factors-either an empirical estimate of the standard deviation of the phase space density (which may be crucially influenced by asymmetric densities) or the attractor diameter (which is itself not a priori known in advance) have certain conceptual problems.
structures resolved by the local vertex properties are preserved [12]. In the following, we will operate with different values of R R for different systems and network-theoretic measures, depending on the computational requirements and the spatial scales of interest.

Spatial distributions of vertex properties
In the following, we will study the interrelationships between local network properties and structural features of the phase space for the three considered chaotic model systems.

Degree centrality.
When considering the degree centrality k v or, equivalently, the local density ρ v for all vertices of the network, a broad range of variability is found ( figure 6). In particular, the behaviour follows the expectation that regions with a high phase space density (for example, the merger of the two scrolls of the Lorenz oscillator) also reveal a high density of vertices and, hence, high degree centralities. Note that the calculation of a recurrence plot depends on the parameter , which should be tailored to the considered system under study and the specific questions one wishes to address. Several 'rules of thumb' for the choice of the threshold have been advocated in the literature [9,11]. It has been suggested that the choice of to achieve a fixed recurrence rate R R is helpful for the estimation of dynamical invariants in many systems [9]. Therefore, this procedure will be adopted here to obtain an overall visualization of the degree centrality k v in phase space, with R R = ρ ≈ 0.03 (which lies within the typical scaling region of the correlation integral). However, as we will see later, for the local clustering coefficients (section 4.3.3) disclosing local fine structures of the phase space density, it is necessary to choose smaller . Figure 7 reflects the spatial distribution of the closeness centrality c v . In good agreement with our previous theoretical considerations on the geometric meaning of this measure (section 3.3.4), we find high values of c v near the centre of gravity of the attractor in phase space, and low values at phase space regions that have large distances from this centre.

Clustering coefficient.
There is no distinct relationship between the local clustering coefficient C v and the degree centrality k v (figure 8), which implies that the spatial patterns of C v are not primarily related to variations of the phase space density itself. Instead, we argue that C v characterizes higher-order properties related to the heterogeneity of the phase space density in the vicinity of a vertex v. For example, in a two-dimensional system, an alignment of vertices along a one-dimensional subspace will produce clearly higher clustering coefficients than a homogeneous filling of the neighbourhood, since a larger number of neighbours of a specific vertex are also mutually connected in this case (see figure 9). The corresponding behaviour is underlined in figure 10 for the Hénon map, where maximum values of C v = 1 can be particularly found at the two tips of the attractor. Generalizing these considerations to arbitrary dynamical systems, we note that C v for a recurrence network quantifies how randomly the vertices are distributed in a specific part of the phase space. In this sense, C v can be considered as an entropy-related quantity. As a more detailed interpretation of this finding, we emphasize that from the theory of spatial random graphs [88] (which can be assumed to yield the lowest possible clustering coefficients among all networks that are embedded in a space with a prescribed dimension m), it is known that in the asymptotic limit N → ∞ and → 0, the possible values of C v are bound between a theoretical lower limit and 1. Moreover, this lower bound systematically decreases with increasing m, which appears reasonable if for a spatial network, we interpret C v as an entropy-related quantity. To be more specific, according to Dall and Christensen [88], this decay follows an exponential function for sufficiently large spatial dimensions.
The presence of distinct structures in the spatial profile of the local clustering coefficient is related to the emergence of specific dynamically invariant objects in the considered model systems. In the case of the Hénon map (figure 10), there is a clear tendency that points that  The finite length segment of the stable manifold is calculated by the method described in [86,87] (with 20 000 iterations).
are close to the stable manifold associated with the hyperbolic fixed point of the system have remarkably higher values of C v . Note, however, that because of finite size effects, this coincidence cannot be found for all corresponding regions of the phase space. For the two continuous systems (figure 11), points close to the trapping regimes of UPOs have higher clustering coefficients. It is, in some sense, trivial to understand the role of UPOs in forming such regimes of higher clustering coefficients. Whenever a trajectory visits the vicinity of a UPO, it is captured in this neighbourhood for a certain finite time, during which the probability of recurrences is increased. Furthermore, once the trajectory is trapped, the local divergence rate becomes smaller (see figure 9). This smaller local divergence rate is captured by the clustering coefficient in terms of higher-order correlations between neighbours of a vertex. Similar to the finite-effect in the Hénon map, the regions with increased clustering coefficients in most cases only coincide with UPOs of lower periods. Therefore, in figure 11, only a few UPOs of low order are shown for comparison. Note that if two UPOs are separated by a distance smaller than in phase space, the clustering coefficient is not able to distinguish between these two structures and, hence, shows a broad band with increased values. In the limit of → 0 (and, hence, N → ∞), also UPOs of higher periods can be detected by C v . Following this argument, C v is a useful measure for detecting phase space regions with a high density of low-order UPOs, which supports corresponding considerations in [42], where the influence of a finite value of has, however, not been properly considered [12].

Betweenness centrality.
Our interpretation of the betweenness centrality in section 3.3.5 implies that b v is a rather sensitive measure of the local fragmentation of the attractor and thus gives complementary information, especially on very small scales. Unfortunately, numerical limitations on the calculation of this measure did not allow us to explore the limit of small neighbourhoods ( → 0). However, from our computations with somewhat larger thresholds (see figure 12), we can already derive some general statements about the behaviour of betweenness centrality for the considered model systems. Firstly, note that regions close to the outer boundaries of the attractor (in contrast to those in the vicinity of the inner boundaries, e.g. of the Rössler oscillator) are not important for many shortest path connections on the recurrence network. Hence, vertices settled in the corresponding parts of the phase space are characterized by low betweenness values. Secondly, if there are pronounced regions with rather few isolated points in between high-density regions (for example, between two UPOs in the Rössler or Lorenz systems), there is an increasing number of shortest paths crossing these vertices, which leads to higher values of b v . In turn, vertices in the vicinity of UPOs (i.e. high-density regions) typically show lower betweenness values. Therefore, betweenness centrality provides a complementary view on the attractor geometry in comparison with the local clustering coefficient C v ( figure 11). However, we have to emphasize that b v alone does not provide a feasible measure for detecting UPOs in dynamical systems [12].

Spatial distributions of edge properties
Similarly to the local vertex properties, we also illustrate the characteristics of different edges in the recurrence networks. Since the resulting structures are more pronounced than for the three model systems considered so far, figure 13 shows the matching index and edge betweenness for one realization of the logistic map x i+1 = ax i (1 − x i ) in the intermittent chaotic regime (see [43]). The presence of intermittent dynamics can be clearly seen from the recurrence plots in terms of extended square recurrence patterns, which hence lead to mutually connected vertices of the associated recurrence networks that correspond to subsequent points in time. Figures 13 and 14 show the complex dependence between phase space distance d i, j , matching index µ i, j and edge betweenness b i, j . For the matching index, the results are consistent with our theoretical considerations presented in section 3.2.5. In particular, for d i, j → 0, we have µ i, j → 1, while for d i, j → 2 , µ i, j → 0. Concerning the temporal evolution during the laminar (intermittent) phase, we recognize that at the beginning, there is hardly any change in the state of the system; hence, d i, j is very small for subsequent points in time (vertices of the recurrence network), which relates to large values of the matching index near 1. As the laminar phases are close to their termination, chaotic variations emerge and rise in amplitude, which leads to a subsequent increase of d i, j and, hence, decrease of µ i, j .
Concerning the edge betweenness b i, j (the spatial pattern of which is very similar to that of the vertex-based betweenness centrality b v due to the spatial proximity of edge and corresponding vertices), the overall behaviour is opposite to that of µ i, j . During laminar phases, we find that since all states are very close to each other, possible shortest connections may alternatively pass through a variety of different edges, leading to low values of the edge betweenness. Close to the termination, there is in turn an increase of this measure. However, the most interesting feature of the edge betweenness is presented by isolated edges with very high values of b i, j , which correspond to rarely visited phase space regions between intervals of higher phase space density. More specifically, the average edge betweenness of vertices in such low-density regions can exceed that of high-density regions by orders of magnitude ( figure 13).

Conclusions
This paper has reconsidered the analysis of time series from complex systems by means of complex network theory. We have argued that most existing approaches for such an analysis suffer from certain methodological limitations or a lack of generality in their applicability. Scatter diagram of the matching index µ i, j against the phase space distance d i, j for the logistic map at a = 3.679 (parameters as in figure 13).
As an appealing solution, we have suggested recurrence networks as a unifying framework for studying time series as complex networks, which is based on a novel approach for the quantitative assessment of recurrence plots in terms of complex network measures. As we have argued, this specific approach is applicable to univariate as well as multivariate time series with or without embedding. In addition, recurrence networks can be applied for studying time series with non-equidistant timescales as well as temporal variations in non-stationary data (by using sliding windows in time) and allow the construction of simple significance tests with respect to the associated network-theoretic measures [43].
Our main achievement is to have provided a thorough reinterpretation of a variety of statistical measures from network theory computed for recurrence networks in terms of phase space properties of dynamical systems. Since all time ordering information is lost in this approach, all complex network characteristics are dynamically invariant, i.e. they are only sensitive to certain properties of the invariant density of the considered dynamical system. From this invariance, it follows that specific measures such as the local clustering coefficient might be used for detecting dynamically invariant objects like UPOs or chaotic saddles. However, this feature also implies that the proposed method cannot be used to distinguish between deterministic (chaotic) and stochastic systems, which is exemplified by our comparison between the Bernoulli map and uniform noise in appendix A. As a possibile means to overcome this potential point of criticism, we emphasize that an additional embedding should change the properties of deterministic systems in a different way than for a stochastic system; hence, studying complex network properties dependent on the embedding dimension might help to solve this interpretation problem. We will further elaborate this idea in future research.
Using widespread statistical characteristics of complex networks such as the 'trinity' of centrality measures (degree, closeness and betweenness) and the clustering coefficient, we were able to provide a detailed interpretation of the corresponding results for recurrence networks in terms of higher-order properties of the phase space density of a dynamical system. In particular, degree centrality relates to the local density, closeness centrality to the average geometrical proximity of an observation to all other observations, and betweenness centrality to the local fragmentation of points in phase space. For the local clustering coefficient C v , the analysis of different model systems has revealed that the resulting values are related to the homogeneity of the spatial filling, which becomes important close to the attractor boundaries and certain dynamically invariant objects such as invariant manifolds or UPOs. The resulting local network properties are qualitatively robust under changes of the recurrence threshold given that the associated recurrence rate R R remains sufficiently small [12]. Additionally, we have presented a rigorous analytical treatment of the local and global clustering coefficients of the recurrence networks of one-dimensional chaotic maps, which perfectly matches our numerical results (see appendix A). A possible further relationship between the local clustering coefficient of a recurrence network and other nonlinear characteristics of the underlying dynamical system such as the local Lyapunov exponent or the point-wise dimension remains a topic for future work.
With respect to existing recurrence plot-based methods of time series analysis, e.g. RQA, we would like to emphasize that our approach yields a complementary view on the phase space properties of the underlying dynamical system. In particular, we note that nearly all of the considered network-theoretical measures have no direct equivalent in traditional RQA, and vice versa. Hence, there may be situations when either one of the two frameworks (i.e. RQA or recurrence networks) provides better results than the other. Examples include the detection of dynamical transitions between periodic and chaotic behaviour in maps [43], continuous dynamical systems [85] or the analysis of protein folding [89,90]. In [85], an additional comparison with the standard approach based on numerical estimates of the maximum Lyapunov exponents has also been presented. In turn, there have been some recent approaches to applying methods of time series analysis to general complex networks (e.g. [91,92]). We would like to underline that due to the duality between the recurrence matrix of a time series and the adjacency matrix of the associated recurrence network, RQA might be another promising candidate for this purpose (as long as we restrict ourselves to measures that are invariant under random permutations of observations), yielding interesting complementary insights on complex networks in a variety of different situations. A more detailed investigation of the potential of a corresponding approach will be the subject of our future research.
Finally, we would like to point out that although our studies in this work have only been supported by results for low-dimensional systems, the proposed methods can also be applied to higher-dimensional systems, e.g. for detecting dynamically relevant phenomena like coherent structures or laminar phases. However, in such cases, we expect some technical problems related to computational demands (associated with certain basic requirements to the number of observations) for calculating different network properties (especially when additional embedding is necessary), the removal of sojourn points, and the visualization of the resulting vertex properties, which could present serious practical challenges for applications of recurrence networks. Nevertheless, an application of recurrence networks to carefully chosen low-dimensional subsets of observables from a high-dimensional system can still yield insightful results, as recently exemplified by a study of a palaeoclimate time series [43]. we obtain the following expressions for the local clustering coefficient C v = C(x v ) if 0.5: 0 x v : x v 1 − : 1 − x v 1: Understanding the global clustering coefficient as the expectation value of the local one (taken over the whole possible range in x), we use the following expression for deriving its value:

A.2. The Bernoulli map
Among all nonlinear maps defined on the unit interval [0, 1], the Bernoulli map x i+1 = 2x i mod 1 has the simplest possible invariant density p(x) ≡ 1. 13 This allows an easy evaluation of the integrals for an analytic computation of the local clustering coefficient, being aware of the restricted integration range. As a result, we find: In particular, for → 0 (and N → ∞), we have C v → 3 4 ∀x ∈ [0, 1], which corresponds to the value of random geometric graphs in one dimension [88]. Hence, one may speculate about this value being a universal limit for the recurrence networks of one-dimensional chaotic maps. Moreover, for x → 0 and x → 1, we have C v → 1 independent of . This behaviour is consistent with our previous observations concerning the effect of sharp attractor boundaries on the local clustering coefficient, for example, in the case of the tips of the Hénon attractor ( which is in excellent agreement with numerical results (see figure A.2). Hence, we argue that deviations from the theoretical value 3 4 occur exclusively due to boundary effects, and can even lead to C = 1 for very large thresholds .

A.3. Logistic map for a = 4
The logistic map at a = 4 is known to have the invariant density .  1]. Note that due to the same distribution of the data, i.e. the same density in the one-dimensional phase space, both empirical curves (circles) are equal and match the theoretical expectations (dashed red lines). However, note that already a two-dimensional embedding would lead to significant differences between both systems (not shown here).
With this, we find the following expressions: and 2 )) , (A.14) systematic tendency towards larger values close to the boundaries. Since the measure of the latter intervals systematically increases with increasing , these boundary effects are again responsible for the systematic increase in the global clustering coefficient C as becomes larger. In turn, the point-wise convergence of C v → 3 4 for → 0 underlines the supposed universality of this value for chaotic one-dimensional maps in agreement with corresponding theoretical predictions for spatial random graphs [88].