Many complex systems exhibit periodic oscillations comprising slow–fast timescales. In such slow–fast systems, the slow and fast timescales compete to determine the dynamics. In this study, we perform a recurrence analysis on simulated signals from paradigmatic model systems as well as signals obtained from experiments, each of which exhibit slow–fast oscillations. We find that slow–fast systems exhibit characteristic patterns along the diagonal lines in the corresponding recurrence plot (RP). We discern that the hairpin trajectories in the phase space lead to the formation of line segments perpendicular to the diagonal line in the RP for a periodic signal. Next, we compute the recurrence networks (RNs) of these slow–fast systems and uncover that they contain additional features such as clustering and protrusions on top of the closed-ring structure. We show that slow–fast systems and single timescale systems can be distinguished by computing the distance between consecutive state points on the phase space trajectory and the degree of the nodes in the RNs. Such a recurrence analysis substantially strengthens our understanding of slow–fast systems, which do not have any accepted functional forms.

Slow–fast oscillations are observed in numerous applications ranging from neuroscience and earth sciences to engineering. In this study, we perform a recurrence analysis of prototypical signals derived from well-established models, namely, the Van der Pol model, a modified form of the Izhikevich model, and the Hodgkin–Huxley model. First, we show a potential pitfall of phase space reconstruction, as the number of slow–fast regions could be exaggerated when the phase space is reconstructed by time delay embedding. We observe that the recurrence network, which represents the high-dimensional phase space of the underlying system, is clustered in certain regions and also exhibits protrusions, in addition to the closed-loop structure typically seen for periodic signals. We argue that these clustering and protrusions effects in the recurrence network arise due to the presence of slow and fast timescales in the system. Additionally, we can detect such features in micro-patterns along the diagonal lines of the corresponding recurrence plot. Finally, we observe similar features on the recurrence plots and recurrence networks of time series of signals acquired from experiments performed on a sub-scale liquid rocket combustor and a model gas turbine combustor during the state of thermoacoustic instability.

## I. INTRODUCTION

The rhythmic beating of the heart,^{1} periodic firing of neurons,^{2} spontaneous oscillations of chemical reactions,^{3} dangerous self-excited oscillations in suspension bridges,^{4} glacial oscillations,^{5} and high amplitude oscillations in aircraft engines and rocket engines^{6} are some examples of the various periodic phenomena we come across in our lives. Most of these phenomena exhibit oscillations at a preferred timescale known as the time period of the oscillation. However, many periodic phenomena are inherently made up of more than one timescale in an oscillation.^{7} Such periodic phenomena are popularly classified as slow–fast oscillations.^{7} Such systems are found across a wide range of applications ranging from medicine,^{8} economics,^{9} physical sciences,^{10} and earth sciences^{11} to engineering.^{12–14}

For example, let us consider the electrocardiogram (ECG) signal, wherein the electrical activity in the heart is recorded using a set of electrodes. A typical cycle of an ECG signal is defined by different processes such as atrial depolarization, ventricular depolarization, and ventricular repolarization.^{15} Each of these processes (designated as P waves, QRS complexes, and T waves in one cycle of the ECG signal) have an intrinsic timescale. A characteristic feature of slow–fast systems is that their periodic waveforms are radically different from those of harmonic oscillators. For the most simple case of a slow–fast system containing two timescales, a slow growth/decay is accompanied by a fast decay/growth. As a result, a slow–fast system could spend more time in the growth or decay phase. To present an example in electrical engineering, the charge and discharge of a capacitor^{16} are characterized by a slow and fast timescale, respectively. In a similar manner, the periodic stick-slip motion of a bowed violin string exhibits more than one timescale.^{17}

In the nonlinear dynamics literature, the term slow–fast systems have also been used to describe multiple timescales that cause periodic amplitude modulations, bursting oscillations, and mixed-mode oscillations.^{19} Hence, we illustrate the slow–fast systems that we discuss in this study along with other slow–fast systems in Fig. 1. Bursting oscillations [Fig. 1(a)] are characterized by epochs of large amplitude periodic oscillations followed by quiescence.^{20} Mixed-mode oscillations are another class of periodic oscillations that exhibit amplitude switching between two or more amplitude states in the signal. In periodically modulated waves [Fig. 1(b)], the amplitude envelope of the signal oscillates at a slow timescale over a fast oscillating signal. In these types of slow–fast systems, the slow timescale corresponds to the modulation of the envelope of the signal, while the fast timescale pertains to the high frequency oscillation in the signal.

However, unlike all these types of periodic oscillations, the slow–fast systems described in this study contain all the slow and fast timescales within one period of oscillation [Fig. 1(c)]. Such slow–fast systems have been long studied under the guise of relaxation oscillators. These oscillators are a class of limit cycle oscillators, which are characterized by a non-sinusoidal periodic waveform.^{21} Relaxation oscillations have been modeled using several models such as the Van der Pol oscillator,^{21} Fitz-Hugh-Nagumo oscillator,^{22} and LEGION.^{23}

Traditionally, slow–fast systems with a pre-established set of governing equations have been solved using conventional methods from the linear theory. A classical technique is to reduce the set of governing equations to the weak or the strong nonlinear limit,^{24} whenever the two timescales are widely separated. Then, the system of equations is solved to obtain the resultant amplitudes and phases of the signal. Apart from this method, various other techniques such as the perturbation theory, the method of multiple timescales, and the method of averaging exist.^{25} However, experimental and other real-world signals rarely have any well-defined functional forms, which can be solved using these methods. Moreover, the timescales in practical systems are seldom widely separated. All these obstacles render the analysis of such signals intractable. At this juncture, the framework of dynamical systems and complex systems theory offers a promising way to understand and characterize the dynamics of complex systems in man-made systems as well as those in nature.^{26,27}

In this study, we first characterize the dynamics of prototypical slow–fast signals obtained from well-established models, i.e., the Van der Pol (VDP) model, a modified form of the Izhikevich model, and the Hodgkin–Huxley model. We use nonlinear time series methods based on the recurrence analysis of the phase space trajectory such as the recurrence plot (RP)^{28} and the recurrence network (RN)^{29} to distinguish the properties of these signals. Following the same methodology, we analyze two high-dimensional slow–fast signals of thermoacoustic oscillations from experiments—the time series of unsteady heat release rate signals from a model gas turbine combustor^{30} and the acoustic pressure signal obtained from a model liquid rocket combustor.^{31}

The rest of the paper is outlined as follows. The methodology used in this study is briefly described in Sec. II. Here, we concisely detail the RP and RN construction techniques adopted in this study. Next, in Sec. III, we begin by analyzing low-dimensional systems and comparing with a single timescale signal (sine wave). Later, we extend the analysis for probing slow–fast dynamics in a high-dimensional model and real-world signals. We demonstrate the robustness of the obtained RNs for different embedding dimensions. We briefly discuss the quantitative features of the VDP model and the modified Izhikevich model. In the supplementary material, we describe the quantitative features of the other slow–fast signals used in this study. Finally, we summarize the key findings in Sec. IV.

## II. METHODOLOGY

The dynamics of nonlinear systems can be understood from the evolution of their phase space trajectories.^{32} Each state point on the phase space trajectory is described as a unique combination of system variables. However, in a physical system, it is almost impossible to acquire all the pertinent system variables required for the construction of the phase space. To circumvent this problem, we can reconstruct the phase space by time delay embedding.^{32} In this method, the phase space is realized by plotting the time series against its delayed versions in an appropriate dimensional space. The optimum time delay and the embedding dimension need to be selected prior to phase space reconstruction. The time delay ($\tau $) is selected such that the delayed vectors are independent of each other. We estimate $\tau $ using the autocorrelation function^{25} (ACF) and use the modified false nearest neighbor method, developed by Cao,^{33} to obtain the optimum embedding dimension ($d$).

Using this method of phase space reconstruction, we can visually unravel the dynamics of nonlinear systems from its phase space attractor only in low dimensions ($d\u22643$). However, a vast number of real-world signals tend to have higher dimensions ($d>3$). As a result, the fundamental property of recurrence of a phase space trajectory is exploited to understand the underlying hidden features of high-dimensional nonlinear systems.^{28,34}

### A. Recurrence plots

Recurrence of state points in the phase space is a fundamental property of deterministic signals. Recurrence plots (RP), introduced by Eckmann *et al.*,^{34} allow us to identify the time instants at which the trajectory of a system roughly revisits the same region in a $d$-dimensional phase space.

For a trajectory $x\u2192(t)$ made of $n$ time instances, the pairwise distances between state points in the phase space can be contained in a distance matrix ($Dij$), as formulated by

Here, $xi\u2192\u2212xj\u2192$ is the Euclidean distance between the two state points, $i$ and $j$, on the phase space trajectory. Next, the distance matrix is binarized by defining a threshold ($\u03f5$) to obtain the recurrence matrix^{35} ($Rij$),

where $\Theta $ is the Heaviside step function and $\u03f5$ is the threshold defining the neighborhood around the state point. Whenever the phase space trajectory recurs within the region defined by the $\u03f5$—size ball, we mark 1 in the recurrence matrix. Non-recurring points are marked by zeros in the recurrence matrix. In our RP, values of one and zero are designated as black and white points, respectively. Thus, the RP is a two-dimensional arrangement of black and white points that exhibits different patterns characterizing underlying dynamics of the signal. For a periodic signal of constant amplitude, we obtain uninterrupted diagonal lines in the RP. Patterns in RPs have garnered the attention of physicists in many instances.^{36} However, understanding such patterns in the RPs of slow–fast systems have not yet been probed, to the best of our knowledge.

One of the methods to select a recurrence threshold ($\u03f5$) is to fix the recurrence rate^{28,37} ($RR$). $RR$ is defined as $RR=1n2\u2211i,j=1nRij$. It estimates the percentage of recurring points in a RP. We observe that a value lower than the optimum $RR$ fails to completely capture the periodicity in the signal and is reflected as broken diagonal lines in the RP. In this study, an optimum value of $RR$ is selected after careful consideration for each slow–fast systems.

### B. Recurrence networks

Recurrence networks^{38} are a class of networks through which high-dimensional systems can be understood. A recurrence network (RN) comprises time instants as nodes, and their links are based on recurrences of state points in the phase space. Similar to the RPs, we can create a $\u03f5$—recurrence network,^{29,38,39} where $\u03f5$ is the optimal recurrence threshold. A value higher than the optimum value results in superfluous connections in the RN, distorting the network topology, whereas a lower $\u03f5$ would not capture the recurrence of trajectories in the phase space, leading to an underdeveloped network topology. The topology of the RN has been found to preserve the phase space of the high-dimensional nonlinear system.^{40}

To construct a RN, we require an adjacency matrix $A$ to be computed from the recurrence matrix $R$ for an $\u03f5$-threshold,

where $\delta $ is the identity matrix of the same size as $R$ and is used to remove self-connections. If the distance between the state space points is within the $\u03f5$-threshold, then, $Aij=1$, and the corresponding two nodes are connected. Otherwise, the two nodes remain disconnected, and $Aij=0$. Once the adjacency matrix is constructed for all pairs of nodes, several network measures can be computed from the RN, quantifying the geometrical structure of the phase space attractor.^{39,41} Using the network properties obtained from the RN, a number of studies have used RN to study the dynamical features of diverse systems.^{40,42–44}

We visualize the RN using the open-source network analysis platform, Gephi.^{45} The geometric feature of RN is attributed by a force directed algorithm known as “Force Atlas” in Gephi, where the connected nodes are attracted to each other, while the disconnected nodes are repelled from each other. An appropriate RN visualization is achieved when the forces are balanced, leading to a static RN. Each node in the RN is color-coded based on degree, a network property.^{46} The degree of a node $i$ ($Ki$) refers to the number of connections node $i$ has to all other nodes in the network and is calculated as

## III. RESULTS AND DISCUSSIONS

We progressively investigate the recurrence properties of slow–fast systems from low-dimensional systems to high-dimensional systems. As case studies for low-dimensional systems, we consider the Van der Pol (VDP) oscillator and the modified signal derived out of Izhikevich’s spiking neuron model.^{18} We consider the Hodgkin–Huxley model^{2,47} as a case for studying high-dimensional prototypical slow–fast signals. We, then, analyze the time series of heat release rate oscillations obtained from experiments in a laboratory-scale gas turbine-type turbulent combustor^{30} and the acoustic pressure signal from a laboratory-scale model multi-element liquid rocket combustor,^{31} during the state of thermoacoustic instability, to understand the recurrence dynamics of slow–fast signals in higher dimensional physical systems.

### A. Recurrence analysis of low-dimensional prototypical signals

Prior to understanding slow–fast systems, we analyze a harmonic signal, namely, a sine wave of unit amplitude and a time period of $2\pi $ [see Fig. 2(a)], which is definitely a single timescale system. The phase space of the sine wave is a circular loop structure [Fig. 2(b)], wherein the phase space trajectory evolves at a uniform speed. Here, the uniform speed of the phase space trajectory is attributed to successive state points on the trajectory separated by equal distances in the phase space. In the corresponding RP [Fig. 2(c)], we observe only equally spaced, non-interrupted diagonal lines with spacing equal to the time period of the oscillation. The corresponding RN topology of the sine wave [Fig. 2(d)] shows a circular loop filled up with same degree nodes.

Now, we start analyzing slow–fast systems where we first consider the VDP system,^{21} which is perhaps the most studied slow–fast system. Its governing equations are

where $\mu $ is referred to as the nonlinearity parameter to obtain relaxation type oscillations. We fix $\mu $ = 2 for the current analysis. The time series of variables, $x(t)$ and $y(t)$ of Eq. (5), are plotted in Fig. 3(a). The corresponding phase portrait exhibits a closed-loop, confirming the periodicity of the time series [Fig. 3(b)]. However, unlike the phase space of the harmonic signal [Fig. 2(b)], we observe that the phase space trajectory evolves at different speeds, giving rise to the slow and fast timescales. The separation between successive state points on the phase space trajectory during the fast epoch is large as compared to that of the slow epoch. As a result, the fast epoch can be visually discriminated from the slow epoch in the phase space. For the VDP system, we observe two epochs of slow oscillations (marked as $S$) and two epochs of fast oscillations (marked as $F$) within a cycle in the original phase space (i.e., a plot between the variables $x$ and $y$ of the system).

In Fig. 3(c), we show the time series of the variable $x(t)$ and its delayed copy $x(t+\tau )$. Here, the delay $\tau $ is obtained by the first zero crossing in the autocorrelation function (ACF). Unlike the original phase space [Fig. 3(b)], in the reconstructed phase space of $x$ [see Fig. 3(d)], we obtain four epochs of slow and fast oscillations (marked as $S$ and $F$, respectively) within a cycle of oscillation. This exercise shows that systems containing slow–fast timescales need to be interpreted carefully based on the technique of phase space reconstruction, since the number of slow–fast regions could be exaggerated with respect to that present in the original phase space. We remark that the reconstructed phase space trajectory of $y(t)$ evolves at a single timescale (not shown here for brevity) and hence does not exhibit any slow–fast features. Hence, we must be wary of slow–fast oscillations in practical scenarios going unnoticed when we are not tracking the appropriate system parameter.

Furthermore, we plot the RP and the corresponding RN for the VDP system from the original phase space and from the reconstructed phase space [see Fig. 4]. The recurrence matrix is constructed by fixing $RR=0.05$. For both RPs [shown in Figs. 4(a) and 4(c), respectively], at a first glance, we observe only diagonal lines, indicating periodic behavior of the system. However, in the corresponding zoomed view of the RP in Figs. 4(a) and 4(c), we identify the presence of momentary thick regions along the diagonal lines of a RP. We attribute these thick regions to slow epochs and the thin regions to the fast epochs in the evolution of the phase space trajectory. We refer to the presence of such distinct black patterns on the diagonal lines in a RP of the periodic signal as micro-patterns of RP. Thus, with the analysis of such micro-patterns, we can distinguish the time instances corresponding to slow regions from the fast regions in the phase space.

The reason behind the occurrence of such a micro-pattern in the RP can be understood from the evolution of the phase space trajectory at slow and fast timescales. When the phase space trajectory evolves at a slower rate in the phase space, it spends relatively more time within an $\u03f5$-threshold as compared to the phase space trajectory for the fast motion. This leads to the thickening of the diagonal lines in the RP. A similar argument can be given to explain the thinning of the diagonal lines whenever the phase space trajectory exhibits fast motion.

Recently, Kraemer and Marwan^{48} reported a tangential motion of phase space trajectories leading to uneven thickening along diagonal lines. They identified the temporal correlations (i.e., preceding and succeeding state points fall within the recurrence threshold) in the data, the presence of noise, and the usage of insufficient embedding dimension as reasons for this effect. Here, we observe thickening of diagonal lines in the RP for the slow epochs in the phase space of the prototypical slow–fast signal (with no noise) embedded using an optimum embedding dimension. As a result, we can attribute the thickening of diagonal lines with the temporal correlations in the slow epoch in the phase space of the slow–fast system.

We observe that the network topologies of both the original and the reconstructed VDP system are similar to the corresponding phase space observed in Figs. 3(b) and 3(d). Thenceforth, the nodes represented in the RNs are color-coded based on the increasing order of their respective degrees. In the corresponding RNs [see Figs. 4(b) and 4(d), respectively], we identify distinct regions that exhibit spatial clustering of high degree (red) nodes among the almost uniform distribution of the low degree (blue) nodes. The spatial clusters of high degree nodes within a cycle represent the region in which the trajectory moves slowly in the phase space, resulting in a higher number of connections in the RN. There are two such regions in the RN constructed from the original phase space [see Fig. 4(b)] and four slow regions in the RN from the reconstructed phase space [Fig. 4(d)], exactly matching their number in the respective phase spaces shown in Figs. 3(b) and 3(d).

Next, we consider another slow–fast signal [see Fig. 5(a)] obtained by modifying the time series of the variable $x$ from the Izhikevich’s spiking neuron model.^{18} First, we solve for the variable $x$ in the set of Eq. (6). The dimensionless parameters $a=0.1$, $b=0.2$, $c=\u221260$, $d=8$, and $I=110$ are used to obtain spiking behavior in $x(t)$,

The parameters $a$, $b$, $c$, and $d$ retain the same meaning as described originally in Izhikevich.^{18} Then, the resulting time series is modified so that enough number of points are present both during the growth and decay phase of the oscillations, to get a connected RN (refer to Sec. S-A in the supplementary material).

We observe that one oscillation in this signal is almost symmetric about the growth and the decay phase [Fig. 5(a)]. The three-dimensional phase space by using $\tau =102$ time steps (obtained from ACF) for this signal is visualized in Fig. 5(b). We find that the three-dimensional phase space attractor is stretched along the three axes, while maintaining a closed-loop structure in the evolution of the phase space trajectory for one cycle of oscillation.

We observe that the RP exhibits continuous equi-spaced diagonal lines, signifying the periodic dynamics of the signal [Fig. 5(c)]. Superimposed on this RP, the micro-patterns exhibit intricate features unique to this slow–fast system. Similar to the VDP system, the thickened portions of the diagonal line [see Fig. 5(c)] correspond to the slow motions in the phase space. A perpendicular line segment occurs amidst two thickened regions along the diagonal line in the RP, whenever the phase space trajectories traversing in opposite directions are spaced within the $\u03f5$-threshold. This is also confirmed by the fast motions of the phase space trajectory at the extremities (or corners) of the phase space in Fig. 5(b), where the phase space trajectory reverses its direction abruptly, akin to a hairpin turn in a mountainous road. Thenceforth, we shall refer to such abrupt reversal in the trajectory, leading to the formation of line segments perpendicular to the main diagonal line as the *hairpin trajectory*. It is important to emphasize that there is no such occurrence of two neighboring phase space trajectories traversing in opposite directions in the VDP system [see Figs. 3(c) and 3(d)]. As a result, we do not obtain any perpendicular lines in the RP of the VDP system.

In Fig. 5(d), the RN for this signal is plotted. Within one cycle, the trajectory is predominantly slow with many nodes having a very high number of connections (red and green). The fast regions in the phase space are present in the protrusions comprising nodes with a low degree (blue). In contrast to the RN of the VDP system, we observe that only nodes with high and medium degrees (red and green colors, respectively) occupy the ring-like structure. However, the protrusions on the ring are predominantly occupied by the nodes with a low degree (blue). This characteristic behavior must arise out of some fundamental difference in these two slow–fast systems, which is being reflected on their respective recurrence properties. Also, we identify that the micro-patterns in the RP and the RN for this prototypical signal are clearly different from the ones obtained for the VDP system.

### B. Recurrence analysis of high-dimensional prototypical signals

We also analyze the recurrence properties of the well-known Hodgkin–Huxley model, which exhibits slow–fast oscillations.^{2,47} It is represented by

Here, $V$ is the potential, $I$ is the current per unit area, and $gi$ is the maximum value of conductance where $i$ corresponds to either one of potassium ($K$), sodium ($Na$), or leak channel ($L$). The gating variables, $\alpha $ and $\beta $, control the activation and inactivation of their respective channels. Variables $m$, $n$, and $h$ are non-dimensional quantities associated with the potassium channel activation, sodium channel activation, and sodium channel inactivation, respectively. These variables acquire values between 0 and 1. In Eq. (7), the constant parameters used are $ENa=115mV$, $EK=\u221212mV$, $EL=10.6mV$, $gNa=120mS/cm2$, $gK=36mS/cm2$, $gL=0.3mS/cm2$, and $Cm=1\mu F/cm2$.

The corresponding steady state values for the gating variables, $\alpha $ and $\beta $, are related to the potential $V$ as

The set of equations is solved using Euler’s method. In Fig. 6(a), we plot the time series of the membrane potential, $V$, obtained for $I=10nA/cm2$ in Eq. (7).

We observe that $V$ exhibits a limit cycle behavior with slow–fast timescales. From the corresponding three-dimensional phase portrait in Fig. 6(b), we see that certain regions are slow (marked $S$), while others are fast (marked $F$). The corresponding RP also exhibits unique micro-patterns on top of the diagonal lines [see Fig. 6(c)]. The hairpin trajectory in its phase space renders a sword-like structure similar to the perpendicular lines observed in the RP of the prototypical spiky signal [Fig. 5(c)].

The corresponding RN exhibits a protrusion made up of high degree nodes and several clusters built of medium degree nodes on top of a ring of low degree nodes [Fig. 6(d)]. Here, hairpins consist of slow epochs. Hence, we obtain protrusions containing high degree nodes, in stark contrast to the RN of the modified Izhikevich model. Altogether, the RN of the Hodgkin–Huxley model contains both features observed in Figs. 4(b) and 4(d) and Fig. 5(d).

After analyzing the phase space dynamics and recurrence properties of these three prototypical slow–fast signals along with a sine wave, we understand that the RN for these slow–fast systems exhibits characteristic features on top of the closed-loop structure, expected for periodic signals. Moreover, the RP of such systems is manifested by unique micro-patterns pertaining to slow–fast dynamics over the diagonal lines.

### C. Recurrence analysis of high-dimensional experimental signals

In order to confirm the aforementioned observations in the slow–fast dynamics of real-world data, we present the results of the investigation of two different time series acquired from experiments in a laboratory scale gas turbine-type turbulent combustor and a model liquid rocket combustor during the state of an oscillatory instability, known as thermoacoustic instability.^{49} Here, thermoacoustic instability is a dynamical regime featured by large amplitude, self-sustained periodic oscillations in the acoustic pressure, $p(t)$, and the heat release rate, $q\u02d9(t)$, along with other dynamical variables of the system. The occurrence of this feedback-driven phenomenon overwhelms the thermal protection systems and compromises the controllability and structural stability of gas turbines and rocket engines.^{49,50}

First, in Figs. 7(a) and 7(b), we consider the time series of heat release rate oscillations ($q\u02d9(t)$) and the corresponding reconstructed phase space obtained during the state of thermoacoustic instability for a gas turbine-type turbulent combustor. We observe that the time series is spikier than a sine wave, exhibiting a clear departure from sinusoidal signals. The spikiness in the signal^{30} [Fig. 7(a)] is attributed to the near instantaneous heat release rate as a result of the impingement of the large-scale coherent vortex structure carrying fuel–air mixture against the walls of the combustor.^{51}

In the corresponding phase space of the heat release rate ($q\u02d9$) signal in Fig. 7(b), we observe a distorted closed-loop structure, indicative of the non-uniform evolution of the phase space trajectory due to the presence of slow and fast timescales. However, such slow and fast timescales are not too separated when compared to the earlier phase space of prototypical signals. In the RP of this signal [see Fig. 7(c)], we see the presence of continuous diagonal lines, indicating sustained periodicity in the oscillations. The corrugations along the diagonal lines arise due to the presence of the slow and fast timescales in the phase space. The RN for this signal [Fig. 7(d)] looks similar to that of VDP as there are clusters of high degree nodes (red) on the ring of medium degree (green) nodes. The clusters pertain to the slow regions in the phase space.

Finally, we investigate the time series of acoustic pressure oscillations ($p\u2032(t)$) in a multi-element model liquid rocket combustor during the state of thermoacoustic instability^{31} [Fig. 8(a)]. We observe that a major portion of the cycle is spent in the slow relaxation phase with a momentary jump in the pressure due to the fast compression phase of the signal. Physically, due to an increase in the speed of sound during the compression phase, the compression side catches up with the expansion side of the pressure wave. This phenomenon, known as wave steepening, results in an abrupt increase in the amplitude of the pressure oscillations.^{52} Under favorable conditions, the steepened wave manifests as a shock wave in the flow-field. Such wave steepened shock waves are commonly observed in the pressure oscillations in the combustion chambers of rockets.^{53,54}

The reconstructed phase space of this pressure signal, shown in Fig. 8(b), is similar to the phase portrait shown in Figs. 5(b) and 6(b), wherein the trajectory moves along the three axes to complete one oscillation cycle. Unlike the usual closed-loop structure of the phase space trajectory of periodic signals observed in the previous slow–fast systems, the phase space of the pressure signal exhibits a peculiar shape like a trefoil-knot.^{55} The geometrical difference of the phase space attractor is attributed to the vast divergence in the slow and fast timescales in the rocket system.

The RP and the RN for the pressure oscillations are plotted in Figs. 8(c) and 8(d), respectively. The RP contains unique micro-patterns arising due to the presence of spikes. On top of the diagonal lines that indicate periodicity of the signal, we observe thick regions divided by a thin region. The thick regions emerge due to the increased trapping of the phase space trajectory in the slow epoch, while the thin regions correspond to the fast spike in the phase space. Due to the hairpin trajectories at the extremities in the reconstructed phase space, the RP of this signal exhibits line segments protruding away from each diagonal line in a periodic manner.

The RN of this signal looks similar to Fig. 5(d), based on its topological similarity. The protrusions in the RN in Fig. 8(d) are made up of high degree (red) nodes, unlike the RN in Fig. 5(d) where the protrusions are made up of low degree (blue) nodes. We believe that suitable measures derived from the RN can be used to benchmark simulations that predict thermoacoustic instability in combustion chambers of liquid rocket engines.

### D. Effect of embedding dimension on the recurrence network of a slow–fast system

We demonstrate the robustness of the topology of the RNs for different embedding dimensions ($d$) for the modified Izhikevich model [Fig. 9(a)] and time series of the acoustic pressure oscillations acquired during the state of thermoacoustic instability in a liquid rocket combustor [Fig. 9(b)]. The corresponding embedding dimensions selected from the modified false nearest neighbors method for these two cases are $d=6$ and $d=10$, respectively.

In general, we observe that the RNs exhibit closed-loop structures characteristic of periodic orbits for the range of $d$ shown. For both cases, the realized RNs for lower embedding dimensions are distorted. With a further increase in $d$, the topology of RN converges and remains largely the same for further increase in $d$. In other words, for higher $d$, we find characteristic features such as the number of protrusions and clustering of nodes to be nearly the same with increasing $d$. However, the topology of RN converges at an earlier $d$ for experimental data, than that estimated by the modified false nearest neighbors method. This observed change in the optimal $d$ from the RN and from the modified false nearest neighbors method is not seen in the case of the prototypical signals. We believe that the presence of noise in the experimental data leads to this deviation. The ability of RN to capture the features of the high-dimensional phase space in slow–fast systems can help us to cross-verify the optimum embedding dimension. A similar variation is confirmed for all the other slow–fast systems discussed here and the respective plots are shown in Sec. S-B of the supplementary material.

### E. Quantitative analysis of recurrence network properties

Next, we quantitatively characterize the recurrence network topology obtained from the reconstructed phase space to unravel the differences between the slow–fast system and a single timescale system such as a sine wave. The conventional approach in networks built from time series is to compute the global network measures such as the mean degree of a network.^{46}

To ensure that no disparity arises due to the length of the time series and the frequency of oscillations, we ensure that both the sine wave and the VDP signals are of the same frequency and amplitude. We find that the mean degree ($Kmean=\u2211i=1nKi/n$) remains the same for both the sine wave and the VDP model ($Kmean=152$). In Fig. 10, we plot the histogram of the probability distribution function of the $Ki$ of each node in the RN of sine waves and the RN of the VDP model. We observe that the distribution for the VDP model has a broader spread compared to a unique value for the sine wave, justifying the presence of multiple timescales in spite of its periodic behavior of both models. This is confirmed for all the other slow–fast systems discussed here, and the respective plots are shown in Sec. S-C of the supplementary material.

To unravel the variation of degree of each node ($Ki$) in the RN, we show the variation of $Ki$ along with the distance ($PDi$) between consecutive points ($x\u2192$) in the reconstructed phase space for the VDP model and sine waves in Fig. 11. Here, the nodes are labeled according to their temporal appearance in their corresponding phase spaces. $PDi$ is calculated using the Euclidean norm as

For every cycle of oscillation in the sine wave, we see that both $PDi$ [Fig. 11(b)] and $Ki$ [Fig. 11(c)] remain invariant. However, for every cycle of oscillation in the VDP model, we observe four oscillations in both $PDi$ and $Ki$, since there are four slow epochs and four fast epochs in the phase space of the reconstructed signal of the VDP model [see Fig. 3(d)]. Moreover, we see the simultaneous occurrence of lower values in $PDi$ and higher values in $Ki$, whenever the phase space exhibits slow motion. During slow motion, the consecutive points in the phase space are located nearby, and hence, $PDi$ is low and $Ki$ is high. Correspondingly, for epochs of fast motion in the phase space, we obtain higher values in $PDi$ and lower values in $Ki$. The four oscillations within a cycle of oscillation manifest as four clusters in the RN, as seen in Fig. 4(d).

In Fig. 12, we plot the variation of $PDi$ and $Ki$ for the time series obtained from the modified Izhikevich ($M.Iz$) model and a sine wave, both of which contain periodic oscillations with a time period of 10.11 s sampled at a temporal step size of 0.01 s. Similarly, we observe significant temporal variations in both $PDi$ and $Ki$ for the slow–fast system, while that of the sine wave remains constant. Also, we find higher values in $Ki$, whenever $PDi$ is low and vice versa. There are five oscillations in $PDi$ and $Ki$ within a cycle of oscillation of the prototypical slow–fast signal. This manifests as five protrusions of low degree nodes in the RN shown in Fig. 5(d). For the other slow–fast signals discussed in this study, we show the corresponding plots in Sec. S-C of the supplementary material.

With the understanding gained from analyzing the various prototypical and experimental systems in this study, we noticed that the dynamics of slow–fast systems can be understood based on their recurrence properties. From the RNs, we observed that some slow–fast systems exhibit protrusions, while other systems display clustering. Each slow–fast system imparts a signature micro-pattern over the diagonal lines in their corresponding RP. It is also interesting to note that even though both real-world systems discussed here operate in a regime of thermoacoustic instability, both systems exhibit different RN topologies due to a difference in their underlying mechanisms that generate such oscillations. Finally, we interpret the RN topology using the temporal variation in the distance between consecutive points in the phase space and the degree of each node in the RN of slow–fast systems.

## IV. CONCLUSIONS

In this study, for the first time, recurrence properties of slow–fast systems are studied by means of recurrence plots and recurrence networks. A systematic approach is adopted by first performing the analysis on prototypical signals before analyzing high-dimensional signals obtained from experiments. We find that slow–fast systems exhibit different recurrence properties compared to periodic systems that operate on a single timescale. We observe that unique features about the slow–fast system can be obtained from the micro-patterns along the diagonal line in the RPs, unlike mere straight lines observed in the RP for harmonic signals. Particularly, we find that hairpin trajectories in the phase space lead to the occurrence of line segments perpendicular to the main diagonal line in the RP. These findings help improve the understanding of the various patterns evident in the RP. Furthermore, we identify characteristic features in the corresponding RN topologies for slow–fast systems. In addition to the closed-ring structure of periodic signals, we also observe protrusions and clustering in the RN for slow–fast systems. Such additional features in the RN result from the temporal variation in the distance between consecutive points in the phase space and the degree of nodes in the RN of slow–fast systems. Such variations are absent in single timescale systems. We believe that this study will have wide ranging applications to understand the dynamics of various diverse systems across natural sciences, medicine, econometrics, and engineering.

## SUPPLEMENTARY MATERIAL

See Sec. S-A in the supplementary material for a detailed description of the need as well as the steps involved in the modification of the Izhikevich’s spiking model output used in this study. In Sec. S-B, we show the effect of embedding dimensions on the topology of the recurrence network for the Van der Pol model, the Hodgkin–Huxley model, and the heat release rate oscillations acquired from the gas turbine-type turbulent combustor. In Sec. S-C, we describe the temporal variation of the distance between consecutive points in the phase space and their corresponding degree in the recurrence network for the modified Izhikevich model, the Hodgkin–Huxley model, heat release rate oscillations acquired from the gas turbine-type turbulent combustor, and acoustic pressure oscillations from the laboratory-scale rocket combustor.

## ACKNOWLEDGMENTS

The authors from IIT Madras gratefully acknowledge the funding provided by Air Force Office of Scientific Research (AFOSR) under Award No. FA2386-18-1-4116 (Grant Program Manager: Lt. Col. Winder). R.I.S. is thankful to Dr. Sankaran (AFRL) for his interest and helpful discussions. We also thank Ms. Godavarthi (UCLA) for helpful discussions on the visualization of the recurrence network. P.K. and I.P. thank MHRD-India and IIT Madras for the support provided with the research assistantship.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*2018 AIAA Aerospace Sciences Meeting*(AIAA, 2018), p. 1185.

*Turbulence, Strange Attractors, and Chaos*(World Scientific,

*Third International AAAI Conference on Weblogs and Social Media*(The AAAI Press, 2009).