ARE PLANETARY SYSTEMS FILLED TO CAPACITY? A STUDY BASED ON KEPLER RESULTS

Our methods for deriving the intrinsic dynamical spacing of planetary systems are as follows. First, we created model populations of planetary systems obeying different underlying distributions of dynamical spacing. Second, we performed synthetic observations of the planetary systems in these populations by the Kepler spacecraft. At this stage we identified which simulated planets were detectable by the Kepler telescope and which were not. Third, we compared the resulting distribution of dynamical spacing of detectable planets from synthetic populations with that of the actual Kepler detections. The actual distribution can be easily obtained from Kepler transit data with an assumed planet radius–mass relationship. Most of these steps are fully explained in Fang & Margot (2012a), and we refer the reader to that paper for details. In the following paragraphs, we summarize the most salient points of our procedure as well as any differences with Fang & Margot (2012a).

Each model population consists of about 10⁶ planetary systems, and we created various model populations that followed different underlying distributions of multiplicity, inclinations, and dynamical spacing. To generate these populations, we needed to restrict the range of physical and orbital properties of the stars and planets that we considered in our simulations. We selected ranges that would adequately overlap those of a Kepler sample that can be considered nearly complete (Howard et al. 2012; Youdin 2011). Stellar properties such as radius R_*, stellar noise σ_*, effective temperature T_eff, surface gravity parameter log(g), and Kepler magnitude K were randomly drawn from the Kepler Input Catalog (see Fang & Margot 2012a). We only considered bright solar-like stars that obeyed the following ranges:

$$\begin{eqnarray} &&4100\,{\rm K} \le T_{\rm eff} \le 6100\,{\rm K}, \nonumber\\ &&4.0 \le \log (g\ [{\rm cm s}^{-2}])\le 4.9, \nonumber\\ &&K \le 15\,{\rm mag}. \end{eqnarray} \tag{ 1 }$$

Planet radii and orbital periods were drawn from debiased distributions, and we obtained these debiased distributions by converting the observed sample of Kepler Objects of Interest (KOI, based on detections up to Quarter 6 released in 2012 February; Batalha et al. 2012) into a debiased sample using calculations of detection efficiencies (see Fang & Margot 2012a). We filtered the KOI sample (and correspondingly limited the parameter space of the synthetic populations described in this paper) to the following orbital period P, planet radius R, and signal-to-noise ratio (S/N) boundaries:

$$\begin{eqnarray} &&P \le 200\,{\rm days}, \nonumber\\ &&1.5 \,R_{\oplus } \le R \le 30 \,R_{\oplus }, \nonumber\\ &&{\rm S/N}{\rm (\rightarrow Q8)} \ge 11.5. \end{eqnarray} \tag{ 2 }$$

These limits were imposed in order to choose a sample of planets with properties unlikely to be missed by the Kepler detection pipeline. For S/N, we required an S/N ⩾ 10 for observations up to Quarter 6, which corresponds to about S/N ⩾ 11.5 for observations up to Quarter 8 by assuming that S/N roughly scales as $\sqrt{N}$, where N is the number of observed transits. This scaling is performed because the S/Ns of observed KOIs have been reported for observations up to Quarter 8 in Batalha et al. (2012), whereas the actual detections have been reported up to Quarter 6 only. Planet masses M were calculated by converting from planet radii R. We used a broken log-linear M(R) prescription obtained by fitting to masses and radii of transiting planets (see Fang & Margot 2012a):

$$\begin{eqnarray} \log _{10} \left(\frac{M}{M_{\rm Jup}}\right) &=& 2.368 \left(\frac{R}{R_{\rm Jup}}\right) - 2.261 \nonumber \\ &&\quad {\rm for } \left(\frac{R}{R_{\rm Jup}}\right) < 1.062, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \log _{10} \left(\frac{M}{M_{\rm Jup}}\right) &=& {-}0.492 \left(\frac{R}{R_{\rm Jup}}\right) + 0.777 \nonumber\\ &&\quad {\rm for } \left(\frac{R}{R_{\rm Jup}}\right) \ge 1.062. \end{eqnarray} \tag{ 4 }$$

Additionally, we repeated all of the methods described in this section by using an alternate mass–radius relationship: (M/M_⊕) = (R/R_⊕)^2.06 (Lissauer et al. 2011b). By adopting this alternate mass–radius equation, our results showed the same best-fit dynamical spacing distribution as defined in Equation (7) with σ = 14.5 (see results presented in Section 2.2). We note that both of these mass–radius equations were obtained by fitting to the sample of planets with known masses and radii. Errors in the mass–radius relationships can potentially affect our results because our determination of dynamical spacing is a direct function of planetary masses. In Figure 1, we investigate how uncertainties in the mass–radius relationship map into dynamical spacing uncertainties. Specifically, we plot histograms showing how the observed dynamical spacing changes if there is a 1σ increase or decrease in mass for the mass–radius equation. For masses lower than nominal, adjacent planets appear to be less closely spaced and so the distribution shifts to the right. For masses higher than nominal, adjacent planets appear to be more closely spaced and so the distribution shifts to the left. The shifts are quantified at the end of Section 2.2.

Zoom In Zoom Out Reset image size

For orbital eccentricities, we adopted circular orbits, as we did in Fang & Margot (2012a). Eccentricities do not directly affect our calculation of dynamical spacing, as we will define below in Equation (5).

Regarding the multiplicity distribution in our model populations, we used a bounded uniform distribution with λ = 1.5–3.5 with increments of 0.25 to assign the number of planets per system. A bounded uniform distribution has a single parameter λ and is defined as follows: first, draw a value N_max (maximum number of planets) from a Poisson distribution with parameter λ that ignores zero values, and second, draw the number of planets from a discrete uniform distribution with range 1–N_max (Fang & Margot 2012a). Thus, each planetary system will have at least one planet. For the inclination distribution of planetary orbits, we used a Rayleigh distribution with σ = 1, 2° as well as a Rayleigh of Rayleigh distribution with σ_σ = 1, 2°. A Rayleigh of Rayleigh distribution has a single parameter σ_σ and is defined as follows: first, draw a value σ from a Rayleigh distribution with parameter σ_σ, and second, draw a value for inclination from a Rayleigh distribution with parameter σ (Lissauer et al. 2011b). These specific multiplicity and inclination distributions were chosen because they yielded fits consistent with transit numbers and transit duration ratios from Kepler detections (see Fang & Margot 2012a). Combinations of these specific multiplicity and inclination distributions add up to a total of 36 possibilities.

The difference between model populations generated in Fang & Margot (2012a) and the model populations generated in this study is the treatment of planetary separations, since here we wish to determine the underlying dynamical spacing of planetary systems. We used a separation criterion Δ to assess the dynamical spacing between all adjacent planets in multi-planet systems, where Δ is defined as (e.g., Gladman 1993; Chambers et al. 1996)

$$\begin{eqnarray} &&\Delta = \frac{a_2-a_1}{R_{H1,2}}, \end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray} &&R_{H1,2} = \left(\frac{M_1+M_2}{3M_*} \right)^{1/3}\frac{a_1+a_2}{2}. \end{eqnarray} \tag{ 6 }$$

In these equations, a is the semi-major axis, R_{H1, 2} is the mutual Hill radius, and M is the mass. Subscripts *, 1, and 2 refer to the star, the inner planet, and the outer planet, respectively. For a two-planet system not in resonance, the planets are required to be spaced with Δ ≳ 3.46 in order to be Hill stable.

In our model populations, adjacent planets in multi-planet systems were spaced according to a prescribed Δ distribution. We used a shifted Rayleigh distribution, which is the same as a regular Rayleigh distribution except shifted to the right by 3.5 (since we require this distribution to provide values of Δ that meet the minimum Hill stability limit). Such a distribution, as we will show, matches the observed sample well. The mathematical form of a shifted Rayleigh distribution f is

$$\begin{eqnarray} &&f(\Delta) = \frac{\Delta - 3.5}{\sigma ^2} e^{-(\Delta - 3.5)^2/(2\sigma ^2)}, \end{eqnarray} \tag{ 7 }$$

and is described by a single parameter σ. In our model populations, we explored values of σ = 10–20 with increments of 0.5 for a total of 21 possibilities. We chose this range of σ values based on the location of the observed Δ distribution (blue histogram in Figure 2) with its approximate peak at about 20 mutual Hill radii. This chosen range of σ values allowed us to explore different distributions of dynamical spacing that spanned a reasonable range of possible model Δ distributions. Increments of 0.5 were chosen as a trade-off between resolution and computational limitations. As will be seen in Section 2.2, the statistically good match between the data and the best-fit model demonstrates that our increments are sufficiently small and have appropriately sampled the possible range of Δ distributions.

Zoom In Zoom Out Reset image size

In order to enforce that adjacent planets are spaced according to the prescribed Δ distribution, we performed the following steps. For each synthetic planetary system, the first planet's orbital period is drawn from the debiased period distribution. If the system's multiplicity is greater than one, for the second planet we draw its separation Δ from the first planet using the prescribed Δ distribution and we also draw a value from the debiased period distribution. If that value is less/greater than the first planet's period, then the second planet will be the inner/outer planet and its exact period will be calculated by satisfying the drawn Δ separation from the first planet. This process repeats if the system has additional planets. These steps are different from Fang & Margot (2012a), where in that study all planetary orbital periods were chosen by drawing them from a debiased distribution. We verified that the periods drawn to match the Δ distribution provide a very close match to the debiased period distribution as well.

After the creation of each model population, we performed synthetic observations of each population's planetary systems by the Kepler spacecraft in order to determine which planets were transiting and detectable (see Fang & Margot 2012a). The transiting requirement was evaluated by picking a random line of sight (i.e., picking a random point on the celestial sphere) and computing the planet–star distance projected on the plane of the sky. The minimum of that distance was compared to the radius of the host star to determine whether or not the planet in our simulations transited. The detection requirement was assessed by calculating each transiting planet's S/N, defined as

$$\begin{eqnarray} &&{\rm S/N} = \left(\frac{R}{R_*}\right)^2 \frac{\sqrt{N}}{\sigma _*}, \end{eqnarray} \tag{ 8 }$$

where the first fraction gives the depth of the transit, N represents the number of transits up to Quarter 6, and σ_* represents stellar noise (Combined Differential Photometric Precision or CDPP; Christiansen et al. 2012). Since CDPP is quarter-to-quarter dependent, we used the median CDPP value over all available quarters. In the calculation of S/N, we took into account gaps between Kepler quarters, the fact that not all stars are observed each quarter, and a 95% duty cycle (Fang & Margot 2012a). If the calculated S/N for a transiting planet met or exceeded the S/N threshold for detection (S/N = 10), then it was considered detectable.

Lastly, we determined the goodness-of-fit between each model population's detected planets and the actual Kepler detections. This was ascertained by comparing the Δ distributions of adjacent planets in their multi-planet systems. We performed a Kolmogorov–Smirnov (K-S) test to assess the fit between the Δ distributions, and this comparison yielded a p-value that we used to evaluate the null hypothesis that the distributions emanate from the same parent distribution. This K-S probability was used to determine how well a particular model matched the observations. We also calculated the goodness-of-fit for multiplicity (by comparing with observed Kepler numbers of transiting systems using a χ² test) and for inclination (by comparing with observed, normalized transit duration ratios using a K-S test) to check that they were consistent with the data (see Fang & Margot 2012a). While we only generated model populations with underlying multiplicity and inclination distributions that are considered to be good fits to the data based on our previous work, this extra step allowed us to confirm that any models with acceptable Δ fits also produced acceptable multiplicity and inclination fits to the Kepler data. We determined which model populations were most consistent with the data by combining (multiplying) the probabilities associated with each one of the three statistical tests that probed multiplicity, inclination, and dynamical spacing. We assumed that these probabilities are independent.

Accounting for all combinations of multiplicity, inclination, and dynamical spacing distributions, in total we created 756 model populations with about 10⁶ planets each. As described earlier, each of these model populations underwent synthetic observations by Kepler as well as statistical tests. The next section presents our results.

We report our results on the dynamical spacing (represented by the criterion Δ) in multi-planet systems based on Kepler data. Using the methods described in the previous section, we determine that our best-fit model for the intrinsic Δ distribution is a shifted Rayleigh distribution (see Equation (7)) with σ = 14.5.

This best-fit distribution is plotted in Figure 2, where we show its probability density distribution (top panel) as well as its cumulative probability distribution (bottom panel). This best-fit distribution has a mean value of Δ = 21.7 with a standard deviation of 9.5. About 50% of neighboring planet pairs have Δ separations larger than 20, and about 90% of neighboring planet pairs have Δ separations larger than 10. This best-fit distribution was obtained by considering shifted Rayleigh distributions with increments in σ of 0.5. The mean values for distributions with σ = 14.0 and σ = 15.0 are Δ = 22.3 and Δ = 21.0, respectively. The combined probabilities for the match to the data using distributions with σ = 14.0 and σ = 15.0 are less than half the combined probability using σ = 14.5.

Our results are valid for the range of stellar and planetary parameters given in Equations (1) and (2), most notably a minimum planet radius of 1.5 R_⊕ and a maximum orbital period of 200 days. It is possible that planets are actually even more closely spaced than this best-fit distribution if there are planets located in intermediate locations with radii less than 1.5 R_⊕. Therefore, our findings about the Δ distribution can be used to represent the spacing of planetary systems confined to the scope of our study, or can serve as an upper limit for the spacing of planetary systems that include planets with smaller radii.

Figure 3 shows the comparison between the Δ distribution of simulated detections from this best-fit model's population and the observed Δ from actual Kepler detections. Note that this figure does not show the underlying Δ distribution that is plotted in Figure 2; instead, this figure shows the distribution of simulated planets that would have been detected, in order to make an appropriate comparison with the observed Δ distribution. The K-S test for comparing these two distributions yields a p-value of 56%, indicating that we cannot reject the null hypothesis that these distributions are drawn from the same parent distribution. In other words, this model is consistent with the observations.

Zoom In Zoom Out Reset image size

Comparison between Figures 2 and 3 shows that the underlying Δ distribution is similar to the observed Δ distribution—both distributions have peaks near Δ ∼ 20. This suggests that the observed Δ distribution is not indicative of a significant population of non-transiting and/or low-S/N planets (within the planet radius and orbital period limits of our study) located inbetween detected planets; otherwise, the underlying Δ distribution would have on average lower Δ values than those of the observed distribution. We caution again that our study cannot rule out the existence of a population of planets with R < 1.5 R_⊕, so the actual underlying Δ distribution could be different from what our results indicate. The similarity between the observed and underlying Δ distributions is due to the rarity of cases where a system has Δ_observed > Δ_underlying (i.e., an undetected planet located in-between two detected planets). The stringent geometric probability of transit means that outer planets are more easily missed than inner planets. As a result, it will be rare to find cases where an intermediate planet is non-transiting and therefore missed, but both the innermost and outermost planets are transiting and detected. For such a case, which rarely occurs, the observed Δ would be greater than the underlying Δ.

We discuss how we expect these results to change if an alternate mass–radius relationship is used. In particular, Figure 1 shows how the observed dynamical spacing distribution changes depending on various choices for the mass–radius relationship. For computational expediency, we chose to evaluate the effects of the mass–radius relationship on the observed Δ distribution. We use this as an approximation for the effects on the underlying Δ distribution, with the justification that the observed and underlying Δ distributions appear similar (see previous paragraph). For the mass–radius relationship shifted down by 1σ, the best-fitting shifted Rayleigh distribution has σ = 16.5. For the mass–radius relationship shifted up by 1σ, the best-fitting shifted Rayleigh distribution has σ = 12.5. As a result, we estimate that our derived dynamical spacing distribution can span σ = 12.5–16.5 due to uncertainties in the mass–radius scaling.