Exoplanets prediction in multiplanetary systems

Abstract

We present the results of a search for additional exoplanets in all multiplanetary systems discovered to date, employing a logarithmic spacing between planets in our Solar System known as the Titius–Bode (TB) relation. We use the Markov Chain Monte Carlo method and separately analyse 229 multiplanetary systems that house at least three or more confirmed planets. We find that the planets in $\sim 53\%$ of these systems adhere to a logarithmic spacing relation remarkably better than the Solar System planets. Using the TB relation, we predict the presence of 426 additional exoplanets in 229 multiplanetary systems, of which 197 candidates are discovered by interpolation and 229 by extrapolation. Altogether, 47 predicted planets are located within the habitable zone of their host stars, and 5 of the 47 planets have a maximum mass limit of 0.1–2 ${\rm M}_{\oplus}$ and a maximum radius lower than 1.25 ${\rm R}_{\oplus}$ . Our results and prediction of additional planets agree with previous studies’ predictions; however, we improve the uncertainties in the orbital period measurement for the predicted planets significantly.

1. Introduction

The number of detected exoplanets is growing rapidly, so that over 4000 exoplanets have been detected and confirmed to date. There are also thousands of other candidate exoplanets that require further follow-up observation. The Kepler spacecraft plays a key role in detecting these systems. The Kepler space mission’s main goal was to discover Earth-size exoplanets in or near the habitable zone (HZ) of solar-like stars and determine the fraction of the hundreds of billions of stars in our galaxy that might possess such planets (Koch et al. 1998).

Photometry with the transit method is the most successful exoplanet discovery method, which has been used by Kepler space mission, Transiting Exoplanet Survey Satellite (TESS), and many ground-based observatories (Borucki et al. 2010; Ricker et al. 2015; Deeg & Alonso 2018). However, this method has its own difficulties. For example, transits are detectable only when the planet’s orbit happens to be almost exactly aligned with the observer’s line of sight. This covers only a small fraction of exoplanets. Furthermore, the planet’s transit lasts for a small fraction of its total orbital period. As a result, it is not very likely to detect planets’ transits, especially those with long orbital periods. This study sets out to predict the existence of additional undetected planets in multiple exoplanet systems.

In our Solar System, there is a simple logarithmic spacing between planets, which has been known for over two centuries as the Titius–Bode (TB) law. Its classical relation is

(1) \begin{equation}a_{n}=0.4+0.3\times2^{n},\end{equation}

where $a_{n}$ represents the semi-major axis of the nth planet in AU. The planet Mercury corresponds to $n=-\infty$ , Venus to $n=0$ , Earth to $n=1$ , and so on (Nieto & Davis 1974). After the discovery of the planet Uranus in 1781 by Frederick William Herschel, it was recognised that the TB law predicted this planet’s semi-major axis (Nieto & Davis 1974).

The discovery of the TB relation motivated many observation programmes to investigate and detect the lost fifth planet, which eventually led to the discovery of the asteroid Ceres (Sawyer Hogg 1948). The predictions made using the TB relation also played a key role in exploring the planet Neptune, but not as accurately as Uranus. Interestingly, the satellite systems of the giant planets also follow a TB relation (Lyttleton 1960; Brookes 1970).

The TB relation was used effectively to predict lost undetected objects in our Solar System. It was believed that this relation could help make similar predictions in detected multiple exoplanet systems, too. The five-planet 55 Cnc system was one of the first multiple exoplanet systems where Poveda & Lara (2008) applied TB to predict the undetected planets. They found that a simple exponential TB relation reproduces the semi-major axes of the five observed planets. They also predicted two additional planets at distances of 2 and 15 AU. Using the 55 Cnc system, Chang (2008) also checked whether the TB relation is enforceable on exoplanetary systems by statistically analysing the distribution of the ratio of periods of two planets in the 55 Cnc system, by comparing it with that derived from the TB relation. Chang (2010) again repeated this calculation for 31 multiple exoplanet systems and concluded that the adherence of the Solar System’s planets to the TB relation might not be fortuitous; thus, we could not ignore the possibility of using the TB relationship in exoplanetary systems.

Moreover, Lara, Poveda, & Allen (2012) showed that like 55 Cnc, 10 other planetary systems (ups And, GJ 876, HD 160691, GJ 581, Kepler-223, HR 8799, Kepler-20, Kepler-33, HD 10180, and Kepler-11) host four or more planets that also obey a similar (but not identical) TB relation.

Bovaird & Lineweaver (2013) (hereafter, BL13) used a sample of 68 multiple exoplanet systems with at least 4 planets, including samples of both confirmed and candidate systems. They identified a sample of exoplanet systems that are likely to be more complete and tested their adherence to the TB relation. They found that most of these exoplanetary systems adhere to the TB relation better than the Solar System. Using a generalised TB relation, they predicted 141 additional exoplanets, including a planet with a low radius ( $R<2.5\,{\rm R}_{\oplus}$ ) and within the HZ of the Kepler-235.

Using the predictions made by BL13, Huang & Bakos (2014) analysed Kepler’s long-cadence data to search for 97 of those predicted planets and obtained a detection rate of $\sim 5\%$ . Considering the possibility that the remaining predicted planets were not detected because of their small size or non-coplanarity, the detection rate was estimated to be less than the lower limit of the expected number of detections. They concluded that applying the TB relation to exoplanetary systems and using its predictive power in BL13 could be questionable.

Bovaird, Lineweaver, & Jacobsen (2015) (hereafter, BL15) used the BL13 method to predict the existence of additional planets in 151 Kepler multiple exoplanet systems that contain at least 3 transiting planets. They found 228 undetected exoplanets and, on average, $2\pm1$ planets in the HZ of each star. They found that apart from the five planets detected by Huang & Bakos (2014), one additional planet in Kepler-271 had been detected around its predicted period. The completeness of observational data was taken into account as an important issue in their analysis. They declared a completeness factor due to the intrinsic noise of the host stars, planet radius, highly inclined orbits, and the detection techniques. Thus, they did not expect all predicted planets to be detected. They estimated the geometric probability to transit for all 228 predicted planets and highlighted a list of 77 planets with high transit probability, resulting in an expected detection rate of $\sim 15\%$ , which was about three times higher than the detection rate measured BL13.

Recently, Lara, Cordero-Tercero, & Allen (2020) used data from 27 exoplanetary systems with at least 5 planets and applied their proposed method to find the reliability of the TB relation and its predictive capability to search for planets. They removed planets from the system one by one and used TB relation to recover them, where they were able to recover the missing one 78% of the time. This number was much higher than when they tried this with random planetary systems, where 26% of planets were recovered. Using statistical tests, they showed that the planetary orbital periods in exoplanetary systems were not consistent with a random distribution and concluded it to be an outcome of the interactions between true planets.

The purpose of this study is to test the adherence of 229 multiple exoplanet systems with at least 3 detected planets to the TB relation and compare their adherence rate with the Solar System’s using Markov chain Monte Carlo (MCMC) (see Goodman & Weare 2010; Foreman-Mackey et al. 2019) as a precise method for regression while considering the reported uncertainties of physical values. We also aim to predict the existence of additional undetected exoplanets and estimate their physical properties, either maximum mass or maximum radius. Moreover, we highlight exoplanets located within the HZ of their host stars. Using a sample of seven multiple exoplanet systems with detected planets after predictions made by BL15, we also aim to determine if the TB relation is a reliable method of identifying undetected member planets. Finally, we compare our predictions with those from BL15.

This paper is organised as follows: Section 2 describes the data and methods used to predict exoplanets and estimate their maximum mass and maximum radius. Section 3 presents our results. Section 4 summarises our main results and conclusions.

Figure 1. The mass–radius distribution of the exoplanets in our sample, separated into five groups based on their detection methods: imaging (black stars), radial velocity (red circles), transit (blue squares), and transit timing variations (green triangles).

2. Data and Method

2.1. Data

We used physical parameters of exoplanets, including planetary orbital period, radius, and mass, and the stellar mass and radius from two exoplanet databases: the NASA Exoplanet Archivea and the Extrasolar Planets Encyclopedia.b We also used conservative and optimistic limits of the HZ of available stars from the Habitable Zone Gallery.c

The total number of multiple exoplanet systems with at least 3 confirmed planets available is 230 systems to date, hosting 818 planets. And 81.5% of these planets have been detected using the transit method and 16.4% using the radial velocity method. The remaining 2.1% of planets have been identified using transit timing variations, imaging, pulsar timing, and orbital brightness modulation methods. Figure 1 represents the mass–radius distribution of exoplanets for five groups of planets, separated according to their detection methods. Of the 230 systems, 142 systems host 3 planets, 59 systems host 4 planets, 20 systems include 5 planets, and 7 systems host 6 planets. TRAPPIST-1 is the only system that contains seven planets, and KOI-351 the only one that contains eight. Figure 2 represents the distribution of the orbital period of member exoplanets in various exoplanetary systems used in this study.

Figure 2. The distribution of the orbital period of 813 member exoplanets in 229 multiple exoplanet systems hosting at least 3 planets (grey bars), 4 planets (red bars), and 5 (or more) planets (blue bars).

The exoplanetary systems 55 Cnc, GJ 676 A, HD 125612, K2-136, Kepler-132, Kepler-296, Kepler-47, Kepler-68, and ups And consist of binary stars. GJ 667 C and Kepler-444 are also triple-star systems. The system Kepler-132 possesses four planets such that two of them (Kepler-132 b and Kepler-132 c) have roughly similar orbital periods of 6.1782 and 6.4149 d. After more detailed studies, it was found that Kepler-132 b and Kepler-132 c cannot orbit the same star (Lissauer et al. 2014). Consequently, we exclude Kepler-132 from our analyses; our final sample of exoplanetary systems thus contains 229 systems.

2.2. Methods

2.2.1. The TB relation and the prediction of additional exoplanets in multiple exoplanet systems

We use the TB relation to predict additional undetected planets for all multiple exoplanet systems, with at least three confirmed planets.

The TB relation can be written in terms of the orbital periods as follows:

(2) \begin{equation}P_{n}=P\alpha^{n}, n=0,1,2,...,N-1,\end{equation}

where $P_{n}$ is the orbital period of the $n{\rm th}$ planet, and P and $\alpha$ are fitting parameters. In analysing the exoplanetary systems, the Solar System is used to guide how well each system adheres to the TB relation. In logarithmic space, the TB relation is written as follows:

(3) \begin{equation} \log P_{n}=\log P + n\log\alpha=b+m\times n, n=0,1,2,...,N-1,\end{equation}

where $b=\log P$ and m= $\log\alpha$ are intercept and slope of the relation, respectively. For a system with N planets, the $\chi^2/dof$ value can be calculated by:

(4) \begin{equation}\frac{\chi^2(b,m)}{N-2}=\frac{1}{N-2}\sum_{n=0}^{N-1}\left[\frac{(b+m\times n)-\log P_{n}}{\sigma}\right]^2,\end{equation}

$N-2$ is the number of degrees of freedom (N planets—2 fitted parameters (b and m)), and $\sigma$ represents the system’s sparseness or compactness. The sparseness/compactness of a system is calculated from $\sigma$ = 0.273 $S_p$ , where $S_p$ represents the average log period spacing between planets as defined by:

(5) \begin{equation}S_{p}=\frac{\log P_{N-1}-\log P_{0}}{N},\end{equation}

where N is the number of planets in the system, and $P_{N-1}$ and $P_{0}$ are the largest and smallest orbital periods in the system, respectively. We plot the $\sigma$ values (see Equation (4)) as a function of sparseness/compactness ( $S_{p}$ ; see Equation (5)) in Figure 3. The black line goes through two points: the origin and the specific $S_{p}$ and $\sigma$ values for the Solar System ( $\sigma$ is the value required for Solar System to yield $\chi^2/dof=1$ in Equation (4)). The black points show the exoplanet systems with no planet insertions, grey points indicate systems before planet insertion, and red points indicate systems after insertions have been made.

Figure 3. $\sigma$ as a function of the average log period spacing between planets, $S_p$ (Equation (5)), of exoplanet systems. The black line goes through two points: the origin (0,0) and ( $S_{p}$ , $\sigma$ ), where $S_{p}$ is the compactness of the Solar System and $\sigma$ is the value required for the Solar System to yield $\chi^2/dof=1$ in Equation (4). The cyan triangle shows the Solar System, and the black dots show the exoplanet systems with no planet insertions (systems with $\chi^2/dof\leq1$ ). Grey dots indicate systems before planet insertions, and red dots indicate the ( $S_{p}$ , $\sigma$ ) of the systems after insertions have been made.

Using a proper value for $\sigma$ , $\chi^2/dof=1$ is adjusted for the Solar System case. If the detected planets in a system adhere to the TB relation better than the planets of the Solar System ( $\chi^2/dof\leq1$ ), we only predict an extrapolated planet beyond the outermost detected planet. If the detected planets adhere worse ( $\chi^2/dof>1$ ), we begin the interpolation process; for the first step of interpolation, one new specific planet, from atleast 5 000 hypothetical planets that have a random period between the innermost and outermost detected planets, is inserted into the system covering all possible locations between the two adjacent planets, and a new $\chi^2/dof$ value calculated for each possibility. The new specific inserted planet is chosen from 5 000 cases when producing the minimum value of $\chi^2/dof$ . Similarly, inserting up to 10 new specific planets for each system step by step (two planets for the second step, three planets for the third step, etc.) covers all possible locations and combinations between 2 adjacent planets. The period uncertainty of an inserted planet ( $e_{\rm ins}$ ) is calculated using the uncertainties of detected planets in the same system (e.g., $e_{1}$ , $e_{2}$ , $e_{3}$ ,…) as follows:

(6) \begin{equation}e_{\rm ins}=\sqrt{e_{1}^{2}+e_{2}^{2}+e_{2}^{3}+...}.\end{equation}

We adopt the highest value of the parameter $\gamma$ , which is the improvement in the $\chi^2/dof$ per inserted planet, for identifying the best combination of detected and predicted planets. $\gamma$ is defined by:

(7) \begin{equation}\gamma=\frac{\left(\frac{\chi_{i}^2-\chi_{j}^2}{\chi_{j}^2}\right)}{n_{\rm ins}},\end{equation}

where $\chi_{i}^2$ and $\chi_{j}^2$ are the $\chi^2$ values before and after inserting of $n_{\rm ins}$ planets.

We analyse each system’s data separately and use the MCMC method to quantify the uncertainties of the best-fit parameters.

By applying the lowest signal-to-noise ratio (SNR) of the detected planets in the same system to the predicted planet’s orbital period, the maximum mass or maximum radius of the predicted planets is calculated.

For transiting detected planets, the maximum radius is calculated by:

(8) \begin{equation}R_{\max}=R_{\min SNR}\left(\frac{P_{\rm predicted}}{P_{\min SNR}}\right)^{0.25},\end{equation}

and for radial velocity detected planets, the maximum mass is calculated by:

(9) \begin{equation}M_{\max}=M_{\min SNR}\left(\frac{P_{\rm predicted}}{P_{\min SNR}}\right)^{\frac{7}{6}},\end{equation}

where $R_{\min SNR}$ , $M_{\min SNR}$ , and $P_{\min SNR}$ are the radius, mass, and orbital period of the detected planet with the lowest SNR, respectively. After calculating the maximum radius or maximum mass of the predicted planets, using the mass–radius relationship established by Bashi et al. (2017) (hereafter, B17) ( $R_{p} \propto M^{0.55}$ and $R_{p} \propto M^{0.01}$ for the small and large planets, respectively), we convert the radius values to mass values, and vice versa.

2.2.2. Dynamical stability and the transit probability of the predicted planets

We use the dynamical spacing criterion ( $\Delta$ ) to analyse how our predicted objects could be stable in their positions in the exoplanetary system. The dynamical spacing $\Delta$ between two adjacent planets with masses $M_{1}$ and $M_{2}$ and orbital periods $P_{1}$ and $P_{2}$ orbiting a host star with a mass of $M_{*}$ was defined by Gladman (1993) and Chambers, Wetherill, & Boss (1996) as follows:

(10) \begin{equation}\Delta=\frac{2M_{*}^{1/3}\left(P_{2}^{2/3}-P_{1}^{2/3}\right)}{\left(M_{1}+M_{2}\right)^{1/3}\left(P_{2}^{2/3}+P_{1}^{2/3}\right)}.\end{equation}

In a planetary system, two adjacent planets are less likely to be stable in their positions when the value of their dynamical spacing is small ( $\Delta \leq 10$ ). However, if they are stable, they are more likely to be in orbital resonance with each other. Using this criterion, we analyse 45 systems in our sample, which have predicted planets within their HZs. We calculate the $\Delta$ value for all pairs in these systems and investigate whether they are in resonance with each other or not (for planets whose mass values have not been reported in catalogues, we use the mass–radius relation of B17 to calculate their masses). To analyse the possibility of resonance, we calculate the period ratios of all adjacent planets considering an arbitrary threshold; the planet pairs that have $x\leq2\%$ are considered to be pairs in orbital resonance, where

(11) \begin{equation}x=\frac{\mid \frac{N_{j}}{N_{i}}-\frac{P_{n+1}}{P_{n}}\mid}{\frac{N_{j}}{N_{i}}},\end{equation}

$N_{i}$ and $N_{j}$ are positive integers with $N_{i}<N_{j}\leq5$ , and $P_{n}$ and $P_{n+1}$ are the orbital periods of two adjacent planets.

The transit phenomenon of a planet is only observable if the planetary orbit plane is close to the line of sight between the observer and the host star. In other words, the pole of the planetary orbit must be within the angle $d_{s}/a_{p}$ , where $d_{s}$ is the stellar diameter and $a_{p}$ is the planetary orbital radius. Then, the geometric transit probability ( $P_{\text{tr}}$ ) of a planet can be estimated using $d_{s}/2a_{p}$ , where $P_{\text{tr}}=0.5\%$ for an Earth-size planet at 1 AU orbiting a Solar-size star (Borucki & Summers 1984; Koch & Borucki 1996). According to the models of planetary systems, multiplanetary systems, like the Solar System, are assumed to be formed out of common protoplanetary discs (Winn & Fabrycky 2015). Therefore, the orbital planes should have small relative inclinations, so that the Kepler multiplanetary systems are highly coplanar BL15. We use this criterion to estimate the transit probability of predicted exoplanets and prioritise how soon the predicted planets can be detected.

3. Results

3.1. The prediction of additional exoplanets in exoplanetary systems

To investigate the probability of the existence of additional planets in exoplanetary systems, we apply the TB relation and the MCMC method to analyse the data of a sample of 229 systems that contain at least 3 confirmed exoplanets. We find that 122 systems adhere to the TB relation better than the Solar System without any need for interpolation. For those 107 systems that adhere to the TB relation worse than the Solar System, we insert up to 10 new additional planets. For example, Figure 4 shows the best linear regression (TB relation) for up to 10 additional planet inserts in the GJ 667 C system. For GJ 667 C, the highest $\gamma$ value is in the fourth step, where the number of inserted planets ( $n_{\rm ins}$ ) is four. Black dots represent the detected planets and red dots the predicted planets. The black line shows the best (mean) scaling relation, and the two dashed lines show the $\pm1\sigma$ uncertainties around this relation. The blue lines are a set of 100 different realisations, drawn from the multivariate Gaussian distribution of the parameters, where for the highest $\gamma$ value: $m=0.1924, b=0.856$ , and ln< $\sigma>=-3.73$ . The scatter covariance matrix is also estimated from the MCMC chain. Figure 5 illustrates the one- and two-dimensional marginalised posterior distributions of the scaling relation parameters for the fourth step of the linear regression for the GJ 667 C system.

Figure 4. The TB relation and steps of linear regressions are applied to the data of system GJ 667 C. Detected and predicted planets are shown with black and red dots, respectively. For each step, the value of $n_{\rm ins}$ represents the number of inserted planets. The values of $\gamma$ and $\chi^2/dof$ are also shown where the highest value of $\gamma$ is in the fourth step as the best combination of detected and predicted planets. The two black dashed lines show $\pm1\sigma$ uncertainties around the best scaling relation (black solid line). The blue lines are a set of 100 different realisations, drawn from the multivariate Gaussian distribution of the parameters (for the fourth step: $m=0.1924, b=0.856$ and ln< $\sigma>=-3.73$ ), and the scatter co-variance matrix is estimated from the MCMC chain.

Figure 5. The one- and two-dimensional marginalised posterior distributions of the scaling relation parameters for the highest $\gamma$ value corresponding to the fourth step (( $n_{\rm ins}=4$ )) of the linear regression of the GJ 667 C system, as shown in Figure 4.

Table 1 represents the data corresponding to Figure 4 and lists the number of inserted planets ( $n_{ins}$ ) (column 1), the values of parameters $\chi^2/dof$ (column 2), $\gamma$ for each step of interpolation (column 3), orbital periods (column 4), and orbital numbers (ON) of inserted planets (column 5).

Table 1. Data corresponding to Figure 4

^a Number of the inserted planet.

^b $^{\rm b}\gamma =(\chi_{i}^2 - \chi_{j}^2) (\chi_{j}^{2}\times n_{ins})$ , where $\chi_{i}^2$ and $\chi_{j}^2$ are the $\chi^2$ values before and after inserting of $n_{ins}$ planets, respectively.

^c The orbital number of the inserted planet.

After interpolating all 107 systems, we find that these systems adhere to the TB relation better than the Solar System or to approximately the same extent. Of these 107 interpolated systems, 50 systems need 1, 33 systems need 2, and the remaining 24 systems need more than 2 additional planets to be inserted. We also predict an extrapolated planet beyond the outermost detected planet for all systems in our sample. We predict the existence of 426 possible additional exoplanets in these systems, of which 197 are predicted by interpolation. It should be noted that six of the predicted planets in Kepler-1388, Kepler-1542, Kepler-164, Kepler-374, Kepler-402, and Kepler-403 have been flagged as ‘Planetary Candidates’ in the NASA Exoplanet Archive.

To verify whether some of the predicted exoplanets are in the HZ of their parent stars, we use conservative and optimistic definitions of the HZ (Kane & Gelino 2012). Among the predicted exoplanets, 47 exoplanets lie within the HZ of their host stars. Twenty-seven exoplanets out of 47 are located within the conservative HZ. Furthermore, 14 exoplanets in the HZ have been predicted by interpolation, and the remaining 33 exoplanets have been predicted by extrapolation. Table 2 presents the 47 predicted exoplanets that lie within the HZ of the host star. Columns 1, 2, and 3 present the id, the host star name, and discovery method (Dis.), respectively. Columns 4, 5, and 6 present the orbital period in days, the distance from the parent star (a) in AU, and the orbital number (ON). The estimated maximum radius and maximum mass in the Earth unit are presented in columns 7 and 8. Column 9 lists the transit probability ( $P_{\text{tr}}$ ). The conservative and optimistic HZ limits for AU are presented in columns 10 and 11, respectively. Kepler-167 is a four-planet system where we predict three interpolated additional exoplanets, in which two planets have an orbital period of 157.0 and 373.5 d, located within the conservative HZ. The host star Kepler-186 is also a five-planet system including two predicted exoplanets with orbital periods of 41.2 and 73.6 d; these two additional planets were found by interpolation and are located within the optimistic and conservative HZ, respectively.

Borucki et al. (2011) classified exoplanets into the following class sizes: Earth-size ( $R_{p}<1.25\,{\rm R}_{\oplus}$ ), super-Earth-size ( $1.25\,{\rm R}_{\oplus}\leq R_{p}<2\,{\rm R}_{\oplus}$ ), Neptune-size ( $2\,{\rm R}_{\oplus}\leq R_{p}<6\,{\rm R}_{\oplus}$ ), and Jupiter-size ( $6\,{\rm R}_{\oplus}\leq R_{p}<15\,{\rm R}_{\oplus}$ ). Following this classification, 5 of our predicted exoplanets within HZ have maximum radii within the Earth-size range, 11 super-Earth-size, 27 Neptune-size, 3 Jupiter-size, and 1 with a maximum radius larger than twice that of Jupiter’s. Using the proposed categories based on planet mass by Stevens & Gaudi (2013), our five predicted planets have maximum masses within the range of Earth ( $0.1\,{\rm M}_{\oplus}-2\,{\rm M}_{\oplus}$ ). In addition, 22 and 19 predicted planets have maximum masses within the range of super-Earth ( $2\,{\rm M}_{\oplus}-10\,{\rm M}_{\oplus}$ ) and Neptune ( $10\,{\rm M}_{\oplus}-100_{\oplus}$ ), respectively. As a result, among our 47 predicted exoplanets within HZ, there are only five exoplanets whose estimated maximum mass and radius are within the mass and radius range of Earth: the fourth and fifth planets of Kepler-186, the second planets of GJ 3138 and Wolf 1061, and the extrapolated planet of YZ Cet.

We use the dynamical spacing criterion ( $\Delta$ ) to investigate our predicted objects’ stability at HZ. We calculate the $\Delta$ values for all adjacent planet pairs in these 45 systems, which host 47 predicted exoplanets within their HZ. We find that when inserting our predicted planets into systems, the average percentage of pairs with $\Delta\leq10$ increases from $\sim 25\%$ to $\sim 38\%$ . Nevertheless, this alone cannot be a reason for instability. In our Solar System, there are pairs of objects whose dynamical spacing values are small ( $\Delta\leq10$ ); however, they are stable in their positions. Neptune and Pluto are one of those pairs; the dynamical spacing value for them is $\Delta\sim7.4$ . We know that Neptune and Pluto are stable because they are in a 3:2 orbital resonance with each other. The same could be true for exoplanetary systems as found by Lithwick & Wu (2012), Batygin & Morbidelli (2013), and Fabrycky et al. (2014). We calculate the ratio of the orbital periods of planet pairs in these 45 systems and find that the number of orbital resonances increases when our predicted planets are considered. Figure 6 shows the number of resonance pairs with and without considering our predicted planets as shown with the solid blue line and dotted black line, respectively. The values for the Solar System are also shown with reference to the orange dashed line. As shown in this figure, the $\Delta\leq10$ pairs are dominated by the pairs in resonance; hence, as a sample, the planets (both detected and predicted) of these 45 systems are more likely to be stable in their positions.

Table 3 lists our predicted exoplanets by extrapolation and the best-fit TB relations. Columns 1 and 2 present the host star name and the system discovery method (Dis.), respectively. Column 3 reports a flag that defines whether the system has already been analysed by BL15 (or BL13) (Y) or not (N). Columns 4–7 present $\chi^2/dof$ , slope (m), intercept (b), and predicted orbital period, respectively. Column 8 reports whether the predicted period values in this paper and BL15 (or BL13) are consistent within error (Y) or not (N). Columns 9–12, respectively, list the orbital number (ON), estimated maximum radius ( $R_{\text{Max}}$ ), and maximum mass ( $M_{\text{Max}}$ ) in the Earth radius and mass unit, and the transit probability ( $P_{\text{tr}}$ ). The columns have been sorted based on the transit probability in descending order.

Table 2. Predicted exoplanets within the HZ of host stars in multiplanetary systems. Columns 1, 2, and 3 present the id, host star name, and discovery method (Dis.). Columns 4, 5, and 6 present the orbital period in days, distance from the parent star in AU, and the orbital number (ON). The estimated maximum radius ( $R_{\text{Max}}$ ) and maximum mass ( $M_{\text{Max}}$ ) in the Earth unit are presented in columns 7 and 8. Column 9 lists the transit probability ( $P_{\text{tr}}$ ). The conservative ( $HZ_{\rm Cons}$ ) and optimistic ( $HZ_{\rm Opt}$ ) HZ limits in AU are presented in columns 10 and 11, respectively

^a Discovery method of the system: ‘Tr’ and ‘RV’ represent transit and radial velocity, respectively.

^b Orbital numbers (ON) followed by ‘E’ indicate the extrapolated planets, and followed by ‘C’ indicate that the corresponding orbital periods have been flagged as “Planetary Candidate” in the NASA Exoplanet Archive.

Table 3. Systems with only extrapolated planet predictions. Columns 1 and 2 present the host star name and discovery method (Dis.). Column 3 reports a flag that defines whether the system has already been analysed by BL15 (or BL13) (Y) or not (N). Columns 4–7 present $\chi^2/dof$ , slope (m), intercept (b), and predicted orbital period, respectively. Column 8 reports whether the predicted period values in this paper and BL15 (or BL13) are consistent within error (Y) or not (N). Columns 9–12, respectively, list the orbital number (ON), estimated maximum radius ( $R_{\text{Max}}$ ), and maximum mass ( $M_{\text{Max}}$ ) in the Earth radius and mass unit, and the transit probability ( $P_{\text{tr}}$ )

^a Discovery method of the system: ‘Tr’, ‘RV’, and ‘Im’ represent transit, radial velocity, and imaging, respectively.

^b Orbital numbers (ON) followed by ‘H’ indicate the predicted planets within the HZ, and ON followed by ‘C’ indicate that the corresponding orbital periods have been flagged as “Planetary Candidate” in NASA Exoplanet Archive.

Table 4 lists the systems with interpolated and extrapolated planet predictions. In this table, we present the $\chi^2/dof$ before and after interpolation in the fourth and fifth columns. $\gamma$ and $\Delta\gamma$ (where $\Delta\gamma=(\gamma_{1}-\gamma_{2})/\gamma_{2}$ ; $\gamma_{1}$ and $\gamma_{2}$ are the highest and second-highest $\gamma$ values for system, respectively) are listed in columns 6 and 7, respectively. The definitions of other columns in this table are the same as in Table 3. Table 4 is also sorted by the transit probability in descending order.

In Figure 7, we show the new systems containing the detected planets (cyan circles) and the predicted (red circles) planets by extrapolation, while the predicted planets within the HZ are also shown as green circles. The sizes of the symbols are scaled based on the planet’s radius. Similarly, Figure 8 presents the detected exoplanets (cyan circles) and the predicted exoplanets (red circles) by interpolation and extrapolation.

Figure 9 illustrates the radius versus the orbital period of the detected (blue points) and predicted (red points) exoplanets. We find similar trends for both detected and predicted exoplanets. As seen, the vast majority of exoplanets have larger radii and shorter orbital periods than Earth (the radius and orbital period of the Earth have been shown with vertical and horizontal dotted black lines) due to observational limitations.

We note that the estimated masses (and radii) of GJ 676 A, Kepler-56, and WASP-47 are excluded because they are higher than the maximum possible limit of a typical planet. Furthermore, due to the lack of stellar parameters, the transit probabilities of HD 31527, HD 136352, and Kepler-402 are not calculated.

3.2. The reliability of planet predictions using the TB relation

To examine the reliability of the predictions made by the TB relation, we look for some planetary systems in the literature that have had new exoplanets detected recently, notably those systems that have been predicted to have additional planets by BL15. We find the following seven systems with new planet detection: Kepler-1388, Kepler-150, Kepler-1542, Kepler-20, Kepler-80, Kepler-82, and KOI-351.

The detected planets in Kepler-150 and Kepler-82 were not predicted by BL15 (Schmitt et al. 2017; Freudenthal et al. 2019), while the rest of the five planets have been detected with orbital periods that agree with their predicted periods. For Kepler-1388 and Kepler-1542, two planetary candidates with orbital periods of 75.73 and 7.23 d have been detected, the detected periods of which are consistent with their predicted orbital periods, $73.0\pm8.0$ and $6.9\pm0.3$ d, by BL15 and also with our predictions $74.2\pm14.3$ and $6.9\pm0.2$ d. BL15 could also predict an additional planet in Kepler-20 with an orbital period of $39.1\pm5.4$ d. This prediction was also confirmed by Buchhave et al. (2016). Shallue & Vanderburg (2018) used deep learning algorithms to classify potential planet signals, where they managed to validate two more new planets in Kepler-80 and KOI-351, which were previously predicted by BL15. Table 5 summarises these results and presents the predicted and detected periods, which are reasonably consistent, especially for Kepler-1542 and KOI-351 where BL15 highlighted them as predictions with a high geometric probability to transit. This demonstrates the potential and capability of such predictions to help search for new planets, using more precise observational data and new detection techniques such as machine learning. However, having seven systems is not statistically significant enough to establish either the reliability or unreliability of the TB relation, and we require more follow-up observations. We must wait for ongoing exoplanet surveys such as TESS to be completed.

Figure 6. Dynamical spacing $\Delta$ and the total number of adjacent exoplanet pairs that are in orbital resonance with each other. The solid blue line shows the number of resonance pairs considering our predicted planets, and the dotted black line shows the number before inserting any predicted planets into systems. The values for the Solar System are also shown for reference via the orange dashed line. The vertical dash-dotted red line corresponds to $\Delta=10$ and separates the less and more stable adjacent planet pairs regimes.

Table 4. Systems with interpolated and extrapolated planet predictions. Columns 1 and 2 present the host star name and discovery method (Dis.). Column 3 reports a flag that defines whether the system has already been analysed by BL15 (or BL13) (Y) or not (N). Columns 4 and 5 present $\chi^2/dof$ before and after interpolation. Columns 6–10, respectively, list the $\gamma$ , $\Delta\gamma$ , slope (m), intercept (b), and predicted orbital period. Column 11 reports whether the predicted period values in this paper and BL15 (or BL13) are consistent within error (Y) or not (N). Columns 12–15, respectively, list the orbital number (ON), estimated maximum radius ( $R_{\text{Max}}$ ), and maximum mass ( $M_{\text{Max}}$ ) in the Earth radius and mass unit, and the transit probability ( $P_{\text{tr}}$ )

^a Discovery method of the system: ‘Tr’, ‘RV’, ‘OBM’, and ‘PT’ represent transit, radial velocity, orbital brightness modulation, and pulsar timing, respectively.

^b $\Delta\gamma=(\gamma_{1} - \gamma_{2})/\gamma_{2}$ , where $\gamma_{1}$ and $\gamma_{2}$ are the highest and second-highest $\gamma$ values for the system, respectively.

^c Orbital numbers (ON) followed by ‘E’ indicate the extrapolated planets, followed by ‘H’ indicate the predicted planets within the HZ, and followed by ‘C’ indicate that corresponding orbital periods have been flagged as ‘Planetary Candidate’ in NASA Exoplanet Archive.

3.3. Comparison between predictions by BL15 and this study

This study uses the TB relation to predict the probability of finding additional planets in multiplanet systems with at least three member planets, similar to the method used by BL15. In principle, we update the BL15 study and utilise the MCMC simulation. The MCMC method is used to analyse systems and quantify the uncertainties of the best-fit parameters (Foreman-Mackey et al. 2019). The main feature of our method is inserting thousands of hypothetical planets with random orbital periods into each of the systems and using the highest $\gamma$ value to achieve the most accurate predictions.

To compare our method and that of BL15, we remove Ceres and the planet Uranus (two objects predicted by the TB law to exist) from the Solar System and insert thousands of random planets into the Solar System, covering all possible locations and combinations between the two adjacent planets, to achieve the minimum value of $\chi^2/dof$ . This leads us to recover the actual Solar System with highest $\gamma$ value, where the predicted orbital periods’ calculated errors are 20% for Ceres and 23% for Uranus. On the other hand, the BL15 method results in an error value of 23% for both Ceres and Uranus. We remove the planets Earth and Mars, which are within the HZ range of the Sun, and apply the TB relation to the system to predict their orbital periods. This recovers the combination of the actual Solar System where the calculated errors of orbital periods are 18% for Earth and 41% for Mars. We perform this process once again using the BL15 method, which gives us the error values of 18.5% for Earth and 42% for Mars. For a better comparison, we repeat this process for other planets in the Solar System and exoplanetary systems, and with different combinations of planet removal.

We also reapply our method to Kepler-1388, Kepler-1542, Kepler-20, Kepler-80, and KOI-351, regardless of the successfully prediction of planets by BL15, to compare the validity of the methods used in predicting the additional planets. As shown in Table 5, our predictions have fewer errors, and sometimes the predicted periods have smaller error bars. These fewer error bars can be interpreted as advantages of using the MCMC method in predicting exoplanets based on the TB relation. Using the paired samples t-test, the p-value is estimated to equal 0.025, which is less than 0.05 and statistically significant. Therefore, we conclude that, on average, there is evidence that applying the MCMC method does lead to more precise predictions than BL15’s method. The difference between the calculated errors with our method and BL15’s is illustrated in Figure 10, where each blue data point belongs to a planet recovery. Five red data points represent those detected planets after the predictions made by BL15.

Figure 7. Orbital periods and scaled radii of exoplanets in multiple planet systems including the predicted exoplanets from extrapolations. The cyan and red circles indicate detected and predicted planets in systems, respectively. The green circles also indicate the predicted planets located within the HZ of their parent stars. The estimated radius of the predicted planet in GJ 676 A is higher than the maximum possible limit of a typical planet. Furthermore, due to the discovery method of HR 8799, the predicted planet’s radius is not calculated. Therefore, GJ 676 A and HR 8799 are excluded.

Figure 8. Orbital periods and scaled radii of exoplanets in multiple planet systems with interpolated and extrapolated planet predictions. The cyan and red circles indicate detected and predicted planets in systems, respectively. The green circles also indicate the predicted planets within the HZ of their parent stars. Due to the discovery methods of KIC 10001893 and PSR B1257+12, the predicted planets’ radii are not calculated. Therefore, KIC 10001893 and PSR B1257+12 are excluded.

Figure 9. The vast majority of the detected and predicted planets have larger radii and tighter orbits than Earth.

Figure 10. Comparison of uncertainty on the predicted orbital periods calculated by BL15 and this study. We remove planets from systems in various combinations and apply the TB relation to recover the orbital periods of removed planets. Each blue data point belongs to a specific combination of removed planets from systems, and five red data points represent those detected planets after the predictions (see Table 5) made by BL15.

Table 5. Systems with detected planets since predictions made by BL15

4. Summary and Conclusions

Using the MCMC method, we apply the TB relation to all available exoplanetary systems with three or more confirmed exoplanets, a total of 229 systems, to examine their adherence to the TB relation in comparison to the Solar System and to predict the existence of possible additional planets. For those systems that adhere to the TB relation better than the Solar System, we extrapolate one additional planet, and for each of the remaining systems that adhere worse, we interpolate up to 10 specific planets between the detected planets and identify new possible additional planets in the systems.

We present a list of 229 analysed exoplanetary systems (of which 123 of them have not been previously analysed by either BL13 or BL15) containing their unique TB relations and a total of 426 additional predicted exoplanets, of which 47 are located within the HZ of their parent stars. We also estimate that five of the predicted planets in HZ have maximum mass and radius limits within the Earth’s mass and radius range.

As an important result, the planets of $\sim 53\%$ of our sample system adhere to a logarithmic spacing relation better than the planets of the Solar System. Therefore, there is a need to work with more comprehensive data of multiplanetary systems to reveal the probable dynamical or gravitational aspects of the TB relation.

We find new planet detection for seven exoplanetary systems after the predictions made by BL15 and compare the detected and predicted orbital periods. We find that both the detected and predicted orbital periods agree very well within errors. Our predictions also agree roughly better than those made by BL15, indicating that using a precise modelling method and measurements could improve our predictions and uncertainties as well. However, to claim the (un)reliability of the TB relation in predicting the presence of additional planets in exoplanetary systems, we require much more data and further follow-up observations of the exoplanetary systems. Thus, we must wait for upcoming new exoplanet surveys or ongoing surveys such as TESS.