Hidden Worlds: Dynamical Architecture Predictions of Undetected Planets in Multi-planet Systems and Applications to TESS Systems

Our goal is to determine the likelihood of finding another currently unknown planet in a multi-planet system using planetary statistics derived from large populations. The formalism of this likelihood is a triple integral of the probability density function over the inclination, period, and planet radius space.

$$L_{p}=\int_{0}^{180}\int_{R_{min}}^{R_{max}}\int_{P_{min}}^{P_{max}}\Pr(P,R,i)\,dP\,dR\,di$$

(1)

Each of these probability density functions describes how likely it is a planet exists with those parameters. For example, $\Pr(P)$ would be the probability density function for finding a yet-unknown planet in a system with period $P$, etc. We make the assumption that the probability function $\Pr(P,R,i)$ is separable in the three dimensions; that is, the probability distributions w.r.t. period, planet radius, and inclination are not dependent on the other two parameters (M18). We explore and justify this assumption in Section 2.2. Thus,

$$\Pr(P,R,i)=\Pr(P)\Pr(R)\Pr(i),$$

(2)

and we will now consider the three components of the likelihood function separately. The period distribution for a system with known planets is:

$$\Pr(P)=\left\{\begin{array}[]{ll}1,&P_{j}=P_{k}\textrm{ and }j=k\\ 0,&P_{j}=P_{k}\textrm{ and }j\neq k\\ f(P),&P_{j}\neq P_{k}\textrm{ and }P_{min}<P_{j}<P_{max}\\ \end{array}\right.$$

(3)

where $P_{j}$ is the period of planet $j$ in the system (known or hypothetical), $P_{k}$ is the orbital period of known planet $k$ in the system, $f(P)$ is defined by a probability density function (PDF), and $P_{min},\,P_{max}$ are the limits on the period. The probability of finding another planet at the orbital period of a known planet is 0, since there’s no stable configuration where this can occur.

The probability distribution for the planetary radius is:

$$\Pr(R)=\left\{\begin{array}[]{ll}1,&R_{j}=R_{k}\textrm{ and }j=k\\ f(R),&(R_{j}=R_{k}\textrm{ and }j\neq k)\textrm{ or }(R_{j}\neq R_{k}),\\ &R_{min}<R<R_{max}\\ \end{array}\right.$$

(4)

where $R_{j}$ is the planetary radius of planet $j$, $R_{k}$ is the radius of known planet $k$, $f(R)$ is defined by a PDF of planet sizes and $R_{min},\,R_{max}$ are the planet radius limits. The posterior probability of finding a planet with radius $R_{k}$ is 1, since we know at least one planet in said system has that planetary radius, but the probability of finding another planet with that radius is determined from the PDF as there are no exclusionary limits on the planet radius.

Similarly, the probability distribution function for the inclination is:

$$\Pr(i)=\left\{\begin{array}[]{ll}1,&i_{j}=i_{k}\textrm{ and }j=k\\ f(i),&(i_{j}=i_{k}\textrm{ and }j\neq k)\textrm{ or }(i_{j}\neq i_{k}),\\ &0<i<180\\ \end{array}\right.$$

(5)

where $i_{j}$ is the inclination of planet $j$, $i_{k}$ is the inclination of known planet $k$, and $f(i)$ is defined by a PDF of planet inclinations. The inclinations are drawn from 0^∘ (face-on prograde) to 180^∘ (face-on retrograde), with 90^∘ being edge-on. The system inclination, which is the plane on which the orbits are centered, is close to 90^∘ for systems with transiting planets. The mutual inclinations between planets in the system are drawn from a Rayleigh distribution with a relatively small scatter centered around the system inclination, which is a parameter to the function that can be inferred from fitting the known planet inclination distribution (Fang & Margot, 2012; Fabrycky et al., 2014).

Once the probability function is determined across all three dimensions, we sample it via the Monte Carlo method. We take $N=10,000$ iterations of each system and inject planets with periods, radii, and inclinations determined by sampling the probability function. The systems with the highest relative likelihood in the shortest periods can then be prioritized for ground-based follow-up observations.

An approximation made in the current version of DYNAMITE is that the period, planet radius, and orbital inclination of the planet are independent of each other. In the following we will explore the justification for this approximation.

Multiple works (e.g., Helled et al., 2016; Carrera et al., 2018) have shown that short-period planets tend to be smaller than planets further out, as the host star evaporates any light gas envelope. This trend is most significant for large planets ($R>4-5R_{\oplus}$; Ciardi et al. 2013, Huang et al. 2016), which we exclude from our sample. Weiss et al. (2018) and M18 show that for Kepler targets, planets within the same system tend to be the same size compared to random draws. Weiss et al. (2018) found a correlation between pairwise inner and outer planet radii, but Gilbert & Fabrycky (2020) showed this correlation is unclear and may be due to scatter within low-monotonicity systems. Dawson et al. (2016) showed that flatter systems with smaller mutual inclinations tend to have smaller planets with closer orbital spacing, whereas systems with larger mutual inclinations tend to include gas giants with larger orbital spacings. Gilbert & Fabrycky (2020) found that the mass scale (and thereby the radius scale) of a system does not correlate to the flatness, but does reproduce the slightly weak but significant correlation of tighter spacing in flatter systems.

Although evidence exists for correlations between fundamental planet parameters, given current data these correlations are relatively weak, and the mechanisms behind these correlations are not yet fully understood. Given that planet properties within a system are only weakly correlated and that there is no well-established, quantified correlation function we could implement, in our study we take the assumption that the parameters are independent. This, we stress, is an approximation, one which can be tested by the predictions given out by DYNAMITE and refined when a more accurate population model can quantify these correlations across all the variables. DYNAMITE would then be updated to include these correlations in the formalism.

We focus our study on the orbital period range from 0.5 to 730 days, for which exoplanet population data is available (see, e.g., M18). The occurrence rates of planets are not well defined interior to 0.5 days and beyond 2 years. It is extremely difficult to find long-period transiting planets due to the amount of time needed as well as the low probability of having the correct inclination to transit; as the simple geometric probability of transiting goes as $R_{*}/a$, the inclination limits for a transiting planet are closer to 90^∘ for a planet at a larger semi-major axis $a$.

In this study we assume planet radius limits from 0.5 to 5 $R_{\oplus}$, as a large majority of the TESS multi-planet systems have planet radii that fall between these boundaries. It has been noted that planets with radii above 3-4 $R_{\oplus}$ are below an observed “break” in the planet radii distribution across all systems, and thus may have their own statistical distribution (H19 and references therein). The Kepler sample does not constrain populations below 0.5 $R_{\oplus}$, as they are often too difficult to detect en masse. Smaller, rocky (super)Earth-like planets are more common in the galaxy than Jupiter-like gas giant planets within 1 AU (e.g., Howard et al., 2012; Mulders et al., 2015), and finding them were direct mission goals for TESS and Kepler, as well as ground-based follow-up observations that are still ongoing and the main reason we prioritize systems likely to contain these planets in our statistical analysis.

We assume the inclination of a planet can be any value between 0^∘ and 180^∘ but is constrained by the PDF of our distribution. However, since the measurement of the inclination is degenerate (i.e., we receive the same signal for planets with inclination $i$ and $180-i$), we assume the inclination of the first planet is $<90$ degrees. Also, for a high number of known transiting planets ($N\geq 4$), we assume the planet inclinations follow the Rayleigh distribution only, whereas for $N=2,3$ we assume an intrinsic scatter where some fraction of systems have inclinations distributed isotropically as per the “Kepler dichotomy” (Johansen et al. 2012, M18) but where the fraction is scaled downwards from the isotropic fraction of singly-transiting planets. For example, if the isotropic fraction of systems with one known planet is 30%, then the isotropic fraction of systems with two planets is 20% and three planets 10%. This is done so as to not have a large break in the inclination distribution between singly-transiting systems and doubly-transiting systems if the isotropic fraction is high; the odds of a misaligned planet in a doubly-transiting system should be higher than a misaligned planet in a 4-transiting-planet system (Ballard & Johnson, 2016).

To ensure that the new systems created after injections are dynamically stable on long timescales, we use the dynamical stability parameter $\Delta$, which is defined as

$$\Delta=\frac{2(a_{2}-a_{1})}{a_{2}+a_{1}}\left(\frac{3M_{*}}{m_{2}+m_{1}}\right)^{1/3},$$

(6)

where $a_{1},a_{2}$ are the semi-major axes for the orbits of the inner and outer planet respectively, $m_{1},m_{2}$ are the inner and outer planet masses, and $M_{*}$ is the stellar mass. The innermost theoretical limit for $\Delta$ is $\Delta_{c}=2\sqrt{3}\approx 3.46$ (Gladman, 1993), but many studies have shown that long-term stability requires $\Delta_{c}\gtrsim 10$ for systems of three or more planets with similar masses/sizes, like those seen in the Kepler dataset (Chambers et al., 1996; Fang & Margot, 2012; Pu & Wu, 2015).

For our minimum threshold for planet dynamical stability, we take the more inclusive value $\Delta_{c}=8$ from H19 in order to allow for a larger area of period space to be tested, as some known systems have smaller $\Delta$ values (e.g., the two planets of Kepler-36; Carter et al. 2012). In most cases for planets discovered via the transit method, the planet masses are unknown, so we use the non-parametric mass-radius relation from Ning et al. (2018) to formulate a range of masses corresponding to the range of planet radii for the injected planets.

We defined two separate period PDFs to study and understand the impact of the assumption on the results of fitting current systems to different model distributions and testing the resulting likelihoods. We first used the period ratio model from M18. The period ratio model assumes a broken power-law for the location in period space of the first planet in the system with the break at 12 days, where the peak is in the first-planet statistics from Kepler (M18). Outwards from the first planet, it calculates a dimensionless spacing between consecutive planets $D_{k}$, defined as

$$D_{k}=2\frac{\mathcal{P}_{k}^{2/3}-1}{\mathcal{P}_{k}^{2/3}+1},$$

(7)

where $\mathcal{P}_{k}=P_{i+1}/P_{k}$. The period ratio model assumes $D_{k}$ follows a lognormal distribution with parameters $\mu_{D}$ and $\sigma_{D}$. Thus, the full period distribution is

$$f(P)=\left\{\begin{array}[]{ll}1,&P=P_{j}\\ (P/P_{break})^{a_{P}},&p_{k}=1\textrm{, }P<P_{break},\\ &\textrm{and }P<P_{1}\\ (P/P_{break})^{b_{P}},&p_{k}=1\textrm{, }P>P_{break},\\ &\textrm{and }P<P_{1}\\ \frac{1}{\sqrt{2\pi\sigma_{D}^{2}}}e^{-\frac{(\log D_{k}-\mu_{D})^{2}}{2\sigma_{D}^{2}}},&p_{k}>1\\ 0,&\Delta_{jk}<\Delta_{c}\\ \end{array}\right.,$$

(8)

where $P_{j}$ is the period of the $j$th known planet in the system (so $P_{1}$ is the period of the innermost known planet), $p_{k}$ is the position of a new planet counting outwards from the star, $a_{P}$ and $b_{P}$ are the power-law parameters, and $\Delta_{jk}$ is the stability criterion measured between known planet $j$ and inserted planet $k$. Values for $\mu_{D}$, $\sigma_{D}$, $a_{P}$, and $b_{P}$ are found in Table 1. In particular, this model tends to insert planets roughly symmetrically in gaps, as it is not sensitive to the planet masses. However, the stability criterion relies on the masses and cuts off the period probability distribution at a further distance from a more massive planet.

Secondly, we used the clustered periods model from H19. To determine the cluster periods (the locations in period space where known planets in a system tend to be found), the clustered period model iterates through a set of test periods and determines the fit of the known planetary periods in the system (scaled by the current test cluster period), to a lognormal distribution,

$$f(P_{i})\propto P_{c}\times\textrm{Lognormal}(0,N_{P}\,\sigma_{P}),$$

(9)

with $N_{p}$ being the number of planets in the cluster and $\sigma_{P}=0.2$ (H19). If planets in the system have scaled periods beyond the 97th percentile of the lognormal distribution that fits the remaining planets best (a variation of $\sim 2.5\sigma$), we assume a two-cluster model where the probability of finding a planet at a certain period can be in one of two clusters. This can be extended out beyond two clusters, but for almost all systems in the TESS multi-planet catalog the two-cluster model is sufficiently complex to describe the orbital hierarchies, as most of the systems only contain two planets. Thus, our full period distribution for the clustered periods model is

$$f(P)=\left\{\begin{array}[]{ll}1,&P=P_{j}\\ P_{c,i}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(\log P)^{2}}{2\sigma^{2}}},&P/P_{c,i}<2.5\sigma\\ &\textrm{ and }\sigma=(N_{P})_{i}\sigma_{P}\\ 0,&\Delta_{jk}<\Delta_{c}\\ \end{array}\right.,$$

(10)

where $P_{c,i}$ is the cluster period for the $i$th cluster (clusters separated by 2.5$\sigma$), $(N_{P})_{i}$ is the number of planets in that cluster, and $\Delta_{jk}$ is the stability criterion measured as before. This model tends to insert planets close to the stability criterion and is more sensitive to the masses of the known and injected planets, as the injected planets are closer in period-space to the lower-mass planets. DYNAMITE is most effective at calculating the likelihood of planets interior to the outermost known planet in a system, as actual planet occurrence is unknown at long periods. DYNAMITE prioritizes systems with gaps and large likelihoods for planets in those gaps.

For our planet radius distribution, we use a variation on the clustered periods and sizes model from H19. Here, the clusters are locations in planet radius space where known planets in a system have similar values. We fit all the known planetary radii in the system to lognormal distributions for a range of cluster radii,

$$f(R_{i})\propto\textrm{Lognormal}(R_{c},\sigma_{R}),$$

(11)

with $\sigma_{R}=0.3$ (H19) and finding the cluster radius or radii where the sets of known planets fit best. The number of clusters is determined by the 97th percentile of the lognormal distribution of the first best-fit cluster; if the planetary radius of a known planet in the system is beyond this value, we conclude that the planet is in a separate cluster and adjust our PDF accordingly. We do not assume that periods and planet radii are clustered together, as that would invalidate our assumption that all three parameters are independent.

Inclinations are taken from a broken distribution similar to M18. In this distribution,

$$f_{iso}=0.38\times\frac{4-x}{3},$$

(12)

where $x=\textrm{min}(4,N)$ and $N$ is the number of planets in the system. A fraction $f_{iso}$ of systems have inclinations distributed isotropically (and when $x=1$, $f_{iso}=38\%$; M18), otherwise they are chosen from a Rayleigh distribution centered on the system inclination $i_{s}$.

$$f(i_{i})\propto\left\{\begin{array}[]{ll}i_{s}+\textrm{Rayleigh}(\sigma_{i}),&i_{iso}>f_{iso}\\ \sin(i_{i}),&i_{iso}\leq f_{iso}\\ \end{array}\right.$$

(13)

with scale factor $\sigma_{i}=2$ (Fabrycky et al., 2014).

The system inclination $i_{s}$ is determined a similar way to the cluster periods and radii, by fitting a Rayleigh distribution with scale factor $\sigma_{i}$ to the known planet inclinations. After constraining the first planet inclination to $<90$ degrees, we permute all the remaining planets across the degeneracy. Whichever permutation of planet inclinations that fit the Rayleigh distribution best for a given $i_{s}$ defines that as the system inclination.

We determine the dynamical stability threshold $\Delta_{c}$ for the known planets in the system compared to an average mass planet for that system (either calculated from the planetary mass values or taken from the average planet radius and the mass-radius relation), and we set the probability of finding a planet at the corresponding periods where $\Delta_{c}<8$ to 0. A full list of parameters in DYNAMITE and the values assigned to them can be found in Table 1.

First, we tested removing a known transiting planet from a system, choosing systems with high multiplicity as they provide the most detailed prior information for our predictions. Kepler-154 (catalog ) is a 6-planet system with 5 smaller planets in a short-period orbital chain (periods from 4 to 62 days) and a gas giant much further out ($P\sim$ 3.4 years; Rowe et al., 2014; Morton et al., 2016; Kawahara & Masuda, 2019). We run the DYNAMITE on the inner transiting planets of Kepler-154 (i.e., assuming the presence of 3.4-year-period gas giant does not affect the interior planets and removing it from the system). Also, we exclude the last verified candidate KOI-435.06 = Kepler-154 f with a period of 9.92 days, as it was the transiting candidate with the lowest confidence/highest false alarm probability (Morton et al., 2016, see Figure 1).

When running DYNAMITE on the four remaining planets, we find that a large percentage (96.7%) of injected planets interior to the outermost planet are found between the first and second planets (Kepler-154 e at 3.93 days and Kepler-154 d at 20.55 days; uncertainties on the periods measured from Kepler are $<10^{-4}$ days). 90.1% of interior injections occur within 5 days of 9.92 days, with a mode at 9.9 days that is consistent with the planet period within the fineness of our model. The radius of Kepler-154 f is measured as $1.50^{+0.34}_{-0.21}R_{\oplus}$, and 67% of DYNAMITE’s predicted planets have a physical radius within 3$\sigma$ of the known radius.

The orbital inclination of Kepler-154 f is not fully constrained, but using the average impact parameter at its period of 9.92 days places it at an inclination of $88.6\pm 0.2^{\circ}$, which is near the calculated plane of the system (as the remaining transiting rocky planets have inclinations of $85.0\pm 0.2^{\circ}$, $88.5\pm 0.3^{\circ}$, $89.5\pm 0.5^{\circ}$, and $89.4\pm 0.4^{\circ}$, as calculated from their impact parameters in Morton et al. (2016)). DYNAMITE predicts $88.9\pm 0.8^{\circ}$ as the orbital inclination of the injected planet, matching the predicted inclination well but with a larger spread, as 43% of the iterations from the model are within 3$\sigma$ of the estimated inclination for Kepler-154 f. As the mean injection value from DYNAMITE matches the estimated value of the orbital inclination very well, we find this agreement to be very strong. Good predictions need only match the inclination within a couple of degrees, to identify whether transits are likely and to constrain the degeneracy on $M\sin i$ for RV planets. A comparison of the known data with the predictions from DYNAMITE is shown in Figure 2.

Kepler-20 (catalog ) is a 6-planet system with all 6 planets in an orbital chain with periods from 3.6 days to 77.6 days (Fressin et al., 2012). The last discovered planet was Kepler-20 g, which has a period of $34.94\pm 0.04$ days and was found not to transit (Buchhave et al., 2016). By removing this planet, we can test the system as it was seen only through transit light curves, or as if there were no radial velocity (RV) data confirming the existence of this planet. We ran DYNAMITE on only the transiting planets of Kepler-20 (see Figure 3).

We find that there should be a very high probability (99.7% of total possible interior injections) of a planet in between the farthest planet (Kepler-20 d at 77.61 days) and second-farthest planet (Kepler-20 f at 19.58 days), the gap in period space where Kepler-20 g was found. Narrowing the injections down to a 10-day window centered on Kepler-20 g’s known period, we still find that 42.4% of interior injections occur in this area, with a mode at 39.0 days. In addition, injected planets located in this gap with the estimated planet radius for Kepler-20 g transit only 22.6% of the time, showing that it is not an outlier because it does not transit. Due to its non-transiting status, a precise planetary radius for Kepler-20 g has not been measured, but we infer a planet radius of $3.46^{+0.30}_{-0.25}R_{\oplus}$ from the non-parametric mass-radius relationship. DYNAMITE determined that the planet radii of this system likely follow a 2-cluster model, and predicts the planet radius of an injected planet in that cluster to be within 3$\sigma$ of the inferred planet radius 55% of the time.

In both tests of this scenario, we find that DYNAMITE predicts the presence of the known transiting or non-transiting planet that was removed, and is more accurate with the period ratio model than the clustered periods model. For Kepler-154, DYNAMITE predicts the location of the planet in period space almost exactly, whereas for Kepler-20, DYNAMITE predicts a longer period than the expected period. However, the presence of Kepler-20 g could still be inferred, and follow-up observations (both transit light curves and precision RV) would prioritize that area in period space for its discovery. The primary function of DYNAMITE is to provide prioritization for follow-up targets and not to predict the parameters exactly. Our tests show that DYNAMITE predicts not only the periods but also the radii and inclinations of planets to levels that are well-suited for guiding observations. In the future a single, combined metric based on these multiple dimensions could be created to allow for, for example, optimizing the assumptions that underpin DYNAMITE.

In this scenario, we take a 4-planet system and remove 2 planets in various configurations. TOI 174, also known as HD 23472 (catalog ), has been confirmed to have 2 sub-Neptune-sized planets ($\sim 2R_{\oplus}$) in $17.667^{+0.142}_{-0.095}$ day and $29.625^{+0.224}_{-0.175}$ day orbits (Trifonov et al., 2019), which are called HD 23472 b and c. TESS has also found evidence for two possibly rocky planets ($\sim 1-1.3\textrm{ }R_{\oplus}$) interior to the known planets at 4.0 and 12.2 days, called TOI 174.04 and TOI 174.03, respectively. Removing both TOI 174.03 and TOI 174.04 from the system, we find that all interior injections occur at orbits shorter than the period of HD 23472 b, and the mode of planet injections for the entire system (interior and exterior) occurs at the Kepler statistical first-planet mode of 12 days (Mulders et al., 2018), showing a very high likelihood of at least one planet interior to HD 23472 b at a 17.7 day period.

When removing TOI 174.04 and HD 23472 b, we find an almost exactly even distribution in the location of interior injections; 50.2% are injected at orbits shorter than TOI 174.04 at 12.2 days, and 49.8% are injected at orbits between TOI 174.04 and the remaining known planet HD 23472 c at 29.6 days, with a mode at 18.5 days. When removing TOI 174.03 and HD 23472 b, we instead find a wide distribution of planet injections with a mode of 12 days. The total probability of injecting one planet interior to HD 23472 c is only slightly more than half of the previous case, however, suggesting that injecting multiple planets here would fit the data better (see Figure 4).

When we remove TOI 174.03 and HD 23472 c, there is a high probability of finding two planets again, but not necessarily in the same positions as currently known. The mode of the injections between TOI 174.04 and HD 23472 b is around 8.5 days, as two period ratios near 2 are more likely than period ratios of $\sim$3 and $\sim$1.5, respectively. The mode of the injections exterior to HD 23472 b, following the Kepler statistics, is $\sim$32 days, but there is still a 46% probability of finding the outermost planet within 5 days of its current period.

We show the different system structures after removal with relative likelihoods in log period space in Figure 5. In all tests of removing two out of four planets in this system, the two-cluster planet radius model failed and the two remaining objects fell into one cluster. When removing one of the TESS planet candidates along with one of the known planets, the cluster centers in between the two, whereas when both TOI 174.03 and TOI 174.04 are removed the cluster shifts towards the planet radii of HD 23472 b and c, and vice versa. This would affect the prioritization to a degree, depending on what types of planets the follow-up observations are looking for. The orbital inclinations of the planets all fell within 1^∘ of each other, so removing any pair of them did not affect the plane of the system significantly.

We also tested removing one planet from the TOI 174 system and perturbing it to see if DYNAMITE is rigorous to changes in the system. We altered the period of each of the remaining planets by a random amount within 3$\sigma$, as well as 5%. Neither perturbation altered the predictions of the mode of the injected planets, within the fineness of the model at 0.5 days; the histograms did shift slightly but insignificantly when perturbations between remaining planets occurred in the same radial direction. The planet radius and inclination predicted values also shifted similarly.

DYNAMITE is also able to provide analysis for planet candidates that are uncertain. TOI 1469 was previously identified as a planetary system around the K3V star HD 219134; it contains 2 rocky super-Earth transiting planets that are analyzed in this sample, and at least 3 non-transiting massive planets exterior to these. This system has been studied extensively as a result, and an unconfirmed planet candidate with its validity in doubt (as it was found in only one of these four studies of this system) was predicted with a period of $94.2\pm 0.2$ days (Motalebi et al., 2015; Vogt et al., 2015; Johnson et al., 2016; Gillon et al., 2017). Running DYNAMITE in an iterative way places injected planets most often at 12.5, 23, and then 45.5 days, adding in an injected planet every iteration; HD 219134 f and d have periods of $22.72\pm 0.02$ and $46.86\pm 0.03$ days, respectively. In addition, DYNAMITE predicts a fourth injected planet at $\sim$ 87 days given the general period ratio sample. These model runs are shown in Figure 6.

As a planet at 12.5 days actually minimizes the standard deviation of the period ratios in this system, if we then also insert a planet beyond HD 219134 f using the mean of the period ratios, as M18 states the period ratios tend to be consistent within a system, then the inserted planet is most often found at a period of 93 days. This inserted planet would only have a 5.1% chance of transiting, given that the closer two planets do not transit. The outermost known planet in the system has a period of $\sim$ 6 years so it would not be found by our analysis, as we only search out to 2 years due to the statistical constraints imposed by the current datasets and models. Based on the planet radius and inclination analysis, the missing planet at 12.5 days would likely be in the smaller cluster (similar in size to the rocky super-Earths interior to it) and have a 75% chance of transiting. However, the transits of a 1.55-Earth-radius planet around a 0.78-solar-radius star would have a depth of $\sim$ 330 ppm and would likely only be visible from space (Gillon et al., 2017).

We tested the sensitivity and robustness of DYNAMITE by altering known multi-planet systems and seeing what was predicted when DYNAMITE was run on these systems. We found that DYNAMITE responded well to the removal of one planet in the system, although the underlying statistical distributions were still noticeable if the original system wasn’t regularly dynamically spaced. Removing multiple consecutive planets caused a drop in the relative likelihood of injecting one planet across the large gap, but we showed that injecting two planets across that gap would recover the planets and near their correct periods as well. The planet radius and inclination were mostly stable, although if one planet radius cluster consisting of two planets was removed from a system, then DYNAMITE would accordingly shift the relative likelihood away from finding planets of that size. All three parameters were robust to changes on the order of a few percent.

A compelling use of the predictive power of DYNAMITE as shown in Section 3 is to analyze recently-discovered multi-planet systems from TESS and find the systems with the highest relative likelihood of containing another detectable planet. TESS has a unique parameter space for exoplanet transit searches, in that it covers the entire sky in one-month sectors with magnitude limits of $m_{T}<13$ for super-Earth transits and $m_{T}<15$ magnitude for precise stellar photometry (Ricker et al., 2015). This method can be used to predict the location of possible unseen planets in these systems when looking through the archival TESS data for low-confidence signals which could be possible planets, as well as which systems are the best targets for follow-up ground-based observations to find new planets. The TESS extended mission will discover many more multi-planet systems that can be probed for these hidden planets.

To create our sample, we took the known multi-planet systems identified by TESS as Objects of Interest (TOIs) located in the ExoFOP-TESS archive. From the total list, we removed 4 systems containing planets or planet candidates with a physical radius larger than 5 $R_{\oplus}$, as these giant planets ($>5-6R_{\oplus}$) tend to have a different distribution in period and multiplicity within a system (Steffen et al., 2012; Dong & Zhu, 2013). This left us with 52 systems containing at least two transiting planets or planet candidates with planet radius less than 5 $R_{\oplus}$.

A majority of the TESS systems have two currently known transiting planet candidates, with a few containing three and only two known systems with four transiting planet candidates. Some systems (e.g., TOI 1469, as shown in Section 3.4) are known to have more non-transiting planets. However, for the purpose of this analysis we used the method on only the transiting planets to gain inclination and planet radius distributions as well as the information on the likely orbital periods for the system.

We ran both orbital period models described in Section 2.4 to fully probe the period space for each system, along with the standard clustered planet radius distribution and Rayleigh inclination distribution with a fraction of isotropic systems. In the period ratio model, if the largest peak in the period distribution interior to the outermost planet is comparable in relative likelihood to the exterior period peak after the last known planet, then we use the inner local maximum as the predicted period. In the clustered period model, we find the specific period inside of one of the clusters that gives the maximum relative likelihood for a new planet. For the planet radius distribution, we simply take the median of the distribution as the predicted planet radius, which is accurate for single clusters and tends to evenly weight the contributions from two clusters in the planet radius.

Given the predicted period, planet radius, and inclination distribution, we then calculate the probability that the planet would transit. In addition, given the predicted planet radius and the host star’s stellar radius, we calculate the most likely absolute transit depth. We use a combination of these two factors to determine which systems should have the highest priority for follow-up, either from more thorough TESS data searches or from ground-based observatories.

Using the period ratio model from M18, we find that most systems have their peak in the period distribution between 10 and 40 days, with a relatively smooth log gradient in both directions away from the center of the peak. The planet radius relative likelihood is high between 1-3 $R_{\oplus}$ and decreasing on both sides. The full gradient across all systems in both dimensions for this model is shown in Figure 7, with the predicted periods and planet radii for each system marked. The full results, including 1$\sigma$ uncertainties calculated from the Monte Carlo iterations, are found in Table 2. The relative likelihood for each system in log period space calculated by the period ratio model is shown in Figure 8, as well as the relative likelihood in planet radius space (consistent between the different period models) in Figure 9. The transit depth and transit probability for the period ratio model are shown in Figure 10.

Using this model, we find for each system that the relative likelihood of finding a planet at a given period tends to peak at the center of a gap between two known planets, especially if the period ratio between those two is $\sim 4$, as well as a factor of $\sim 2$ exterior to the outermost known planet, as shown in Figure 8. The likelihood measured vs separation in mutual Hill radii peaks at 20, similar to the analysis by Gilbert & Fabrycky (2020). For large period ratios between planets, the normalization of the Monte Carlo iterations to 1 injected planet smooths out the distribution and significantly lowers the local maximum in the relative likelihood. The planet radius relative likelihoods show a relatively wider spread in the distribution between $\sim 1-5R_{\oplus}$.

The clustered periods model from H19 tends to predict planets closer to their host stars, with peaks between 1 and 20 days (see Figure 11). As a majority of these systems have a planet within 10 days, this model prefers to place another planet within that cluster. Due to the varied positions of the second planet in the system, the clustered periods model has non-negligible probabilities of forming a second narrow cluster or combining the two planets into one cluster, which shows through in the period predictions. The log period and planet radius parameter space for this model is shown in Figure 11), with the predicted periods and planet radii for each system marked. The full results, including 1$\sigma$ uncertainties calculated from the Monte Carlo iterations, are found in Table 3. The relative likelihood for each system in log period space calculated by the clustered periods model is shown in Figure 12, and the transit depth and transit probability for this model are shown in Figure 13.