The multiplicity distribution of Kepler's exoplanets

TL;DR

This study infers the true multiplicity distribution of Kepler's exoplanet systems using Bayesian methods, revealing many hidden planets and challenging the need for multiple population models.

Contribution

It demonstrates that a single Zipfian distribution can explain Kepler's observed multiplicities, providing an alternative to dichotomous population models.

Findings

A single Zipfian model fits Kepler data well.

Approximately 2,393 hidden planets are inferred in the studied systems.

About half of single-transiting systems likely host additional planets.

Abstract

The true multiplicity distribution of transiting planet systems is obscured by strong observational biases, leading low-multiplicity systems to be overrepresented in the observed sample. Using the Kepler FGK planet hosts, we employ approximate Bayesian computation to infer the multiplicity distribution by comparing simulated catalogs to the observed one. After comparing a total of ten different multiplicity distributions, half of which were two-population models, to the observed data, we find that a single-population model following a Zipfian distribution is able to explain the Kepler data as well as any of the dichotomous models we test. Our work provides another example of a way to explain the observed Kepler multiplicities without invoking a dichotomous planet population. Using our preferred Zipfian model, we estimate that an additional $2393_{-717}^{+904}$ planets likely reside in the 1537 FGK Kepler systems studied in this work, which would increase the planet count by a factor of $2.22_{-0.36}^{+0.46}$. Of these hidden worlds, $663_{-151}^{+158}$ are expected to reside in ostensibly single-transiting-planet systems, meaning that an additional planet(s) is expected for approximately 1-in-2 such Kepler systems.

Figures

Tables

Table 1. Table 1: Observed multiplicities in the final subset of 1966 KOIs considered in this work. Taking the sum of each multiplicity by its count yields 1966, as expected.

Multiplicity, $m$	Counts, $n_{\mathrm{obs},m}$
1	1225
2	218
3	76
4	15
5	1
6	2
7	0
8	0
9	0
10	0

Table 2. Table 2: One-sigma credibility intervals of the model parameters for each of the ten multiplicity models. Recall that β 𝛽 \beta is defined differently for each model, as described in Section  2.5 .

Model	${\alpha }_{\mathrm{small}}$	${\alpha }_{\mathrm{big}}$	$R_{\mathrm{crit}}$ [$R_{\oplus }$]	${\sigma }_R$ [$R_{\oplus }$]	$\beta$	${\sigma }_I$ [∘]	$f$
Constant	${0.22}_{-0.29}^{+0.42}$	${3.90}_{-1.49}^{+1.78}$	${2.69}_{-0.57}^{+0.61}$	${0.29}_{-0.20}^{+0.29}$	${4.70}_{-0.98}^{+0.54}$	${4.15}_{-0.72}^{+0.68}$	-
Uniform	${2.40}_{-0.49}^{+0.63}$	${2.61}_{-0.98}^{+1.58}$	${2.53}_{-0.37}^{+0.46}$	${0.25}_{-0.16}^{+0.20}$	${0.51}_{-1.77}^{+1.25}$	${5.80}_{-0.49}^{+0.46}$	-
Poisson	$-{0.72}_{-0.13}^{+0.22}$	${4.52}_{-1.63}^{+3.43}$	${3.73}_{-1.00}^{+1.35}$	${1.20}_{-0.64}^{+0.27}$	${5.07}_{-1.35}^{+1.36}$	${0.02}_{-0.02}^{+0.05}$	-
Exponential	${0.53}_{-0.36}^{+0.66}$	${4.55}_{-0.79}^{+1.06}$	${2.49}_{-0.28}^{+0.31}$	${0.19}_{-0.13}^{+0.21}$	${1.65}_{-0.27}^{+0.44}$	${2.62}_{-0.90}^{+0.95}$	-
Zipfian	${0.60}_{-0.36}^{+0.60}$	${4.72}_{-0.78}^{+1.00}$	${2.53}_{-0.29}^{+0.27}$	${0.14}_{-0.09}^{+0.16}$	${0.86}_{-0.29}^{+0.28}$	${2.04}_{-0.69}^{+0.76}$	-
Di-Constant	$-{0.17}_{-0.24}^{+0.23}$	${3.42}_{-0.75}^{+0.98}$	${2.83}_{-0.32}^{+0.33}$	${0.14}_{-0.10}^{+0.20}$	${4.90}_{-0.37}^{+0.40}$	${2.69}_{-0.55}^{+0.33}$	${0.72}_{-0.04}^{+0.04}$
Di-Uniform	${1.88}_{-0.46}^{+0.58}$	${3.94}_{-0.95}^{+1.25}$	${2.55}_{-0.32}^{+0.35}$	${0.19}_{-0.14}^{+0.22}$	${1.98}_{-0.66}^{+0.64}$	${4.03}_{-0.51}^{+0.57}$	${0.65}_{-0.08}^{+0.06}$
Di-Poisson	${0.14}_{-0.28}^{+0.29}$	${5.22}_{-0.75}^{+1.30}$	${2.61}_{-0.29}^{+0.33}$	${0.14}_{-0.09}^{+0.17}$	${5.16}_{-1.18}^{+1.21}$	${0.80}_{-0.33}^{+0.30}$	${0.55}_{-0.05}^{+0.04}$
Di-Exponential	${0.65}_{-0.46}^{+0.69}$	${4.53}_{-0.69}^{+0.95}$	${2.49}_{-0.27}^{+0.30}$	${0.15}_{-0.11}^{+0.19}$	${1.51}_{-0.30}^{+0.40}$	${2.14}_{-0.97}^{+1.11}$	${0.44}_{-0.12}^{+0.11}$
Di-Zipfian	${1.08}_{-0.54}^{+0.86}$	${4.47}_{-0.82}^{+1.14}$	${2.53}_{-0.29}^{+0.33}$	${0.17}_{-0.11}^{+0.18}$	${0.15}_{-0.82}^{+0.57}$	${2.87}_{-0.88}^{+0.95}$	${0.33}_{-0.22}^{+0.23}$

Table 3. Table 3: AIC scores and estimated uncertainties for the ten different models used to describe the Kepler exoplanet multiplicity distribution.

Multiplicity Model	AIC
Constant	$73.75\pm 1.52$
Uniform	$64.17\pm 0.82$
Poisson	$62.72\pm 0.80$
Exponential	$58.96\pm 0.90$
Zipfian	$54.52\pm 0.76$
Di-Constant	$55.92\pm 0.20$
Di-Uniform	$55.02\pm 0.25$
Di-Poisson	$55.87\pm 0.15$
Di-Exponential	$55.39\pm 0.57$
Di-Zipfian	$54.55\pm 0.97$

Table 4. Table 4: AIC scores for twenty fake data sets fitted using two models. Boldened numbers indicate the favored model, which equals the true model in 20/20 cases.

Experiment	AIC (Poisson)	AIC (Zipfian)
Truth = Poisson
1	45.24	56.57
2	45.61	53.93
3	46.56	56.25
4	48.16	61.57
5	51.97	63.72
6	44.78	58.66
7	49.44	54.01
8	44.85	57.28
9	44.58	53.27
10	46.10	52.38
Truth = Zipfian
1	67.99	54.22
2	96.51	58.96
3	68.23	54.99
4	91.57	64.30
5	104.78	60.79
6	61.67	54.25
7	86.42	63.48
8	72.13	59.85
9	85.57	56.84
10	91.90	58.16

Table 5. Table 5: One-sigma credible intervals of the underlying planet multiplicity for Kepler FGK stars with periods 6.25 < ( P / days ) < 400 6.25 𝑃 days 400 6.25<(P/\mathrm{days})<400 and sizes 0.5 < ( R / R ⊕ ) < 32 0.5 𝑅 subscript 𝑅 direct-sum 32 0.5<(R/R_{\oplus})<32 , as computed from the marginal posterior of the favored non-dichotomous Zipfian model. Quoted scores are defined as the percentage of FGK planetary systems with at least one planet in the quoted period and radius range.

Underlying multiplicity	Credible interval
1	${36.53}_{-10.45}^{+10.08}$%
2	${20.23}_{-1.68}^{+1.34}$%
3	${13.52}_{-1.79}^{+1.40}$%
4	${9.59}_{-2.20}^{+2.16}$%
5	${6.83}_{-2.01}^{+2.30}$%
6	${4.80}_{-1.65}^{+2.12}$%
7	${3.27}_{-1.29}^{+1.77}$%
8	${2.17}_{-0.94}^{+1.36}$%
9	${1.45}_{-0.68}^{+1.02}$%
10	${1.02}_{-0.51}^{+0.82}$%

Equations

$$\mathrm{Pr}(\mathrm{detection}|\mathrm{MES})$$

$$\mathrm{Pr}(R)\propto \{\begin{array}{ll}(\log R-\log R_{\min }{)}^{-{\alpha }_{\mathrm{small}}}& \text{if }{R}_{\min }<R\le R_{\mathrm{crit}},\\ (\log R-\log R_{\mathrm{crit}}{)}^{-{\alpha }_{\mathrm{big}}}& \text{if }{R}_{\mathrm{crit}}<R<R_{\max }.\end{array}$$

$$R_H$$

$$\frac{a_2-a_1}{R_H}>{\Delta }_{\mathrm{crit}},$$

$$\mathrm{Pr}(\Delta I)$$

$$\mathrm{Pr}(n_m=n_{\mathrm{obs},m})$$

$$\rho ({\mathbf{n}}_{\mathrm{obs}},{\mathbf{n}}_{\mathrm{sim}})$$

$$\mathrm{Pr}(m|{\beta }_{\mathrm{const}})=\{\begin{array}{ll}1& \text{if }m={\beta }_{\mathrm{const}},\\ 0& \text{otherwise },\end{array}$$

$$\mathrm{Pr}(m|{\beta }_{\mathrm{const}},f)=\{\begin{array}{ll}f& \text{if }m=1,\\ 1-f& \text{if }m={\beta }_{\mathrm{const}},\\ 0& \text{otherwise },\end{array}$$

$$\mathrm{Pr}(m|{\beta }_{\mathrm{uniform}})=\{\begin{array}{ll}\frac{1}{(m_{\max }+1)-{\beta }_{\mathrm{uniform}}}& \text{if }{\beta }_{\mathrm{uniform}}\le m\le m_{\max },\\ 0& \text{otherwise }.\end{array}$$

$$\mathrm{Pr}(m|{\beta }_{\mathrm{poisson}})\propto \{\begin{array}{ll}\frac{{\beta }^m}{m!}& \text{if }1\le m\le m_{\max },\\ 0& \text{otherwise },\end{array}$$

$$\mathrm{Pr}(m|{\beta }_{\exp })\propto \{\begin{array}{ll}{\beta }_{\exp }^m& \text{if }1\le m\le m_{\max },\\ 0& \text{otherwise},\end{array}$$

$$\mathrm{Pr}(m|{\beta }_{\mathrm{zipf}})\propto \{\begin{array}{ll}{m}^{-1-{\beta }_{\mathrm{zipf}}}& \text{if }1\le m\le m_{\max },\\ 0& \text{otherwise},\end{array}$$

$$\mathrm{AIC}$$

Full text

††pubyear: 2019††pagerange: The multiplicity distribution of Kepler’s exoplanets–The multiplicity distribution of Kepler’s exoplanets

1 Introduction

From our vantage point within the Solar System, it seems natural to expect that stars should be accompanied by multiple planets. Around stars similar to the Sun, planetary systems are common, with numerous studies of the Kepler sample converging on an occurrence rate of at minimum one planet per star (Petigura et al., 2013; Foreman-Mackey et al., 2014; Burke et al., 2015; Hsu et al., 2018). In the majority of Kepler systems, there is just a single transiting planet detection (Thompson et al., 2018) but the strong observational biases plaguing transit surveys (Kipping & Sandford, 2016) mean that one might reasonably expect many of these to in fact be multi-planet systems (“multis”) yet to be revealed.

Measuring the multiplicity distribution is crucial, as it is a vital clue to the origins and evolution of detected systems. For example, hot Jupiters rarely reside in multi-planet systems (Wright et al., 2009; Steffen & Agol, 2005; Gibson et al., 2009; Latham et al., 2011; Steffen et al., 2012) barring exceptional cases like WASP-47b (Becker et al., 2015; Weiss et al., 2017). This is often interpreted as evidence for late inward migration from beyond the snow line, leading to scattering of interior planetesimals (Beauge & Nesvorný, 2012; Spalding & Batygin, 2017; Heller, 2018; Dawson & Johnson, 2018). At the other extreme, systems like Kepler-11 (Lissauer et al., 2011a) and TRAPPIST-1 (Gillon et al., 2017) pack half a dozen planets within the orbit of Mercury, which suggests that disk migration and resonant trapping may guide the evolution of such systems (Quillen, 2006; Mustill & Wyatt, 2011; Ormel et al., 2017; Tamayo et al., 2017; Papaloizou et al., 2018).

Unfortunately, the multiplicity distribution is not directly observable from transit surveys like Kepler, because of the extreme biases inherent to the technique. Nevertheless, the unparalleled volume and homogeneous detection biases of the Kepler planets still make it arguably the best resource for the task. If we imagine a system of multiple planets around a star, it is likely that only a subset of the planets (if any) will be detected by a transit survey such as Kepler, because it is guaranteed neither that all of the planets will be aligned with our line of sight, nor that all will transit at high enough signal-to-noise to be detected. Furthermore, Zink et al. (2019) investigate the detection efficiency of the Kepler pipeline and find that it drops in multi-planet systems—specifically, the detection efficiency is higher for the first transiting planet discovered around a star than for subsequent transiting planets in the same system.

Unveiling the true multiplicity distribution from the observed one is therefore a challenging task that needs to account for both geometric and detection biases. Specifically, the mutual inclination distribution between planets should be taken into account in any exploration, since it combines with the underlying multiplicity distribution to produce the observed catalog (Tremaine & Dong, 2012; Brakensiek & Ragozzine, 2016).

One of the first attempts to model the Kepler multiplicity distribution is presented by Lissauer et al. (2011b), who test a Poisson multiplicity model and find that the Kepler catalog is best fit when the mean of this model is equal to 5.5 planets, and the mutual inclinations of the planets are low. However, they find that this best-fitting model significantly underpredicts the number of single-planet systems observed by Kepler. The case for low mutual inclinations in particular has been reproduced in numerous studies (Fang & Margot, 2012; Figueira et al., 2012; Weissbein et al., 2012; Fabrycky et al., 2014). Ballard & Johnson (2016), who study the Kepler M-dwarfs, suggest that the under-prediction of single-planet systems can be resolved by introducing a dichotomous population, with one component being a dynamically cold set of multis and the second being a population of singles or highly-mutually-inclined multis.

This Kepler dichotomy, if it extended to FGK stars, may also explain the under-prediction of singles observed by Lissauer et al. (2011b). This has motivated follow-up efforts to determine if there are fundamental differences in the stellar hosts between singles and multi-planet systems. In particular, Munoz Romero & Kempton (2018) search for, but ultimately find no evidence of, metallicity differences among the hosts of the two types of systems, which might be expected if giant planets were responsible for the dynamically hot population.

The lack of any significant metallicity difference has led some to question the dichotomous hypothesis, even as it explains the observed M-dwarf multiplicities better than a single-population model. The drop in detection efficiency for subsequent planets detected in multi-planet systems noted by Zink et al. (2019) hints at another explanation for the overabundance of observed singles, that the Kepler pipeline simply fails to detect subsequent planets around some percentage of “singles." Indeed, after accounting for this effect, Zink et al. (2019) find that a modified Poisson model fits the observed multiplicities of Kepler GK stars well.

Another alternative explanation to a dichotomous model is a flatter inclination distribution in the inner parts of multi-planet systems than previously assumed—because this inclination distribution works together with the underlying multiplicity distribution to create the observed multiplicities, it is important to model both (Tremaine & Dong, 2012). Bovaird & Lineweaver (2017) show that by adopting a flat-disk model, rather than a “flared" disk, they can match the observed Kepler multiplicities without invoking a dichotomous model.

In this work, we aim to address the question of multiplicity by presenting a comparison of several plausible multiplicity distribution models, including both single and dichotomous populations.

We structure this paper as follows. In Section 2, we introduce the methods, models and inference approach of this work. In Section 3, we present the results, visualizations and analysis of the fits. Finally, we place our work in a broader context in Section 4.

2 Methods

2.1 Input catalog

We downloaded the Kepler DR25 Kepler Objects of Interest (KOIs) catalog via the NASA Exoplanet Archive (NEA; Akeson et al. 2013), with several filters applied. First, we selected only KOIs for which the “Disposition using Kepler Data” was reported as “CANDIDATE”. Second, we required that the NEA-reported surface gravity of the star satisfied $\log g>4$ and that the stellar mass was $0.8<(M_{\star}/M_{\odot})<1.2$, in order to focus on FGK dwarfs. Finally, we filtered for planetary candidates which satisfied $6.25<(P/\mathrm{days})<400$ (the same range considered by Petigura et al. 2013 and Foreman-Mackey et al. 2014) and $0.5<(R_{P}/R_{\oplus})<32$. This led to a population of 1966 KOIs, of which the majority were dispositioned as “CONFIRMED”.

These 1966 KOIs define our observed data set, $\mathcal{D}_{\mathrm{obs}}$, which comprises three key pieces of information. First, the observed multiplicity distribution, which is simply the occurrence tally of multiplicities from 1 to 10 and is reported in Table 1. Second, the list of maximum a-posteriori probability orbital periods, of which there are 1966 elements. Third, the list of maximum a-posteriori probability planetary radii, of which again we have 1966 entries.

After compiling $\mathcal{D}_{\mathrm{obs}}$, we also queried all Kepler target stars for stars which match the filters imposed above. For this, we took the Mathur et al. (2017) DR25 catalog of stellar parameters, which listed 197,096 stars, and cross-matched these with the CDPP6 values (combined differential photometry on a 6-hour timescale; see Christiansen et al. 2012) as obtained from the Mikulski Archive for Space Telescopes (MAST). The mean CDPP6 across all quarters of a given star was saved as the representative CDPP6. In some rare instances, these values were not available on MAST and thus these stars were dropped, leaving us with 196,792 stars. We then applied the same cuts for $\log g$ and $M_{\star}$ as described in the previous paragraph, leaving us with 108,429 FGK dwarfs. This catalog of stars will be used later in Section 2.3.

2.2 Tackling completeness

Our objective is to infer the multiplicity distribution from the Kepler catalog. At a very basic level, this objective is challenged by the incompleteness of the Kepler catalog itself - just because a star has a planet doesn’t mean Kepler is guaranteed (or even likely) to see it. A great deal of attention has been paid to this issue in connection to estimating the underlying planet occurrence rate from Kepler, and so although our objective is distinct, it is useful to briefly review the approaches used in such studies a source of guidance.

The simplest form of incompleteness to deal with is the geometric transit probability, which decreases with increasing planet orbital radius. In estimates of Kepler planet occurrence rates, this can be most easily accounted for by simply dividing apparent occurrence rates by $R_{\star}/a$ (the geometric transit probability), under the assumption of close-to-circular orbits (e.g. see Howard et al. 2012).

The second, and more challenging, component to completeness is detection efficiency. The simplest solution is to limit one’s analysis to a parameter subset where one assumes that the completeness is approximately unity (e.g. Howard et al. 2012; Fang & Margot 2012). This naturally comes at the expense of a smaller sample size. In order to expand the sample to lower signal-to-noise ratio (SNR) events, it is necessary to estimate the detection efficiency in more detail. A typical approach is the so-called “inverse detection efficiency method” (IDEM)111As dubbed by Foreman-Mackey et al. (2014)., where each planet is assigned a detection efficiency score and ultimately the true occurrence rate is inferred by dividing by both the transit probability and the detection efficiency. Detection efficiencies are typically estimated by injection and recovery exercises (e.g. see Petigura et al. 2013; Dressing & Charbonneau 2015). As an example, for a given choice of orbital period and planetary radius, the associated detection efficiency for a particular star may be computed using the KeplerPORTs software (Burke et al., 2015; Burke & Catanzarite, 2017).

In the case of planet occurrence rate estimation, the simplest strategy to account for detection efficiency is the IDEM approach. However, recently Hsu et al. (2018) argue that this approach leads to systemic biases in the inferences since the efficiencies are drawn from estimated planet properties, which are themselves uncertain. Instead, they use a forward-model to inject a population, filter it through a realistic detection efficiency model, and then compare the surviving population to the observed population with some distance metric and Approximate Bayesian Computation (ABC). This approach is shown to more faithfully infer the occurrence rate of the injected population and so a similar approach is adopted here. Rather than comparing occurrences, our work ultimately is interested in the frequency of various multiplicities, but the same approach can be employed (and is discussed further in Section 2.4).

In principle, it should be possible to define a detection efficiency model unique to each star using KeplerPORTs. However, when conducting Bayesian inference, detailed calculations of these efficiencies for every star and at every step in period and radius comes at high computational cost. A simpler yet still accurate approach is to use a global Kepler detection efficiency model, for which one inputs the so-called multiple event statistic (MES) of a planetary candidate and the model returns a detection probability. This is appropriate since we primarily care about the ensemble rather than individual systems. One example of a global detection efficiency model comes from Christiansen et al. (2016), who use the actual Kepler detection pipeline to inject and recover and planets and find that, for FGK stars, the average detection efficiency is well-approximated by either a cumulative gamma function or a logistic function. We use the latter in this work since it is faster to compute. It is given by

[TABLE]

where $b_{l}$, $c_{l}$ and $d_{l}$ are coefficients defined in Christiansen et al. (2016) using Kepler DR24. An update to the cumulative detection fraction versus MES is presented in Thompson et al. (2018) for DR25, who find a similar distribution which we use in this work.

The MES (see Jenkins 2002) is a statistic that measures the combined significance of all the observed transits in the detrended, whitened Kepler light curves, assuming a linear ephemeris. In practice, it is not feasible to generate very large populations of synthetic planets (required for Bayesian Monte Carlo work), inject their transits, detrend, whiten, and thus compute MES in the same way as the real Kepler pipeline.

Instead, we use the transit SNR as a proxy for the MES in what follows. We note that the two are not equal—specifically, the MES depends on the goodness-of-fit between a transit search template and a transit signal mediated by stellar noise, while the SNR does not depend on the template—but they are, to first order, proportional to each other (Burke & Catanzarite, 2017). By combining the SNR with Equation (1), we are able to estimate the detectability of any synthetic KOI.

2.3 The forward model

Our model works by first choosing a random star from the filtered stellar catalog described in Section 2.1. We then inject a planetary system around it composed of $m$ planets, where $m$ is always less than or equal to $m_{\mathrm{max}}=10$ and is drawn from a chosen multiplicity distribution as described later in Section 2.5. Each of the planets is assigned a random period drawn from a log-uniform distribution from 6.25 to 400 days. A log-uniform distribution was chosen since it both provides a reasonably close match to the observed marginalized period distribution reported by Foreman-Mackey et al. (2014), and is the same assumption used in previous multiplicity studies, such as Ballard & Johnson (2016).

Next, the innermost planet in the system is assigned a random radius drawn from a double-sided power law (DSPL) distribution, described by

[TABLE]

We normalize the DSPL distribution such that the two sides meet at $R_{\mathrm{crit}}$ and integrate to unity over the interval $R_{\mathrm{min}}<R<R_{\mathrm{max}}$. The terms $R_{\mathrm{min}}$ and $R_{\mathrm{max}}$ are fixed to 0.5 $R_{\oplus}$ and 32 $R_{\oplus}$ respectively, but the parameters $R_{\mathrm{crit}}$, $\alpha_{\mathrm{small}}$ and $\alpha_{\mathrm{big}}$ are treated as unknown shape parameters to be inferred.

To reflect the observed “peas-in-a-pod" covariance of planet radii (Weiss et al., 2018), the radii of subsequent planets in the system are drawn from a Gaussian distribution centered at the innermost planet’s radius. The scale parameter of this distribution, $\sigma_{R}$, is treated as another free parameter to be inferred. This parameter is able to extend out to very large values, thereby accounting for the possibility of no correlation (Zhu, 2019).

The simplified double-sided power law radius distribution is designed to capture the turn-over in planet occurrence seen at around mini-Neptune radii, reported in numerous studies (Fressin et al., 2013; Petigura et al., 2013; Foreman-Mackey et al., 2014). It does not, however, describe the radius valley reported by Fulton et al. (2017). This effect was only revealed by substantial improvements to the precision of measured stellar radii, and we argue that it is not influential enough to significantly affect our study which focusses on multiplicity.

The radius distribution used here essentially represents a set of nuisance parameters which is marginalized over in the final results.

Having generated $m$ proposal planets around the star, with periods and radii drawn from the distributions described above, we next check whether the system is dynamically stable. Following the same approach as Ballard & Johnson (2016), we test for Hill stability using Equation (3) of Fabrycky et al. (2014). Specifically, we define the mutual Hill radius between planets “1” and “2” as

[TABLE]

where the Hill stability criterion is satisfied if

[TABLE]

where $M$ and $a$ refer to the masses and semi-major axes of the planets. To estimate $M$ for each planet, we use the maximum a-posteriori probability forecaster mass-radius relation derived by Chen & Kipping (2017). We compute semi-major axes are from periods using the stellar mass and Kepler’s Third Law.

The critical separation is $\Delta_{\mathrm{crit}}=2\sqrt{3}$ for neighboring planets, and for three-or-more planets, Fabrycky et al. (2014) require $\Delta_{\mathrm{inner}}+\Delta_{\mathrm{outer}}>18$ for neighboring inner and outer pairs of planets.

If the proposed planetary system violates Hill stability, we use the same star and same multiplicity but make a new realization of the periods and radii for the planetary system. We allow this process to repeat up to 1000 times, after which we abandon the star and draw a new star from the KIC catalog.

After this point we have generated a stable multi-planet system. Next, we need to calculate how many of these planets actually transit. The innermost planet, labelled with subscript “1”, has a transit impact parameter, $b_{1}$, of $(a_{1}/R_{\star})\cos I_{1}$, where $I_{1}$ is the orbital inclination angle222Note that a near-circular orbit is assumed here and throughout.. Inclination is isotropically distributed and thus we adopt a uniform distribution for $\cos I_{1}$, which in practice means that we draw a random real number for $\cos I_{1}$ from $\mathcal{U}[0,1]$, where $\mathcal{U}$ denotes a uniform distribution. The other planets in the system are assumed to have inclinations perturbed away from this angle by an angle $\Delta I$, representing their mutual inclinations within a flared disk (see Bovaird & Lineweaver 2017 for a flat-disk model). $\Delta I$ is drawn from a Rayleigh distribution characterized by a scale parameter $\sigma_{I}$:

[TABLE]

For each planet in the system, we calculate the impact parameter $b_{j}=(a_{j}/R_{\star})\cos{(I_{1}+\Delta I_{j})}$. Any planet for which $b_{j}<1+(R_{j}/R_{\star})$ is treated as a transiting planet and is saved. Systems with zero transiting planets need not be considered further and are discarded, leading us to draw a new star from the KIC catalog.

At this point, we now have a simulated system of at least one transiting planet orbiting a chosen KIC star. The final component of our forward model is to simulate what fraction of transiting planets in the system would actually be detectable. To do this, we first assign each planet a random transit epoch. Next, we query which quarters that particular KIC star was observed by Kepler for, since many stars were not observed in every quarter due to spacecraft rotation and loss of CCDs during the mission.

Using our simulated ephemeris for each planet, we can now calculate how many transits of each planet would have been observed by Kepler. We estimate the SNR of each planet using Equation (10) of Kipping & Sandford (2016), multiplied by the square-root of the number of observed transits. Finally, the detection probability is computed using Equation (1) from Christiansen et al. (2016). To decide if the transiting planet is detectable or not, we make a random Bernoulli draw, with probability equal to this computed detection probability.

This process culminates in a set of $m$ simulated detected transiting planets around a particular star. At this point, we loop back to the beginning of the forward model and keep going until 1966 detected planets have been generated (since this represents the size of the observed sample, $\mathcal{D}_{\mathrm{obs}}$, that we will ultimately compare to).

The forward model therefore ultimately yields a simulated data set, $\mathcal{D}_{\mathrm{sim}}$, with the same elements and form as $\mathcal{D}_{\mathrm{obs}}$. Further,$\mathcal{D}_{\mathrm{sim}}$ is clearly dependent upon the simulation’s choice of multiplicity, radius and inclination distributions - which are characterized by model parameters $\bm{\theta}=\{\beta,\alpha_{\mathrm{small}},\alpha_{\mathrm{big}},R_{\mathrm{crit}},\sigma_{R},\sigma_{I}\}$ (where $\beta$ is a stand-in term(s) describing the multiplicity distribution, described below in Section 2.5, and the other terms have been previously defined).

We note that since non-detections are discarded and not counted, our approach does not enable an estimate of the underlying planet occurrence rate.

2.4 Comparison to observations

A single run of the forward model described in Section 2.3 generates a population of simulated detected transiting planets described by $\mathcal{D}_{\mathrm{sim}}$. Our task is now to infer the parameters of the forward model, $\bm{\theta}$, which would have generated the observed Kepler systems, $\mathcal{D}_{\mathrm{obs}}$, by comparing $\mathcal{D}_{\mathrm{sim}}$ to $\mathcal{D}_{\mathrm{obs}}$. In particular, we are interested in inferring the parameters of the multiplicity distribution, $\beta$.

Conventional Bayesian inference might proceed using hierarchical Bayesian modeling (HBM), where the multiplicity distribution is described by some parameterized form and then each system’s true multiplicity is treated as a free parameter drawn from this overall distribution - giving rise to a large number of unknown variables to solve for (see Hogg et al. 2010 for an astronomer’s introduction to HBMs). In this case, the likelihood function used for inference would be well-defined as the product of the likelihoods for each individual system.

Hierarchical models allow for rigorous inference but typically come at great computational expense. Instead, we seek to learn the multiplicity distribution by comparing some distance metric which quantifies how closely the simulated population matches the observed population - thereby ignoring the individual systems and treating the population as an ensemble. By using one or more distance metrics to quantify goodness-of-fit, we are thus conducting what is typically referred to as approximate Bayesian computation, or ABC (see Ishida et al. 2015; Hahn et al. 2017; Hsu et al. 2018; Witzel et al. 2018 for recent applications in astronomy).

The three key ingredients for ABC inference are a forward model which generates $\mathcal{D}_{\mathrm{sim}}$, prior distributions for the model parameters, $\mathrm{Pr}(\bm{\theta})$, and a distance function $\rho(\mathcal{D}_{\mathrm{obs}},\mathcal{D}_{\mathrm{sim}})$ which quantifies how well the simulated distribution resembles the observed sample.

Although our primary goal is to learn the multiplicity distribution, we elect to define a distance metric which considers the agreement between the simulated and observed multiplicities but also the agreement between the simulated and observed radius distribution. This is because the two cannot be assumed to be independent: planetary radii determine planetary masses, which in turn determine their stability and whether they could reside in a high-multiplicity system. We therefore infer not only the multiplicity model parameters, but also the radius distribution parameters. Further, the mutual inclination distribution strongly influences the fraction of planets observed to transit and thus is also a parameter we should expect to constrain and be covariant with the other model terms.

Our goal here, of course, is not to infer the true radius distribution of the Kepler catalog, nor the distribution of mutual inclinations among its multi-planet systems. However, the known inter-relationships between these terms necessitates that we have some reasonable description of them and that we freely explore them in conjunction with the multiplicity distribution. At the end, we can simply marginalize over these “nuisance” terms in our final calculation of the multiplicity distribution.

Having established that we require a distance metric which incorporates both the multiplicities and radii, let us consider the multiplicity component first. Previous works have most commonly invoked a Poisson likelihood function in comparing a simulated multiplicity to the observed value (Weissbein et al., 2012; Ballard & Johnson, 2016). This essentially asserts that probability distribution for the observed number of $m$-planet systems, $n_{m}$, is a Poisson distribution with a mean rate given by $n_{\mathrm{sim},m}$. The Poisson model is well-motivated for inference based on counting statistics, such as this, and thus is adopted in this work too. Accordingly, the probability of observing a particular number of $m$-planet systems, $n_{\mathrm{obs},m}$, is given by

[TABLE]

The Poisson likelihood function is defined by a product of the above over all $m$ (=multiplicities), and this function certainly describes how close a simulated set of multiplicities, $\mathbf{n}_{\mathrm{sim}}$, resembles the observed set, $\mathbf{n}_{\mathrm{obs}}$ - thereby providing a suitable distance metric:

[TABLE]

where $\mathbf{n}_{\mathrm{obs}}$ and $\mathbf{n}_{\mathrm{sim}}$ represent vectors containing all the $m$-indexed observed and simulated population multiplicities, respectively. The $\mathbf{n}_{\mathrm{obs}}$ vector is fixed and given by Table 1. Meanwhile, the $\mathbf{n}_{\mathrm{sim}}$ vector, which needs to be counted up after each forward model call, is directly controlled by the choice of forward model parameters $\bm{\theta}$, which we ultimately wish to infer. It should be noted that our choice of distance metric here, shown in Equation (7), does not decrease as the distributions approach one another, but rather increases. Accordingly, it could perhaps be better thought of as an inverse distance metric although we’ll continue to refer to it as distance metric in what follows, with the only important consequence being that our task is to maximize $\rho$, rather than minimize $\rho$.

We now turn our attention to the component of the distance metric which characterizes the planetary radius distribution. The objective here is not to fit each and every planetary radius - which would be more in line with an HBM. Instead, we wish to simulate a population whose statistical properties broadly match those of the observations. A straightforward approach for accomplishing this is to the use the Kolmogorov-Smirnov (K-S) test, following on from the approach adopted by Fang & Margot (2012).

We therefore compute the K-S $p$-value between the observed radii and the simulated set as our radius distance metric, since this follows the behaviour of the multiplicity component in terms of being a term we seek to maximize. We multiply this by the multiplicity distance metric given by Equation (7) to define an overall distance metric, $\rho(\mathcal{D}_{\mathrm{obs}},\mathcal{D}_{\mathrm{sim}})$. The two components are equally weighted under this definition.

A variety of sampling techniques are suitable for ABC inference (Beaumont, 2019), and in this work we elect to use the Markov Chain Monte Carlo (MCMC) approach (Marjoram et al., 2003; Marin et al., 2012). We sample the model parameter space of $\bm{\theta}$ with Gaussian proposals where the acceptance criterion is chosen such that the probability of accepting a proposal is $\rho_{\mathrm{proposal}}/\rho_{i}$, where $i$ denotes the current index in the chain. In this way, improvements in the distance metric (which recall equates to an increase in the “distance” under our definition) are always accepted. This means that samples near the beginning of the chain, prior to convergence, can often have poor distance scores. We remove these burn-in samples by only including samples in the chain past the first instance exceeding the median distance metric. We demand that 50,000 accepted samples are achieved for each model, with the final chains inspected to verify convergence and mixing.

Because of the somewhat subjective nature of choosing an appropriate distance function in ABC, there is no formal guarantee the model posterior will converge to the true posterior distribution. For this reason, it is important to test the accuracy of our ABC inference framework through fake data generation and recovery simulations, which we present in Section 3.2.

2.5 Proposed multiplicity models

In Section 2.3, we described how exoplanetary radii could be described using a DSPL distribution, mutual inclinations with a Rayleigh distribution, and orbital periods with a log-uniform distribution. However, we did not propose a specific form for the multiplicity distribution itself - which we turn our attention to here.

Specifically, we here describe ten choices of multiplicity distribution: five single-population models, each parameterized by a single free parameter $\beta$, and five corresponding “dichotomous" models. The dichotomous models are mixture models of two populations: (1) a fraction $f$ of single-planet systems and (2) a fraction $1-f$ of multi-planet systems, distributed according to one of the models parametrized by $\beta$.

2.5.1 Constant model

One of the first multiplicity models proposed in the literature is presented by Ballard & Johnson (2016), who initially adopt a simple approach where every system has the same multiplicity, $\beta_{\mathrm{const}}$, which is treated as a free parameter (we refer to this as the constant model):

[TABLE]

Ballard & Johnson (2016) conclude that this model is unable to provide a good fit to the observed multiplicities of Kepler M-dwarf systems and thus expand upon the constant distribution in that same work to include a second component of single-planet systems ($m=1$), which represent a fraction $f$ of all systems:

[TABLE]

Following the terminology used by the authors, we refer to this as the “dichotomous” constant model.

Simple models are often attractive since one might plausibly purport that the laws that govern planetary architectures (or whatever other phenomenon one is considering) are fundamentally simple themselves. However, we suggest here that the constant model is almost certainly too simple a model to be a realistic description of the exoplanet multiplicity distribution. It is rather implausible to suppose that every system should have an identical number of planets, and even after including a second population of singletons, this still yields a highly unnatural distribution composed of two distinct peaks at $m=1$ and $m=\beta_{\mathrm{const}}$, with zero probability that systems have multiplicities other than this. For this reason, we felt motivated to consider other models in addition to the constant formalism.

2.5.2 Uniform model

A simple improvement to consider would be to adopt a discrete uniform distribution, where every multiplicity is just as likely as any other, above some minimum multiplicity, $\beta_{\mathrm{uniform}}$:

[TABLE]

This can be similarly be extended to a dichotomous model by assuming that some fraction of planetary systems, $f$, belong to a separate population of singles.

2.5.3 Truncated Poisson model

Arguably, a more natural model is a Poisson multiplicity distribution, which appears to be the most commonly adopted law (e.g. see Lissauer et al. 2011b; Fang & Margot 2012; Gaidos et al. 2016; Bovaird & Lineweaver 2017). This might be expected if the multiplicity of an exoplanetary system were the result of a constant rate of generating planets within a fixed interval of time or space. One may write that the probability of forming an $m$-planet system would thus be

[TABLE]

where we drop normalization terms which do not depend on $m$. For $\beta_{\mathrm{poisson}}\geq 2$, this implies a peaked, non-monotonic distribution at a specific multiplicity, unlike the uniform case. In practice, we reject any trial $m$ equal to zero or exceeding $m_{\mathrm{max}}$, i.e. we truncate the distribution. The Poisson can again be extended to a dichotomous model as was done before.

2.5.4 Exponential model

Another previously adopted law is that of a discrete exponential distribution (e.g. Bovaird & Lineweaver 2017) which imposes that the multiplicity, $m$, follows

[TABLE]

which can again be extended to a dichotomous model as above.

2.5.5 Zipfian model

Finally, we consider a Zipfian distribution, which represents a discrete power-law given by

[TABLE]

and a corresponding dichotomous model constructed as above.

Zipf’s Law (Zipf, 1935) is known to be an excellent approximation for word frequency versus rank for human languages, and even some animal communications (Doyle et al., 2011). On this basis, it might seem like a peculiar choice to use when modeling the exoplanet multiplicity distribution, but Zipf’s Law also appears in much wider array of problems, such as the population frequency of cities (Auerbach, 1913). Zipf’s Law has been argued to be a natural by-product of models with many underlying latent variables, somewhat analogous to the arguments behind the Central Limit Theorem (Belevitch, 1959; Aitchison et al., 2016), and thus on this basis would seem a very reasonable model to propose for exoplanets too - despite the fact it has seemingly not been used in the past for this purpose.

3 Analysis

We run the forward model described in Section 2.3 ten times, once for each choice of multiplicity model described above, and fit for the free parameters $\bm{\theta}=\{\beta,\alpha_{\mathrm{small}},\alpha_{\mathrm{big}},R_{\mathrm{crit}},\sigma_{R},\sigma_{I}\}$ (and $f$, for the dichotomous models) via ABC, as described in Section 2.4.

In these fits, we adopt a uniform prior on $R_{\mathrm{crit}}$ between $R_{\mathrm{min}}=0.5R_{\oplus}$ and $R_{\mathrm{max}}=32R_{\oplus}$, and a uniform prior on $f$ between 0 and 1. We also adopt (improper) priors insisting that $\alpha_{\mathrm{small}}$ and $\alpha_{\mathrm{big}}$ be greater than $-1$ and that $\sigma_{R}$ and $\sigma_{I}$ be positive.

One-sigma credibility intervals for these parameters in each of our ten model fits are presented in Table 2, and an example posterior distribution, for the single-population Zipfian model, is presented in Figure 1.

3.1 Model comparison

To compare models, we use the Akaike Information Criterion, AIC (Akaike, 1974). There are two major reasons behind this choice. First, the AIC does not require that one of the models being tested is the correct model, it merely asks which of the models is the closest approximation to the truth (unlike the BIC; Schwarz 1978). Second, the AIC does not functionally depend on the sample size, which is somewhat ill-defined in our problem since our inference employed a likelihood approximation. We therefore calculate the AIC for each model by first finding the most probable realization from the 50,000 posterior samples, as defined by the distance function, and then using

[TABLE]

where $k$ is the number of free parameters used by each model. Since we constructed our distance metric $\rho$ as a product of two likelihood-like terms, $\rho$ approximates the likelihood in the AIC calculation above. For the non-mixture models, $k=6$ since we have free parameters $\beta$, $\sigma_{I}$, $R_{\mathrm{crit}}$, $\sigma_{R}$, $\alpha_{\mathrm{small}}$ and $\alpha_{\mathrm{big}}$. The dichotomous models add one extra free parameter, the fraction of single-planet systems $f$.

We assign uncertainties to our AIC scores through a bootstrapping procedure. First, we split the chain up into $S$ segments. For each segment, we compute the AIC, and estimate its standard deviation as 1.4826 multiplied by the median absolute deviation. We repeat this procedure, varying $S$ from 2 to 50 in unity steps, and for the non-mixture models use the median score across all experiments. For the dichotomous models, we find that the scatter tends to decrease as we approach small $S$, and thus we fit a simple quadratic model of scatter versus $S$ to estimate the scatter at $S=1$. These uncertainty estimates, along with the overall AIC scores, are compiled in Table 3.

Amongst the non-mixture models, Table 3 shows that the Zipfian distribution is preferred, favored over the next-best model (exponential) with an odds ratio of $e^{(58.96-54.52)/2}=9.21$. Since all of these models have the same number of parameters, this preference is purely driven by the much improved distance metrics. The constant model is found to be the worst description of the multiplicity distribution, disfavored versus a Zipfian model by a factor of 15,000.

Amongst the dichotomous models, the field is much more level, with all five models roughly equally favorable to each other and also to the single-population Zipfian distribution model. As a check on this, we also tried computing the Savage-Dickey ratio by evaluating the posterior density at $f_{\mathrm{single}}=0$. Only the di-exponential and di-Zipfian models had enough samples around this region to reliably estimate the single-population model:dichotomous model odds ratio, yielding ratios of 0.055 and 1.12, respectively. These are broadly consistent with the AIC results of approximately equal weights ($=e^{-\Delta\mathrm{AIC}/2}$), demonstrating that the AIC approach is a suitable approximation for this model selection problem.

On this basis, we conclude that the simpler hypothesis of a single population model is not significantly rejected by the current data. A single Zipfian distribution appears quite capable of describing the Kepler exoplanet multiplicity distribution for FGK hosts.

In Figures 2 and 3, we plot the underlying simulated population of planetary systems (inset figures) and the “detected" subset of these systems, for all ten multiplicity models, to investigate the effect of detection biases on this subset and to compare it to the real Kepler detections. Despite very different underlying multiplicity models, the detected sample is qualitatively similar in all ten cases: strongly peaked at $m=1$, and falling off at higher multiplicities. It is not surprising in this light that the Zipfian model performs best of the five single-population models, as it has this general shape already, and that the dichotomous models (which by definition include a peak at $m=1$) perform equally well.

Consequently, we conclude overall that the current Kepler data prefer multiplicity models which peak at $m=1$, but have little distinguishing power between such models. Choosing among them thus becomes a question of prior beliefs about the underlying planet distribution—e.g., is there theoretical support for two planetary system formation pathways?

In what follows, we investigate more fully the single-population Zipfian multiplicity model, on the basis of its simplicity.

3.2 Testing the inference framework

The use of ABC and also AIC model selection are both approximate tools and thus one might reasonably question how robust they really are. To test this, we decided to generate a total of 20 fake data sets where the true multiplicity distribution is known and test how well we can recover that true distribution using the same machinery used thus far.

We generate the first ten mock populations assuming a Poisson multiplicity distribution. Every population has the same input parameters, chosen to be close to the inferred solutions in our earlier fits, specifically $\alpha_{\mathrm{small}}=0.33$, $\alpha_{\mathrm{big}}=5.0$, $R_{\mathrm{crit}}=2.5$ $R_{\oplus}$, $\sigma_{R}=0.05$ $R_{\oplus}$, $\beta_{\mathrm{Poisson}}=12.0$, $\sigma_{I}=2.0^{\circ}$. However, these ten fake “observed” data sets, $\mathcal{D}_{\mathrm{obs}}^{\prime}$, are slightly different to each other due to the stochastic nature of the forward simulation. The second ten are generated in the same way except we switch to a Zipfian distribution (replacing $\beta_{\mathrm{Poisson}}$ with $\beta_{\mathrm{Zipf}}=1.0$).

We fit each of these twenty data sets with two models: a Poisson and a Zipfian. Thus, we should be able to test whether AIC scoring is able to pick out the correct model in each case - a basic assumption upon which the previous subsection rests. Second, we can test whether the inferred parameters (in cases where the correct model is regressed) are compatible with the input values. In this way, we can provide a detailed assessment of the validity of our inference framework. To save computational time, we use 10,000 post-burn-in steps for each MCMC fit.

The AIC results, summarized in Table 4, show that the correct model is identified in 20 out of 20 cases. In general, the Zipfian model appears to be more flexible and gets closer to describing the Poisson model than vice versa, likely as a result of the very harsh selection functions which push the distributions towards an ostensibly monotonic form. Nevertheless, the AIC scoring system appears to be a reliable tool for identifying the best model.

Comparing the actual parameters which result, as shown in Figure 4, we find good agreement between the results and the injected truths: over the 10 fits to the 6 parameters, the recovered parameter is in $1\sigma$ agreement with the injected parameter in 60 of 60 cases for the Poisson trial and 54 of 60 cases for the Zipf trial. We also inspected the distances metrics versus parameter samples and verified that the distance metrics approach their maximum around the true injected values, as expected. This establishes that the ABC inference framework is able to accurately recover the correct parameters, as well as being suitable for model selection via the AIC.

3.3 Properties of the preferred Zipfian distribution

Given that the non-dichotomous Zipfian multiplicity model is the favored model of this work, it is worthwhile to consider the parameters inferred from this model. We show a corner plot of the joint posteriors in Figure 1, where the converged, unique nature of the inferred solution is evident.

It is also instructive to compare the multiplicity of the generated systems in our fits, versus the apparent multiplicity of these same systems after being filtered through our mock Kepler pipeline. This is is shown in Figures 2 and 3, where one can see how the true sample (inset figures) is considerably diminished as a result of detection bias.

3.3.1 Interpreting the Zipfian slope

When interpreting the inferred value of $\beta_{\mathrm{Zipf}}$, it is worth highlighting that $\beta_{\mathrm{Zipf}}=-1$ leads to a precisely uniform distribution, $\beta_{\mathrm{Zipf}}<-1$ leads to distributions whose probability density monotonically increases with increasing multiplicity, and vice versa $\beta_{\mathrm{Zipf}}>-1$ leads to distributions which decrease with increasing multiplicity.

With a shape parameter of $\beta_{\mathrm{zipf}}=(0.86^{+0.28}_{-0.29})$, our fit strongly favors a distribution which decreases with increasing multiplicity. This is evident from Table 5, which presents the relative frequency of each multiplicity as determined from the Zipfian fit. Since Mercury is less than half an Earth radius and Mars is beyond the period threshold used in this work, the Solar System is a 2-planet system in our framework - a configuration found in $(20.2^{+1.3}_{-1.7})$% of our simulated FGK planetary systems. Packed, compact systems are rare, with 6-planets or more constituting $\sim 13\%$ of our simulated FGK planetary systems.

3.3.2 Probability of additional planets in known systems

One can also see that although the fraction of observed one-planet systems represents 80% of all systems (see Table 1), in reality only 37% of system are truly single (see Table 5). Another way to think about this is that the Zipfian model finds that 37% of the simulated systems are genuinely single and all of these must yield a planet with the correct geometry and detectability to have been “detected” by the simulated Kepler survey (else they would not have been included in the final simulated catalog since our code would have not saved the realization). Since 100% of the 37% truly single planet systems appear as singletons in the final catalog, $80-37=43\%$ of the detected planets are multiple planet systems for which only one planet was detected to transit. Thus, of the 80% of ostensibly single planet systems, 37/80=46% are indeed genuinely single and the other 43/80=54% are yet-to-be-revealed multi-planet systems.

Accordingly, radial velocity follow-up of single transiting FGK Kepler systems has an a-priori $54.1_{-12.3}^{+12.9}$% chance of detecting new planets (after correctly propagating the uncertainties) with periods and radii in our specified range. Given that there are 1225 single-planet systems in our sample, that equates to $663_{-151}^{+158}$ hidden planets in the single-planet systems.

3.3.3 Total number of missing planets

By calculating the total number of planets generated in the simulated systems, we find that the Zipfian model predicts a total of $4359_{-717}^{+904}$ planets residing around the 1537 FGK systems with known detections. Since only 1966 known planets reside around these stars, that means that there are $2393_{-717}^{+904}$ hidden planets - which are expected to be dynamically stable. Discovering these planets, perhaps through radial velocity follow-up, could increase the planet count around these stars by a factor of $2.22_{-0.36}^{+0.46}$.

4 Discussion

The principal finding of this work is that the observed multiplicities of the Kepler FGK transiting systems can be well-explained without invoking a dichotomous population model. Specifically, we find that a Rayleigh mutual inclination distribution with a Zipfian multiplicity distribution (the latter of which appears to have never been tried before) is able to well-reproduce the observed catalog. This is not to say that dichotomous models are disfavored—indeed, the single-population Zipfian and the five dichotomous models perform equally well—only that invoking a dichotomous population is not necessary to explain the detected Kepler multiplicities. Furthermore, we find that the Kepler data do decisively prefer multiplicity models peaked at $m=1$ over those peaked at higher multiplicities.

Bovaird & Lineweaver (2017) also suggest that the dichotomous model may not be necessary by considering an alternative inclination distribution. Since inclination and multiplicity both affect the final catalog (Tremaine & Dong, 2012), then it is certainly plausible then either (or both) of these effects are able to explain the observed multiplicities without invoking dichotomy. Zink et al. (2019) note, furthermore, that the Kepler pipeline’s decreased detection efficiency for multi-planet systems could also explain the overabundance of Kepler singles. Finally, although our work centers on FGK stars, we highlight that Gaidos et al. (2016) also find that a dichotomous distribution may not be necessary by changing the underlying models in the case of M-dwarfs.

It is curious that Zipf’s Law (Zipf, 1935), most commonly associated with linguistics, works well for exoplanet multiplicities. Zipfian laws are argued by Aitchison et al. (2016) to be natural outcomes of systems involving a large number of latent variables, and this may represent another example. Extending our analysis to M-dwarfs, particularly from TESS, will provide a good test as to whether the Zipfian model can persist in the face of new data.

Using our preferred model, we are able to make predictions about the numbers of missing planets. For example, we predict that 7 or more planet systems are rare, with just 7.9% of detected systems being so packed. This is in sharp contrast to Mulders et al. (2011), who recently estimated that 42% of Sun-like stars have nearly coplanar planetary systems with 7 or more exoplanets. Although our numbers are not measuring precisely the same quantity, it would be difficult to reconcile the Mulders et al. (2011) value with our estimates given the stark paucity of such systems in our observed sample.

Our model does predict a large number of missing planets, $\simeq 2400$ around the 1537 host stars considered, of which some $\simeq 660$ reside in ostensibly single-transiting-planet systems. It may therefore be possible to test the predictions of these models by conducting radial velocity follow-up of the Kepler field in the future to measure the true multiplicities.

Acknowledgments

ES, DK, & MC acknowledge support from the Columbia University Data Science Institute “Seed Funds Program”. Thanks to members of the Cool Worlds Lab for useful discussions in preparing this manuscript. We are grateful to the anonymous reviewer for their constructive feedback, and to Jessi Cisewski-Kehe for guidance on ABC.