Analyzing the Stability of Non-coplanar Circumbinary Planets using Machine Learning
Zhihui Kong, Jonathan H. Jiang, Zong-Hong Zhu, Kristen A. Fahy, Remo Burn
1 Analyzing the Stability of Non-coplanar Circumbinary Planets using Machine Learning
Zhihui Kong , Jonathan H. Jiang , Zong-Hong Zhu , Kristen A. Fahy , Remo Burn Department of Astronomy, Beijing Normal University, Beijing, China 2.
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, USA 3.
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
Copyright: © 2020 All rights reserved.
Keywords: Planet, Binary Star, Machine Learning, Stability * Corresponding to: [email protected]
Abstract
Exoplanet detection in the past decade by efforts including NASA’s Kepler and TESS missions has discovered many worlds that differ substantially from planets in our own Solar system, including more than 400 exoplanets orbiting binary or multi-star systems. This not only broadens our understanding of the diversity of exoplanets, but also promotes our study of exoplanets in the complex binary and multi-star systems and provides motivation to explore their habitability. In this study, we analyze orbital stability of exoplanets in non-coplanar circumbinary systems using a numerical simulation method, with which a large number of circumbinary planet samples are generated in order to quantify the effects of various orbital parameters on orbital stability. We also train a machine learning model that can quickly determine the stability of the circumbinary planetary systems. Our results indicate that larger inclinations of the planet tend to increase the stability of its orbit, but change in the planet’s mass range between Earth and Jupiter has little effect on the stability of the system. In addition, we find that Deep Neural Networks (DNNs) have higher accuracy and precision than other machine learning algorithms. Introduction
Previous research suggested that sun-like stars are common in binary systems, where almost half of sun-like stars have a companion (Raghavan et al. 2010; Moe & Di Stefano 2017) and long-period ( ≤ 𝑑 ), giant ( 𝑅 𝑝 > 6 𝑅 ⊕ ) circumbinary planets appear to be as common as those orbiting single stars (Armstrong et al. 2014). However, since the highly successful transit method is not suitable for observing long-period planetary systems (Gongjie, Holman & Tao 2016; Kipping & Sandford 2016), only a small number of circumbinary planets have been detected so far. The Open Exoplanet Catalogue (OEC) and NASA Exoplanet Archive have listed 25 circumbinary planetary systems, where the planet essentially orbits the center of mass of the stars (P-type orbits). Current and future observation missions (e.g. TESS and JWST) and application of the new methods (e.g. Transit Timing Variations) have potential to discover more complex systems of exoplanets. Study on the stability of planetary orbits in the binary system can help us to understand the complex formation of planetary systems, and to find habitable planets in the binary systems. Assessing the stability is typically done either thorough N-body simulations or using the analytic formula developed from Rabl & Dvorak (1988) by Holman & Wiegert (1999), who regressed a second-order polynomial in their dynamical simulations. N-body simulation method, in principle, uses regressions of three-or-more body systems that check the resulting dynamical stability of each trial set of parameters, however, this requires high computation cost. Using second-order polynomial method, as noted in Holman & Wiegert’s work, orbital resonance makes the ‘Islands of Instability’, in which the binary or multiple stars’ gravitational influence on each other changes their orbital paths, creating a “region of instability” over time, i.e. if a circumbinary planet moves into this region, chances are it’s going to get kicked out. However, the second-order polynomial method can’t evaluate the stability of circumbinary planets. Lam & Kipping (2018) use the machine learning (ML) method to predict successfully the ‘Islands of Instability’ and provides a new way to study the stability of circumbinary planet. In this study, we improved the circumbinary planet model based on Lam & Kipping’s work. We focused on the stability of non-coplanar circumbinary planetary systems for planets with non-zero mass, which was not included in the previous work. Although, there are theoretical indications that planetary inclinations should be preferentially coplanar with the binary’s plane (Foucart & Lai 2013), the actual exoplanetary orbits may not lie in the same plane, with influences possible from protoplanetary disk alignment, planet-planet scattering (e.g. Chatterjee et al. 2008) and other sources of orbital evolution (e.g. Kley & Haghighipour 2014). There is a study (D.J. Armstrong et al. 2014) that suggested the rate of occurrence of circumbinary planets and planetary inclination have a strong connection and that the discovery rate will increase significantly for planets with non-zero inclination orbits. Therefore, studying the stability of non-coplanar circumbinary planetary systems is important. In Section 2, we describe the circumbinary planet model and how we generated a large suite of training data using numerical simulations. In Section 3, we describe the network architecture of the deep neural networks (DNNs). In Section 4, we show the effect of planetary inclination and mass on the stability of circumbinary planetary systems and the performance of our ML algorithm. Finally, we discuss possible applications of our method and future areas for development in Section 5. Model and Training Dataset 2.1
The model of circumbinary planet system
To study the P-type orbit stability, we use the N-body simulation package, REBOUND, which first developed by Rein & Liu (2012). We use IAS15 integrator which is a non-symplectic integrator with adaptive time-stepping (Rein & Spiegel 2015). In the REBOUND environment, we set up a circumbinary planetary system. The central binary consists of mass m A and m B with total mass m A + m B = 1 m ⊙ and the mass fraction of the binary is μ , where μ≡ m B / ( m A +m B ) (1) The binary orbit has semimajor axis, a b , and the orbital eccentricity is e b . The Keplerian orbit of the planet with mass m p around the center of mass of the binary is defined initially by two orbital elements, its semimajor axis a p , inclination i . We ran 10 REBOUND experiments with mass of the planet ( m p = m ⊕ or m p = m jupiter ), with each simulation lasting 10 binary period. For each simulation, we uniformly sampled the mass ratio, μ∈ [0, 0.5], the initial binary eccentricity, e b ∈ [0, 0.99], and the inclination of the planet, i ∈ [0, π ]. The initial eccentricity of the planet e b is set to zero. Table 1 illustrates the details. Table 1:
Initial conditions for circumbinary planetary system
Binary Parameters 𝑚 & + 𝑚 ( = 1 𝑚 ⊙ 𝑚 & > 𝑚 ( Planet Parameters <=> ≤ 𝑎 ? ≤ 6𝑎 <=> 𝑠𝑡𝑒𝑝 𝑤𝑖𝑑𝑡ℎ: 0.01 𝑚 ? = 1𝑚 ⊕ 𝑜𝑟 1𝑚 DE?=FGH
Inclination: 𝑖 ∈ [0, 2𝜋]
Step width:
OPQ
Sample size: R 𝑇 = 10 T 𝑇 U Stability criterion
The circumbinary planetary system can become unstable by the planet being either ejected from the system, or colliding with one of the stars. The former instance can be tested by evaluating whether the radial component of the planet’s velocity, v p , exceeds the escape velocity of the binary pair, v esc , indicating that the planet will escape to infinity. The latter can be tested using the orbital radius of planet. If the planet’s orbit passes the interior to that of the binary orbit then they will collide. We perform these two tests every quarter with a binary period and label the simulations as “instable” and terminate them. Next, there are two typical examples of the circumbinary systems which were labelled as ‘stable’ and ‘unstable’. For instance, there are some systems with the different inclination and the same initial parameters, as e b = 0.4, μ = 0.1, m p = m ⊕ and the semimajor axis of the planet is 0.93 × a HW99 , Holman & Wiegert(1999), which is given by 𝑎 𝐻𝑊99 = Y \ + Y \ 𝑒 +(−2.22 ± 0.11)𝑒 ‘ + (4.12 ± 0.09)𝜇 +(−4.27 ± 0.17)𝑒𝜇 + (−5.09 ± 0.11)𝜇 ‘ +(4.61 ± 0.36)𝑒 ‘ 𝜇 ‘ (2) Figure 1:
Simulated orbits for four different inclinations.
The red asterisk shows the binaries’ star and the blue dot is the planet; the light-gray region is the orbit trails of planet. On the top, there are unstable systems after a few binary periods; On the bottom, there are stable systems after 10 planetary periods. Figure 1, shows the systems after 10 planetary periods. The red asterisk shows the binary stars and the blue dot is the planet; the light-gray region is the orbit trails of planet. Figure 2 shows the speed of the planets in the systems changing over time. Through the above samples, we know that the inclination of planets will affect the stability of their orbit. For example, the top panels ( i= i= unique REBOUND simulations and recorded their label as ‘stable’ or ‘unstable’. In the next section, these samples were used as training data for ML. Figure 2 : Evolution of two systems with different initial planetary inclinations of 60 o and 120 o : The x-axis is the time-step and the y-axis is the planet’s velocity squared. The unit is the square of Third Cosmic Velocity, 𝑉 T , and the square of escape velocity is 0.59 𝑉 T in these systems. Machine learning models 3.1
Deep neural networks
In past decades, techniques of ML have been rapidly developing, impacting many areas in industry, such as self-driving cars, speech-recognition systems, manufacturing, energy harvesting, and more. At the same time, scientists have increasingly become interested in the potential of ML for fundamental research, such as physics. To some extent, both ML and physics share their methods and goals. Both ML and physics are concerned about the process of gathering and analyzing data to design models that can predict the behavior of complex systems (Giuseppe Carleo et al. 2019).
Deep neural networks (DNNs) have been demonstrated to be powerful tools in predictive applications, exoplanetary studies is no exception (e.g. Kipping & Lam 2017; Lam & Kipping 2018). In our work, we used Pytorch which used a syntax similar to Python NumPy to describe models.
Network architecture
In supervised learning, we are given a set of n samples of data; let us denote one such sample X j ∈ R p with j=1, … ,n . For each sample X j , we are further given a label y j ∈ R d , most commonly 𝑑 =1. The goal of supervised learning is finding a function f so that when a new sample X new is presented without its label, then the output of the function f ( X new ) approximates well the label (Giuseppe Carleo et al. 2019). In our work, we used the DNNs to predict the stability of non-coplanar circumbinary planets. In the previous section, we have prepared the training dataset. Then we take a set of N-dimensional ( μ, e, a_p, i ) inputs and pass the data through ‘hidden’ layers of neurons, in which non-linear activation functions transform the data. The output of the network is compared with the true label (stable or unstable), and the error calculated by some loss function. Then we used the back-propagation algorithm to adjust the weights. Finally, we can find the minimum value of the loss function through iterative operations. Neural network architecture can be determined by the following parameters: the number of hidden layers, the number of neurons in these hidden layers, the activation functions, the loss function, and dropout. In our work, the net’s architecture is set up by 4 hidden layers of 256 neurons each, using ReLU (rectifier linear unit) activation function, CrossEntropyLoss loss function, and the SGD optimizer. The ReLU activation function, which lives in the hidden layers, is defined as f ( x ) ≡ max(
0, x ) (3) The SGD optimizer reduces computational cost at each iteration, that is a good estimate of the gradient descent. Dropout is a technique to avoid overfitting, an issue endemic to models as complex as DNNs . Feature selection
In section 2, We ran 10 REBOUND experiments for each mass of the planet ( m p = m ⊕ or m p = m jupiter ), hence we have two sets of training dataset, and each group has four parameters and one label (‘stable’ or ‘unstable’). Then we select μ, e, a p , i as features to train our DNNs. Before starting, we need to preprocess the data. For the continuous value, like μ, e, a p , we need put them on a common scale. The standardization is defined as x ← ( x − u ) / σ (4) Where u is the average, and σ is stand deviation. For the discrete value, i , we replace them with a one-hot encoding. For example {dog, cat, chicken} can be written using one-hot encoding as {(1,0,0), (0,1,0), (0,0,1)}. Finally, we have 22 features to train our DNNs. Iterative learning
When training a DDNs, one generally aims to pass the dataset through the network iteratively until the weights of the neural net are satisfactorily tuned, usually to the minimum value of the loss function. To achieve this, we put the dataset into the training set and validation set, with 9:1, that we use to tune parameters, select features, and make other decisions regarding the learning algorithm. Next, we increased the number of training rounds, measured in ‘epoch’, until the loss on the validation set stopped significantly decreasing. An ‘epoch’ is simply a forward pass through the network and the subsequent back-propagation for all training examples. We set the ‘epoch’ to 100, and there are some other parameters as: batch size is 512, learning rate is 0.5. The next section, we will discuss the result of the REBOUND and DNNs. Results 4.1 The stability of non-coplanar circumbinary planetary systems
We have discussed a simple example of the stability of non-coplanar circumbinary systems in the section 2. Figure 2 illustrates that the inclination of planet would affect the stability of planetary systems. Then, we selected 10 samples from the m p = m ⊕ , and drew a figure to describe the distribution of variables, as in Figure 3. Figure 3 illustrates that the samples of encounter are mainly concentrated in the area with the smaller μ , a p , i value and the larger e value; the samples of escape are mainly concentrated in the area with the smaller a p value and the larger e value; the samples of stable are mainly concentrated in the area with the smaller e value and the larger a p , and i value. Next, we select the stable samples from all 10 samples with m p = m ⊕ and draw a heatmap at Figure 4 with the variables of a p and e . Figure 3:
There are 10 data in the figure with m p = m ⊕ . There are three kernel density estimate (KDE) subgraphs of μ, e, a p , and one histogram of i . There are two colors to distinguish the orbital status, seeing the legend. Figure 4:
There are stable samples with m p = m ⊕ . In the figure, different colors represent the sample density, and the black dotted line emphasizes that the “Islands of instability” may result in uneven density transition. The black dotted line of the Figure 4 marks the significant change in density of the region that is called ‘Islands of instability’ (Holman & Wiegert 1999). These islands are driven by mean motion resonances. Lam & Kipping (2018) use the ML method to successfully predict these islands but their circumbinary planetary systems are coplanar and ignore the mass of the planet. In our work, Figure 4 illustrates this phenomenon, ‘Islands of instability’, also exists in the non-coplanar circumbinary planetary systems with non-zero planetary mass and unlike Lam & Kipping’s the model the border is not very clear.
Figure 5:
This is a joint histogram using hexagonal bins of 606926 stable samples with 𝑚 ? = 𝑚 ⊕ and hexagon darker description of the greater the density of the samples. The unit of inclination is radian and the unit of 𝑎 ? is AU. We then continue to use these stable samples to draw a joint histogram using hexagonal bins with the variables of a p and i , Figure 5. Figure 5 illustrates that the higher a p value, the more stable orbits of the planets, because planet far away from center of the binary is subject to smaller gravitational perturbations by the binary stars. In addition, we have found the smaller the angle between the planetary plane and binary plane, the smaller the stability boundary. Especially when the planet in retrograde orbit appears more stable (the upper-left corner of the Figure 5). Figure 6:
There are R data in the figure with 𝑚 ? = 𝑚 DE?=FGH . There are three kernel density estimate (KDE) subgraphs of 𝜇, 𝑒 𝑎 ? , and one histogram of 𝑖 . Next, let’s look at these samples of m p = m jupiter , Figure 6 shows that the distribution of these variables look very similar to Figure 3. This suggests that, in our model, the planet’s mass range is between that of Earth and Jupiter, with little effect on the stability of the system. In the section 4.1, we find the orbital resonance will make the ‘Islands of Instability’, which can’t be fitted by second-order polynomial (Holman & Wiegert 1999). Therefore, we train ML algorithms to predict the stability of non-coplanar circumbinary planetary systems with m p = m ⊕ and m p = m jupiter . Normally, the performance of ML is related to the size of trainset data, the more data the higher accuracy but using N-body simulations to generate the trainset data is computationally expensive. Hence, we choose the samples of 10 and 10 to train the Deep Neural Networks (DNNs). Once training is complete, we can use DNN to quickly predict more samples, that is an important advantage of ML algorithms. Next, we split the dataset 1⁄10 into test and training sets. The trainset data is used to train the structure of the neural network, so we save the model parameters of DNNs, and then using the test-set test the performance of DNNs. In this work, we used the F1 score to evaluate. Here is the definition: F1= ( ) / ( P+R ) (5) Where P is the Precision,
P=TP /( TP+FP ); and the R is the Recall, R=TP /( TP+FN ). We tallied the
TP, FN, FN, TN in the Table 2. We have also calculated the Accuracy, acc = ( TP+TN ) / (
TP+FP+FN+TN ). Table 2
Truth Prediction stable unstable stable TP FN unstable FP TN
In addition to our own design of the DNNs model, we also compared the performance of several other ML algorithms in our sample, Table 3 shows the results of different ML algorithms. The result, in the Table 3 illustrates that DNNs can well predict the stability non-coplanar circumbinary planetary systems with planets of non-zero mass.
Table 3:
The performance of different machine learning algorithms
Model Accuracy Precision Recall F1 score
DNNs
CatBoost Classifier
Extreme Gradient Boosting
Random Forest Classifier Summary and Discussion
With the development of transit timing variations (TTVs), more and more multi-planetary systems have been found. This has greatly expanded our understanding of the diversity of exoplanets, and has led to the scientific study of complex planetary systems. Recently, many studies have focused on non-coplanar planetary systems. (Doolin & Blundel 2011; Chen et. al 2020). Specially, Chen’s work used the model similar to ours, that show non-coplanar circumbinary planetary orbits have nodal precession. Contrarily, we used a more accurate integrator IAS15. In this way, we can study the transition regions of orbital stability and instability. Our results show that there are indeed ‘Islands of instability’, possibly caused by mean motion resonances (MMRs), in the non-coplanar circumbinary planetary systems with the non-zero plantary mass. This is significant to study the MMRs to form planets in circumbinary formations. In addition, the results in the section 4.1 also show that planets in the retrograde orbit appear more stable, this is significant in understanding the migration of planetary orbits. The excitation of planets to retrograde orbits can occur due to (a) interactions with other planets in the system or (b) due to interactions with external stars. The first pathway (a) for planet formation around single stars can be explored using planetary population synthesis models (Benz et al. 2014, Mordasini et al. 2018). The recent simulations presented by Emsenhuber et. al 2020 integrate N-body interactions of up to 100 growing protoplanets and results in the formation of 35 retrograde planets in 1000 systems after 20 Myr of integration. The fraction of systems with retrograde planets in their simulations is only 1.5% and the typical retrograde planet has a mass comparable to mars. The second production method (b) of retrograde planets is more promising to lead to observable results: Stock et al. 2020 model the stage of planet - star interactions using realistic conditions of stellar clusters for planetary systems similar to the Solar system. They find that encounters with stars lead to the excitation of retrograde orbits in approximately a percent of Solar system giant planet analogues. These more massive planets are more easily observable. We would like to point out that the stars in the stellar cluster of Stock et al. 2020 do not include primordial binary stars. Therefore, binary systems are too rare in their simulations and this shortcoming should be improved on in the future. Another important aspect is that we used the ML method to predict the stability of non-coplanar circumbinary planetary systems. Previous research has shown that ML are well-suited for predicting the stability of multi-planetary systems (Tamayo et. al 2016; Kipping & Lam 2017; Lam & Kipping 2018). ML has the following advantages: 1. ML is very fast compared to N-body simulation. In our work, we need a dozen days to generate the training data by N-body simulation, and training DNNs model, however, takes about 30 minutes with laptop computing power. In addition, it takes a just few minutes to use the trained DNNs to predict more samples. 2. ML can describe the transition regions of orbital stability and instability better than polynomial function presented in Holman & Wiegert (1999). In the section 4.2, we show the performance of DNNs, that can predict the orbital stability with great F1 score (98.25%). In the future, we can study the process of planet formation with non-coplanar orbits, or the long-term stability of non-coplanar circumbinary planetary systems with MMRs. In fact, these works have already begun on multi-planet systems with one center star (Sotiris Sotiriadis et. al 2018), and even a study which also used ML (Daniel Tamayo et. al 2020). Acknowledgement
Authors JHJ and KAF were supported by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. We also acknowledge NASA ROSES XPR program for support. Authors ZK and ZZ thanks the support by the National Natural Science Foundation of China under Research Grants No. 11633001, 11920101003 and 12021003, and the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000. We thank Xiaoming Jiang for helpful comments and discussions during development of the this study.
Data availability:
The simulated data underlying this article are described in the article. For additional questions regarding the data and code sharing, please contact the corresponding author at [email protected].
References
Armstrong D. J., Osborn H., P., Brown D., J., A. et al. 2014, MNRAS, 444, 1873, doi: 10.1093/mnras/stu1570 Benz W, Ida S, Alibert Y. et al. 2014, doi: 10.2458/azu_uapress_9780816531240-ch030 Chatterjee Sourav, Ford Eric B., Matsumura Soko, et al. 2008, ApJ, 686, 580, doi: arXiv:astro-ph/0703166 Cheng Chen, Alessia Franchini, Lubow Stephen H. et al. 2019, MNRAS, 490, 5634, doi: 10.1093/mnras/stz2948 Deepak Raghavan, Harold A McAlister, Todd J Henry, et al. 2010, ApJS, 190, 1, doi: 10.1088/0067-0049/190/1/1 Emsenhuber Alexandre, Mordasini Christoph, Burn Remo. Et al. 2020, A&A, 32, 21, arXiv:2007.05562 Foucart Francois, Lai Dong. 2013, ApJ, 764, 106, doi: 10.1088/0004-637X/764/1/106 Gabl G., Dvorak R. 1988, A&A, 191, 385. Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer et al. 2019, RvMP, 91, 045002, doi: 10.1103/RevModPhys.91.045002 Gongjie Li, Matthew J. Holman, Molei Tao. 2016, ApJ, 831, 96, doi: 10.3847/0004-637X/831/1/96 Holman Matthew J., Wiegert, Paul A. 1999, AJ, 117, 621, doi: 10.1086/300695 Kipping David M., Sandford, Emily. 2016, MNRAS, 463, 1323, doi: 10.1093/mnras/stw1926 Kipping David M., Lam Christopher. 2017, MNRAS, 465, 3495, doi: 10.1093/mnras/stw2974 Kley Wilhelm, Haghighipour Nader. 2014, A&A, 564, A72, doi: 10.1051/0004-6361/20132323513