Recurrence Network Analysis of Exoplanetary Observables

TL;DR

This paper applies recurrence network analysis to exoplanetary data, successfully identifying different dynamical regimes in synthetic and real observations, thus providing a new approach to studying planetary dynamics.

Contribution

It introduces a novel application of recurrence network analysis to exoplanetary data, bridging synthetic models and real observations for dynamical regime detection.

Findings

Different dynamical regimes identified in synthetic data

Analysis of real astronomical observations aligns with previous studies

Recurrence network measures effectively characterize planetary dynamics

Abstract

Recent advancements of complex network representation among several disciplines motivated the investigation of exoplanetary dynamics by means of recurrence networks. We are able to recover different dynamical regimes by means of various network measures obtained from synthetic time series of a model planetary system. The framework of complex networks is also applied to real astronomical observations acquired by recent state-of-the-art surveys. The outcome of the analysis is consistent with earlier studies opening new directions to investigate planetary dynamics.

Tables

Table 1. Table 1: The embedding parameters of various time series in SJS system used to construct RNs. Initial conditions in Fig  1 (a) are considered with two different length, the original 950 data point signals and longer one containing 3500 measurements. The time delay parameters d , 𝑑 d, τ , 𝜏 \tau, and ϵ italic-ϵ \epsilon are the embedding dimension, time lag, and threshold, respectively. The threshold in each case is obtained from recurrence matrix with fixed R ​ R = 0.1 . 𝑅 𝑅 0.1 RR=0.1.

		Regular		Chaotic		Resonant
Data		$950$	$3500$	$950$	$3500$	$950$	$3500$
RV-Syn333Synthetic data are obtained directly from numerical integration.	$d$	$4$	$6$	$10$	$9$	$5$	$5$
	$\tau$	$2$	$2$	$3$	$2$	$3$	$2$
	$ϵ$	$2.1e-4$	$1.84e-4$	$3.3e-4$	$2.5e-4$	$1.5e-4$	$2.1e-4$
RV-Spl444Spline interpolation has been made on synthetic (Syn) time series after adding noise and removing some data points as described in the text.	$d$	$10$	$11$	$14$	$13$	$11$	$11$
	$\tau$	$3$	$5$	$3$	$2$	$3$	$3$
	$ϵ$	$5.65e-6$	$4.8e-6$	$5.72e-6$	$5.0e-6$	$5.7e-6$	$4.6e-6$
TTV-Syn	$d$	$8$	$9$	$18$	$17$	$7$	$9$
	$\tau$	$4$	$4$	$5$	$4$	$10$	$19$
	$ϵ$	$251.8$	$292.8$	$4922.05$	$6814.92$	$340.35$	$1019.59$
TTV-Spl	$d$	$9$	$12$	$15$	$16$	$13$	$17$
	$\tau$	$3$	$3$	$12$	$14$	$11$	$16$
	$ϵ$	$34.05$	$53.51$	$634.87$	$987.25$	$47.38$	$374.74$

Equations

$$A_{i,j}=A_{i,j}(ϵ)=R_{i,j}(ϵ)-{\delta }_{ij},$$

$$\Delta P(\mathcal{M},x)=P(\mathcal{M})-P(x),$$

$$\begin{array}{rl}{S}_1({\mathcal{M}}^{\ast })& :=\{(a,e)|\mathcal{M}(a,e)\le {\mathcal{M}}^{\ast }\},\\ S_2({\mathcal{M}}^{\ast })& :=\{(a,e)|\mathcal{M}(a,e)>{\mathcal{M}}^{\ast }\},\end{array}$$

$$\begin{array}{rl}{S}_1^{\prime }(Q_{\alpha }(x))& :=\{(a,e)|x(a,e)\le Q_{\alpha }(x)\},\\ S_2^{\prime }(Q_{\alpha }(x))& :=\{(a,e)|x(a,e)>Q_{\alpha }(x)\},\end{array}$$

Full text

††preprint: AIP/123-QED

The dynamical variety of known extrasolar planetary systems stimulates their stability analysis as an important characteristic of the system. The common procedure to perform such an investigation is to obtain the parameters of the best-fit planetary model and then integrate the equations of motion numerically. In this work we propose a method in order to describe system dynamics that is based on the observed uni-variate time series and the topology of complex network reconstructed from the signal itself. We show that the procedure is capable to distinguish the ordered and chaotic motion in a synthetic 2-planet system with a given level of significance. We also find that different kind of data sets, produced by the currently used detection methods, works well within the framework of the method. The analysis of real-world observations also provides results consistent with former studies. The technique is computationally very efficient since it does not require the phase space trajectories and, therefore, costly n-body simulation can be avoided. Thus this new strategy can be used as a complementary tool to extract the dynamical behavior of extrasolar planetary systems.

I Introduction

Dynamical stability is one of the most relevant physical characteristics (beside the bulk density, atmospheric composition, interiors of a planet, etc.) to describe extrasolar planetary systems. In order to study the dynamical evolution of a system the initial conditions of the numerical integration are essential. In planetary science this input includes the orbital elements for a given epoch. Nowadays the ground and space-based observations provide a large amount of extremely precise radial velocity (RV) and transit timing data. These information can be transformed by using comprehensive statistical methods Foreman-Mackey et al. (2013) into the desired initial conditions.

The majority of known exoplanetary systems harbour more than one planet resulting in complex non-Keplerian dynamics Fabrycky (2010). The possible instability of the system is an internal attribute of the deterministic dynamics when three or more celestial objects participate in the motion. To explore the chaotic nature one has to have the integrated trajectories and/or the structure of the relevant part of the phase space assisted by a chaos detection method (Lyapunov exponent, MEGNO ($\mathcal{M}$) – Mean Exponential Growth of Nearby Orbits Cincotta and Simó (2000), etc.). These requirements can be fulfilled by having the best-fit planetary model incorporating the orbital elements and the masses as it has already been shown Rivera and Lissauer (2001); Armstrong et al. (2015); Batygin, Deck, and Holman (2015). The results of these studies indicate the heterogeneity of exoplanetary dynamics. Interestingly, there are also candidates with chaotic properties Deck et al. (2012); Panichi, Goździewski, and Turchetti (2017) beside the larger number of stable resonant configurations.

We propose an alternative method to perform stability analysis of exoplanetary systems that requires only a scalar time series of the measurements, e.g. RV, transit timing variation (TTV) 111As a consequence of the mutual gravitational interaction of the planets one might observe slight deviations in orbital periods called transit timing variation., or astrometric positions. The fundamental concept of Poincaré recurrences in closed Hamiltonian systems and the powerful techniques of nonlinear time series analysis Kantz and Schreiber (2003) combined with complex network representation Donner et al. (2010) allow us to investigate the underlying dynamics without having the equations of motion. That is we can ignore the orbital elements as initial conditions and other system parameters from the entire analysis. Moreover, the procedure completely disregards the use of numerical n-body integration which is an essential and sometimes fairly time consuming part of the stability analysis. We present that this new scheme, first time applied in exoplanetary research, works well with signals acquired by the currently available observational techniques (see above).

II Method

A scalar time series carries relevant dynamical information if the measured physical quantity is coupled to other state space variables of the system. Based on this fact Takens’ embedding theorem Takens (1981) ensures that the phase space trajectory can be reconstructed from the measured data. The procedure is called time delay embedding and requires two parameters, the time lag ($\tau$) and the embedding dimension ($d$). There are standard techniques to find the appropriate values of these parameters, e.g. false nearest neighbours help to estimate the embedding dimension Kantz and Schreiber (2003) while the first minimum of the mutual information function yields the time delay Semmlow and Griffel (2014).

Once the reconstructed trajectory is ready to use, one can search for Poincaré recurrences, the events when the trajectory $\mathbf{x}$ comes back at time $t_{j}$ into an $\epsilon$-ball around its earlier position at time $t_{i}$ ($t_{i}<t_{j}$). The idea of recurrence plot (RP) Eckmann, Oliffson Kamphorst, and Ruelle (1987) is a 2D binary matrix representation of these events. More precisely, the matrix element $R_{i,j}(\epsilon)=1$ if $||\mathbf{x_{j}}-\mathbf{x_{i}}||<\epsilon,$ and 0 otherwise. The $||.||$ can denote various norms such as Euclidean, Maximum, Manhattan. In this work the Maximum norm is used through the whole analysis. The $\mathbf{R}$ matrix is symmetric by its nature. First Refs. Zbilut and Webber (1992); Webber and Zbilut (1994) 222The visual interpretation of the matrix $\mathbf{R}$ can achieved by plotting a dot when the matrix element is 1 and leave it empty otherwise. proposed the method of recurrence quantification analysis (RQA) which identifies various measures based on the diagonal and horizontal/vertical texture of an RP. The advantage of RQA in dynamical analysis has been demonstrated in many applications from climatology through neuroscience to astrophysics. The idea of recurrence networks (RN) Donner et al. (2010, 2011) revolutionized the RP representation of dynamical systems using the fact that an RP can be thought of as the adjacency matrix of a complex network embedded in phase space. Quantitatively

[TABLE]

where $A_{i,j}$ is the adjacency matrix, $R_{i,j}$ the recurrence plot, and $\delta_{ij}$ the Kronecker delta. The role of $\delta_{ij}$ is to exclude the loops. The nodes of these graphs are the state space vectors that are connected by edges if their proximity is smaller than a threshold $\epsilon.$ The main advantage of the complex network framework is that the temporal correlations can be avoided since the dynamical feature of the underlying system is preserved via its topology, thus, the explicit time ordering is not relevant. The standard quantities known in graph theory (such as degree centrality, average path length, transitivity, etc.) are also applicable to RNs and they are capable to describe different dynamical regimes and dynamical transitions Zou et al. (2019). This fact will be utilized in the following.

III Results

III.1 The model system

First a well-defined model is investigated in order to present the reliability of the method. Consider the Sun-Jupiter-Saturn (SJS) two-planet system in our own Solar System with their actual orbital elements. One can explore, for example with the REBOUND package Rein and Tamayo (2015), the stability of the system by varying, say, Saturn’s initial semimajor axis ($a_{\mathrm{Sat}}$) and eccentricity ($e_{\mathrm{Sat}}$) in a grid and quantify the dynamics. Figure 1(a) shows the result of dynamical stability of the system in the ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) initial condition plane for 100x100 different trajectories characterized by the chaos indicator MEGNO ($\mathcal{M}$). The values of $\mathcal{M}\approx 2$ indicate stable motion while larger values (in this case truncated at 8) represent chaotic dynamics. One can observe that for higher eccentricities chaotic nature dominates the phase space. Nevertheless, there are tight green lobes penetrating into the red area. These formations are the MMRs between the two planets. The most remarkable resonance, 5:2, around $a_{\mathrm{Sat}}\approx 9.6$ extends up to very high eccentricities. The stability map of the SJS system is based on numerical integration that cover ca. 1000 orbital periods of Jupiter for each initial condition pairs. This astronomically fairly short time frame has been chosen in order to mimic real observations. The compact extrasolar planets, for instance the TRAPPIST-1 Gillon et al. (2017), Kepler-412 Deleuil et al. (2014), make thousands of revolutions on human time scales producing time series with desired length to be analyzed: $\sim$1000 periods of the innermost planet.

In order to apply recurrence network analysis (RNA) to SJS system, synthetic scalar time series were acquired from the full phase space trajectories in the following way. Radial velocity data is the x-component of the Sun’s velocity vector while transits are taken when the planet crosses the positive x-axis imitating the position of a real occultation. In case of the RV signals the trajectories were integrated until 1050 Jupiter orbits and 950 data points were sampled in order to set the sampling frequency larger than the mean motion of the Jupiter. For transits exactly 950 events have been taken. This also implies that the trajectories examined via transits have somewhat different length, namely, 950 Jupiter transits might take shorter or longer time in different dynamical regimes.

Figure 1(b) and (c) depict the network measures average path length ($\mathcal{L}$) for TTV signals and transitivity ($\mathcal{T}$) for RV data, respectively, in the ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) parameter plane. The phase space reconstruction and the network analysis has been done by the TISEAN Hegger, Kantz, and Schreiber (1999) and PYUNICORN Donges et al. (2015) packages, respectively. The similarity is evident between the panels (a) and (b), (c). In panel (b) for TTV signals the larger the $\mathcal{L}$, the more chaotic the dynamics. This result is in good agreement with general findings in discrete Hamiltonian systems like the standard mapZou et al. (2016). The RNA of RV time series produces similar results for the measure $\mathcal{T}$, panel (c). In this case transitivity closer to 1 yields weaker instability according to the general picture that regular trajectories show smaller divergence in phase space. This fact results in a more coherent topology of the network that produces higher transitivity. One can also notice that the stability map of TTVs is more sensitive to the edge of chaos inducing more details between the green and red domains. This is also important because weak chaos, also known as stickiness Contopoulos (2002) in Hamiltonian dynamics, emerges just at the edge of stability islands. The stickiness effect was also identified by using complex network measures in low dimensional conservative systems Zou et al. (2016). The method of RNA applied to single variable time series is, thus, able to reconstruct the stability map what is obtained from numerical integration and phase space trajectories. Moreover, it is also remarkable that achieving this outcome the length of the measured data must not be longer than that of the phase space trajectory.

Table 1 contains the embedding parameters for the reference trajectories introduced in Figure 1. In addition, time delay embedding has been done for longer time series containing 3500 data points in all three cases. This analysis verifies that 950 RV and TTV measurements provide a reliable RNA. In order to demonstrate that RN measures ($\mathcal{L},\mathcal{T}$) are suitable to distinguish regular and chaotic orbits, further statistical analyses are performed.

First, a point-wise difference based on the cumulative distribution functions (CDF) of $\mathcal{M}$ and ($\mathcal{L},\,\mathcal{T}$) in the ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) parameter plane is designed Zou et al. (2010)

[TABLE]

where $P(x),$ $x\in\{\mathcal{L}_{\mathrm{RV}},\mathcal{T}_{\mathrm{RV}},\mathcal{L}_{\mathrm{TTV}},\mathcal{T}_{\mathrm{TTV}}\},$ is the corresponding value of CPD at each combinations of 10000 ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) pairs. An example of CDF differences, $\Delta P(\mathcal{M},\mathcal{T}),$ for RV data is depicted in Figure 2. The difference is close to zero almost in the entire parameter plane. Remarkable deviation from zero can be observed in chaotic region.

Another more sophisticated quantitative comparison of RN measures can be obtained as follows. Let us define two disjoint subsets of MEGNO distribution defined by a critical value of $\mathcal{M}^{*}$

[TABLE]

with group size $n_{1}$ and $n_{2}=n-n_{1}$ ($n=10000$), respectively, see Figure 3(a). Then, take the $\alpha$-quantile $Q_{\alpha}(\mathcal{M})$ with $\alpha=n_{1}/n$ of the distribution $\mathcal{M}$ for a given $\mathcal{M}^{*}$ and construct the same division for the corresponding ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) parameters based on various RN measures

[TABLE]

where $x$ again incorporates certain network characteristics. This analysis allows one to compare the difference of two distributions based on distinct measures. Figure 3(b) and (c) depict one particular example of two subsets of ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$) pairs associated to $S_{1}(\mathcal{M}^{*}=8)$ and $S^{\prime}_{1}(Q_{\alpha}(\mathcal{T}_{\mathrm{RV}})),$ $\alpha\approx 0.57.$

The relative frequency $p$ of those parameter pairs that do not belong to the same group based on the two different measures indicate the false detection of dynamical nature. Slightly varying the chaos indicator close to the border of regular and chaotic feature, $2\leq\mathcal{M}^{*}\leq 8,$ the relative frequency of ”grouping errors” can be quantified. Figure 4 collects the frequencies of the average path length and transitivity obtained from RV and TTV signals representing that the classification error remains under 10% in all cases.

Although, the RN measures give qualitatively reasonable success for either observables, it is clear that their values depend on the choice of the threshold $\epsilon$. To avoid this weakness two possibilities are known. First, one can derive dynamical invariants (e.g. Rényi entropy, correlation dimension, maximal Lyapunov exponent) from recurrence plots. However, to get reliable feedback about the dynamics by these invariants one has to have much longer time series. And this is against our will, to make stability analysis based on real observables that contain only several hundreds/thousands of data points. The second option is to perform hypothesis tests with surrogate time series Schreiber and Schmitz (2000); Luo, Nakamura, and Small (2005). Since in planetary dynamics the signals show predominantly quasi-periodic variations, the null hypothesis in surrogation method we want to test should be that the original data comes from quasi-periodic process. Pseudo-periodic twin surrogates (PPTS) Carrión and Miralles (2016) combine the power of Pseudo-periodic Small, Yu, and Harrison (2001) and Twin surrogate Thiel et al. (2006) methods and fit perfectly to test quasi-periodicity in planetary dynamics.

More concretely, time series generated by the PPTS method provide different RP matrices for periodic/quasi-periodic orbits and the chaotic ones preserving the phase space structure. This ensures that these surrogates are relevant to test the null hypothesis that the original signal comes from quasi-periodic process. The question, however, still remains open: What are the dynamical properties of the underlying system in case of rejecting the null hypothesis. In view of the fact that planetary dynamics is deterministic, high dimensional and definitely nonlinear and by declining the null hypothesis of PPTS, meaning that the motion is not quasi-periodicity, the dynamics should be chaotic according to certain confidence level.

The hypothesis tests in Figure 5 comprise also the robustness of the RNA against measurement noise and missing data points. In order to reproduce a realistic astronomical observation, first Gaussian white noise with zero mean and unit variance is superimposed onto the synthetic time series in the SJS system with Signal-to-Noise ratio $\approx$6-7. Then, in a systematic way a random point is chosen along with an interval around this point which with random length is deleted from the data. This procedure is repeated until the desired amount of information (15%) is missing from the original signal. Based on Lekscha and Donner (2018), spline interpolation is performed on the unevenly spaced data before the time delay embedding and RNA is applied.

To obtain a 1% level of significance for a one-sided test 99 surrogate time series must be generated Kantz and Schreiber (2003) and then the RN measures calculated from the original time series has to be compared to those acquired from the surrogates. Usually we do not know whether the distribution of RN measures is Gaussian, therefore, a rank-based statistics is preferred instead of a normal. A rank-based statistics of the RNA specifies either to keep or reject the null hypothesis. Figure 5 illustrates the hypothesis tests for three different kind of dynamics marked by blue triangles in Fig. 1(a), and two observables (RV, TTV) in the SJS system. Let us concentrate on the first column, panels (a) and (d). The analysis here is devoted to show that the leftmost blue triangle ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$)=(7.2,0.02) corresponds to regular quasi-periodic motion. The red solid lines represent the RN measures $\mathcal{L}$ and $\mathcal{T}$ gained from the synthetic time series while the blue triangles symbolize the same RN measures based on the generated 100 PPTSs. Both plots indicate that $\mathcal{L}$ and $\mathcal{T}$ associated to the original time series fall into the ensemble of blue triangles. Thus, one can keep the null hypothesis, i.e. the original signal comes from quasi-periodic dynamics.

Due to strong gravitational perturbation taking place primarily at large eccentricities the system is destroyed, i.e. one of the planets escapes. Therefore, an initial condition close to the border of the regular part has been considered where the motion is chaotic while bounded for the integration time, ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$)=(8.0,0.2). Executing the RNA one finds the results portrayed in panels (b) and (e). The scheme is the same as before, in turn, it can be clearly seen that the red lines are located outside the blue triangle zoo suggesting that the null hypothesis can be rejected. That is, the original time series was generated by chaotic dynamics. For the RNA we pick up a third, resonant, initial condition in the stability map at position ($a_{\mathrm{Sat}},e_{\mathrm{Sat}}$)=(9.6,0.4). The reassuring results of the hypothesis tests are summarized in panels (c) and (f).

In order to verify that the spline interpolation does not cause any artificial effect during time delay embedding, the same analysis has been made on synthetic, i.e. noiseless uniformly sampled, signals as well. Figure 6 presents the outcome of hypothesis tests that perfectly match those based on noisy and scanty analysis in Figure 5.

III.2 Real-world measurements

Recently the number of known extrasolar planetary systems outstandingly increased due to the cutting-edge technology. Extremely precise light curves with transit timing measurements might shed light on the dynamical diversity of the observed systems as shown above. The catalog of the full-cadence data set of the Kepler mission Holczer et al. (2016) contains a large number of dynamically interesting systems especially from the TTV point of view. Some of them have been extensively studied by means of stability Deck et al. (2012); Panichi, Goździewski, and Turchetti (2017). One of these systems is Kepler-36 Carter et al. (2012) wherein two planets (Kepler-36b and c) orbiting the central Sun-like star nearly in 7:6 MMR. The orbital separation of the two ”Super-Earths” is fairly small, 0.013 au. Due to the tightly packed configuration, the mutual gravitational interaction causes large TTV, see Figure 7. The irregular dynamics of the system was proposed first by Deck et al. Deck et al. (2012), further analysis revealed Panichi, Goździewski, and Turchetti (2017) that the system might be close to the border of the 7:6 MMR and the dynamics is governed by the stickiness of the resonance.

Recurrence network analysis has been carried out for the TTV data of Kepler36b and c. It should also be noted that the data points in Kepler-36 time series are less with a factor of 8-10 than those in the SJS system. That is, the stability analysis is based on $\sim$100 orbital periods of the planets. Moreover, the number of missing transit points is somewhat larger than in the synthetic SJS, roughly 13% for inner planet and 25% for the outer one. The 99% significance level hypothesis test involving RN measures $\mathcal{T}$ and $\mathcal{L}$ is presented in Figure 8 and 9. Results depicted in Fig. 8 lower panel show that in the case of the more massive outer planet the dynamics is regular. In contrast, for the inner planet, upper panel, the null hypothesis can be rejected in accordance with the rank-based statistics described above. The outcome of the RNA implies that the system’s behaviour is irregular. However, average path length for the same system in Fig. 9 stipulates regular dynamics for both planets. Comparing the results with Figure 1, the statistical description of Kepler-36 fits to the picture of stable chaos appearing at the border of the resonances. Moreover, it also supports the view published in Deck et al. (2012); Panichi, Goździewski, and Turchetti (2017) based on different methods of dynamical analysis.

IV Conclusion

We conclude that the RNA method is suitable to reproduce the dynamical behaviour of a planetary system with reliable significance. This result is based on the fact that the topology of recurrence networks preserves the underlying dynamics of the system producing the time series under study. Thus, the RN measures, $\mathcal{L}$ and $\mathcal{T},$ are adequate to distinguish regular and chaotic nature of the motion. However, one should keep in mind that these measures are geometric characteristics rather than dynamical invariants of the dynamics. Therefore, the above mentioned contrast between quasi-periodic and chaotic behavior can only be made by appropriate hypothesis tests.

Further advantage of the present idea is that it uses directly the measured data set and requires neither Monte-Carlo simulation to achieve the best-fitting planetary model nor costly n-body integration. Consequently, the operation needs substantially shorter time to achieve the result. For example, the analysis of a time series with 950 data points in addition with the 100 surrogates requires no more than 10 minutes on a medium desktop machine.

We emphasize that 950 measurements for both signals (RV, TTV) seem to be a sufficient amount of data to restore the dynamical operation given by the same length of numerical integration. Extending the data set with the length of the integration time makes change to the embedding parameters, nevertheless, the recurrence network analysis provides the same results for longer signals. It has to be noted that this method depends on the measured time series which is the past of the planetary dynamics and has a limited length evidently shorter than those obtained from numerical integration. Therefore, the conclusions drawn from the analysis should be treated in place and interpreted correctly, especially taking into account the time scales.

Moreover, the future surveys will serve further measurements to the existing ones generating longer time series that are going to be the basis of more trusty dynamical analysis. The proposed scheme can be generalized to more than two planets in a system 555T. Kovács, MNRAS (2019) in preparation, thus, it can be used as a supportive method to the current efficient stability investigations.

Acknowledgements.

Special thank is devoted to Tamás Tél his valuable comments. The author is also grateful to the anonymous referee(s) for important remarks. This work was supported by the NKFIH Hungarian Grants K119993, PD121223. The financial support of Bolyai Research Fellowship is also acknowledged.