THE EXOPLANET CENSUS: A GENERAL METHOD APPLIED TO KEPLER

The probability that a planet on a randomly oriented circular orbit of semimajor axis a transits its host star, with radius R_* and mean density ρ_*, is

$$\begin{eqnarray} \eta _{\rm tr} &=& {R_\ast \over a} = \left(3 \pi \over G \rho _\ast P^{2} \right)^{1/3}\nonumber \\ &=& 0.051 \left(10 \; {\rm days}\over P \right)^{2/3} \left(\rho _\odot \over \rho _\ast \right)^{1/3} \end{eqnarray} \tag{ 1 }$$

for R ≪ R_*. Eccentricity tends to increase the detectability of planets at fixed a (Burke 2008).¹ The duration of Kepler transits imply a mean eccentricity ≲ 0.1–0.25 (Moorhead et al. 2011). The ≲ 10% effect on η_tr is ignored here.

For simplicity, we approximate the mean stellar density (or more specifically the mean ρ^−1/3_*) as solar. The densities of planet hosts can be estimated from precise light curve parameters, though corrections for the eccentricity apply (Tingley et al. 2011) and the required precision is not universally achieved. Furthermore, planet hosts have a lower density than the target star population, on average. This bias arises from the higher transit probability around lower density stars. Uncertainties in the average stellar density modestly and only weakly affect our results due to the weak (−1/3) power-law dependence and the restriction to solar-type stars.

Transits with large impact parameters, 0.9 ≲ b ⩽ 1, are detected relatively rarely due to their short duration and unusual shape. This selection bias can be included as either a modification to η_tr or η_disc, but of course not both. Observed impact parameters can be used to modify the transit efficiency to η_tr = b_maxR_*/a for an appropriate b_max < 1 (Catanzarite & Shao 2011). Alternately, the dependence of η_disc on the transit impact parameter can be precisely modeled.

One might (incorrectly) expect multi-planet systems to require special values of η_tr. When mutual orbital inclinations are small, finding one planet does increase the odds that others will be found. However, since the planetary systems as a whole are randomly oriented, low mutual inclinations also make it easier to miss an entire planetary system. The mutual inclination of planets within systems has no effect on the overall detection rate (aside from sampling noise, as always).

Precisely quantifying the discovery efficiency of the Kepler survey is difficult and ongoing work. Our study relies on HKep11 estimates of discovery efficiency for a subsample of bright solar-type stars with relatively high η_disc. The sample consists of N_* = 58, 041 stars that satisfy the effective temperature, surface gravity, and Kepler magnitude cuts:

$$\begin{equation} 4100\; {\rm K}\le T_{\rm eff} \le 6100\; {\rm K} \end{equation} \tag{ 2a }$$

$$\begin{equation} 4.0 \le \log g /({\rm cm} \,{\rm s}^{-2}) \le 4.9 \end{equation} \tag{ 2b }$$

$$\begin{equation} {{K}_P} \le 15 {\rm \,mag}. \end{equation} \tag{ 2c }$$

Section 2.3 discusses the planet candidates in this sample.

The discovery efficiency of transiting planets with a signal-to-noise ratio (S/N) >10 on a log-uniform grid covering 0.68 < P/(days) < 50 and 1 < R/R_⊕ < 16 is reported in Figure 4 of HKep11. HKep11 provide η_discN_* over this range in the bottom left of each cell. The resulting values of η_disc are plotted in Figure 2.

Zoom In Zoom Out Reset image size

For bins with no detected planets, HKep11 do not report a discovery efficiency. This omission is not because planet detections are required to estimate efficiencies. Rather HKep11 applied efficiencies only to bins with detections, which differs from our approach of applying the efficiencies across all parameter space and to unbinned data. The missing efficiencies in Figure 4 of HKep11 can be readily obtained by interpolation.

Instead of merely interpolating the reported Kepler efficiencies, we are motivated by their smooth variation to find an analytic fit. Our best fits are for power laws, not to η_disc itself, but to the related

$$\begin{equation} r_{\rm disc}\equiv {\eta _{\rm disc}\over 1 - \eta _{\rm disc}} \,, \end{equation} \tag{ 3 }$$

the ratio of discoverable to non-discoverable planets. Since r_disc ranges from zero to arbitrarily large values, a power-law fit will automatically satisfy the η_disc ⩽ 1 restriction.

Figure 2 shows that the broken power law is an excellent fit to reported discovery efficiencies. The power-law fit is

$$\begin{equation} r_{\rm disc} = 0.324 \left(R \over R_\oplus \right)^{6.34} \left(P \over {\rm day} \right)^{-1.07} \end{equation} \tag{ 4 }$$

below a break at

$$\begin{equation} R < R_{\rm b} = 2.2 \left(P \over {\rm day}\right)^{0.17} R_\oplus \,. \end{equation} \tag{ 5 }$$

The efficiency at the break is quite high, η_disc = 0.98 (and remarkably constant, probably reflecting the high S/N threshold). Note that η_disc drops to 50% (r_disc = 1 in Figure 2) at R = 0.54R_b, rising roughly from 1 to 2 R_⊕ over the periods considered. When discovery efficiencies are small, η_disc ≈ r_disc ≪ 1, and the power law in Equation (4) applies directly to η_disc.

While the scalings for η_disc > 0.98 are of little practical concern for our study, we report the fit above the break for completeness:

$$\begin{equation} r_{\rm disc}(R > R_{\rm b}) = 5.46\left(R \over R_\oplus \right)^{2.74} \left(P \over {\rm day} \right)^{-0.460} \,. \end{equation} \tag{ 6 }$$

This fit describes a small fraction of noisy stars, but does not affect our results because η_disc ≈ 1 in this regime.

The relevant fit in Equation (4) agrees well with theoretical estimates. For S/N-limited detections, the combined efficiency should scale as

$$\begin{equation} \eta _{\rm disc}\eta _{\rm tr} \propto P^{-5/3} R^6; \end{equation} \tag{ 7 }$$

see, e.g., Equation (7) of Gaudi (2007). The small R limit of Equation (4) gives η_discη_tr∝P^−1.74R^6.34, rather good agreement.

The efficiencies summarized here only apply to planets detected during Q2 of the Kepler mission (BKep11). Estimates will need to be updated for future data releases. Moreover, future estimates will likely be improved by more detailed studies that include extraction of simulated transits from the actual Kepler data analysis pipeline. Accurate estimates of the Kepler discovery efficiency for all relevant parameters—especially stellar T_eff—are the most crucial ingredient for determining the frequency of Earth-like planets.

Tables 1 and 2 of BKep11 list the properties of 1235 planet candidates and their host stars. For the η_disc estimate to be valid, we can only consider the subset of planets that satisfy the conditions set by HKep11, that (a) host stars are within the solar sample set by Equation (2), (b) P < 50 days, and (c) the transit detection has S/N > 10. The last condition is discussed further below. Applying these filters reduces the number of planet candidates from 1235 to 567, for a detected planet-to-star ratio of just under 1%.

Our "full" planet sample considers

$$\begin{equation} 0.5 \; {\rm days}< P < 50 \; {\rm days} \end{equation} \tag{ 8a }$$

$$\begin{equation} 0.5 \,R_\oplus < R < 20 \,R_\oplus \,. \end{equation} \tag{ 8b }$$

Adding these radius and the minimum period cuts only removes five planets (from 567 to 562). Instead of using the BKep11 reported values for R, we compute R from the reported values for the host star radius, R_*, and R/R_*.²

We now discuss relevant differences between our planet sample and that of HKep11. HKep11 conservatively restricted attention to R > 2 R_⊕, because the discovery efficiency is high for these planets (except at longer periods). We extend our analysis down to R = 0.5 R_⊕. This extension includes all small planets in the solar sample, the smallest having R = 0.58 R_⊕. This extension involves (a) trusting the HKep11 efficiencies down to the 1 R_⊕ level to which they were reported and (b) extrapolating the efficiencies down to 0.5 R_⊕. There are 19 planets in the sample with 0.5 < R/R_⊕ < 1.0. Both the smooth variation of the reported efficiencies and their agreement with theory (shown above) increase our confidence in making the extrapolation. We also report results for samples of larger planets that are unaffected by this potentially risky extrapolation.

HKep11 used different S/N values than the publicly reported BKep11 S/N values on which we rely. HKep11 define S/N based on only one quarter of data (Q2, Kepler's first full quarter). By contrast, BKep11 report S/N values up to Q5, so that roughly four times more data have been collected. Our goal is to translate the BKep11 S/N values to the Q2 S/N values for which the η_disc estimate is valid. However, this translation is not so simple.

With perfect noise scaling, the BKep11 S/N values should be roughly twice as high as those used by HKep11. Table 2 of HKep11 presents noise scalings for four planets with ∼2.5 R_⊕ radii and ∼40 day periods. The ideal scaling roughly holds for three objects, but the S/N is saturated and only rose modestly for the fourth. It is unclear what the general noise scaling is and how it varies with size and toward shorter periods.

We proceed cautiously by performing our analyses on both an S/N >10 and an S/N >20 planet sample (as defined by BKep 11). These cases cover the complete range of possibilities between saturated noise and perfect scaling. Fortunately, the adopted S/N threshold has a modest effect on results presented here.

We briefly note a discrepancy between the planet counts reported by BKep11 and HKep11. HKep11 report a sample of 438 planet candidates (of the 1235 released by BKep11) that (a) have hosts that satisfy the stellar parameter cuts of Equation (2), (b) have planet parameters P ⩽ 50 days and 2 R_⊕ ⩽ R ⩽ 32 R_⊕, and (c) satisfy the HKep11 definition of S/N > 10. If we apply cuts (a) and (b) only to the BKep11 tables, but allow all S/N, there are 378 candidates. That number increases to 393 using the reported R instead of calculating it as (R/R_*) × R_*. Thus, there are at least 60 (or 45) planets that HKep11 say should appear in the BKep11 tables, but they do not.

Applying any S/N threshold to the BKep11 tables can only make the discrepancy larger. We made no choices that would cause us to miss planets in the BKep11 tables. Cuts were applied inclusively (including any values that fall on a boundary) and all the planets in any system were counted. We conclude that a discrepancy exists between the planet and/or stellar parameters appearing in the current versions of HKep11 and BKep11. Speculatively, parameters may have been refined over time. This mismatch does not affect our results as long as the estimates of the discovery efficiencies in HKep11 are sufficiently accurate.

For a survey of N_* stars, we divide planet detections into bins indexed by $\ell$, with $N_{\ell}$ planet detections in a given bin. For generality the meaning of the bin is here kept arbitrary. If the average detection efficiency per bin is $\eta_{\ell}$, then the best estimate of the NPPS in each bin is

$$\begin{equation} f_\ell = {N_\ell \over \eta _\ell N_\ast }\,, \end{equation} \tag{ 9 }$$

the number of detected planets per star divided by the efficiency.

The total planet occurrence, $\sum_\ell} f_{\ell}$, depends on bin size when $\eta_{\ell}$ is not constant. For arbitrarily small bins, there is only one planet per bin, and the efficiency is evaluated for the parameters of each planet and its host. This small bin limit is the only truly assumption-free way to estimate the planet fraction and suggests that binning of data is never required.

A maximum likelihood analysis using Poisson statistics can reproduce the estimates of $f_{\ell}$ given by Equation (9) and develop intuition for the general method. The number of expected detections in all bins is

$$\begin{equation} N_{\rm exp}= N_\ast \sum _\ell \eta _\ell f_\ell \,. \end{equation} \tag{ 10 }$$

The likelihood of the detections is a product of the Poisson distributions for each bin,

$$\begin{equation} \tilde{L} = \left[ \prod _\ell {(N_\ast \eta _\ell f_\ell)^{N_\ell } \over N_\ell !} \right] \exp (-N_{\rm exp}) \,. \end{equation} \tag{ 11 }$$

We are interested in the values of $f_{\ell}$ that maximize the likelihood (and its logarithm). As a consequence, we can ignore any constant multiplicative factors, where constant means independent of the $f_{\ell}$ values. This simplification results in a likelihood function,

$$\begin{equation} L = \left[ \prod _\ell f_\ell ^{N_\ell } \right] \exp(-N_{\rm exp}) \,. \end{equation} \tag{ 12 }$$

The efficiencies now appear only in N_exp, a feature that has greater significance for the parameterized fits discussed below.

The maximum likelihood values of $f_{\ell}$ are the roots of ∂ln (L)/∂$f_{\ell}$ = 0, which reproduce the expected result of Equation (9).

We now describe a general method to fit a parameterized PLDF to unbinned data. We consider a PLDF that depends on planet properties, ${{\bm x}}$ (an N_x component vector) and stellar properties ${{\bm z}}$ (with N_z components). Usually, the elements of these vectors are logarithms of measured quantities such as the orbital period or stellar mass.

The probability that a star with properties ${{\bm z}}$ has a planet with properties ${{\bm x}}$ that lie in a volume³ $d {{\bm x}}$ of phase space is

$$\begin{equation} df = {\partial f ({{\bm x}},{{\bm z}}) \over \partial {{\bm x}}^{N_x}} d{{\bm x}}\,. \end{equation} \tag{ 13 }$$

Defined this way, the integral $f({{\bm z}}) = \int df$ represents the NPPS with stellar properties ${{\bm z}}$ over the range of planet properties, ${{\bm x}}$, covered by the integral.

We express the differential distribution,

$$\begin{equation} {\partial f ({{\bm x}},{{\bm z}}) \over \partial {{\bm x}}^{N_x}} \equiv C g({{\bm x}},{{\bm z}}; {{\bm \alpha}})\,, \end{equation} \tag{ 14 }$$

as an amplitude, C, times a shape function, g. A functional form for g must be assumed, and this function can include any number of shape parameters, ${{\bm \alpha}}$. The main result of the maximum likelihood analysis is the best-fit values of ${{\bm \alpha }}$. As a simple example, the values of ${{\bm \alpha}}$ can be power-law exponents.

The total number of planets around N_* stars (indexed by j) is

$$\begin{equation} \hspace{-3.8pc}N_{\rm tot} = C \sum _j^{N_\ast } \int g({{\bm x}},{{\bm z}}_j;{{\bm \alpha}}) d{{\bm x}} \end{equation} \tag{ 15a }$$

$$\begin{equation} = N_\ast C\int \mathcal {F}_\ast ({{\bm z}}) g({{\bm x}},{{\bm z}};{{\bm \alpha}}) d{{\bm x}}d {{\bm z}}, \end{equation} \tag{ 15b }$$

where this expectation value applies over the range of ${{\bm x}}$ covered by the integral.⁴ In Equation (15b), the sum over individual stars is replaced by an integral over the stellar distribution function $\mathcal {F}_\ast ({{\bm z}})$, which is normalized to one as $\int \mathcal {F}_\ast ({{\bm z}}) d{{\bm z}}= 1$. This replacement approximates individual stellar properties as obeying a continuous distribution. Such an approximation is not required. Indeed direct summation over all stars is always possible (and preferred when feasible). We proceed with the integral notation to minimize indices.

Consider a survey with a net detection efficiency $\eta ({{\bm x}},{{\bm z}})$. As described for the case of transit surveys in Section 2, the net efficiency is the product of (a) any number of uncorrelated detection probabilities (all ⩽ 1) that arise from individual selection effects and (b) a false-positive efficiency (⩾1) that gives the ratio of all detections (planets and false positives) to actual planets.

The expected number of planet detections (including false positives when relevant) is then

$$\begin{equation} N_{\rm exp}= N_\ast C F, \end{equation} \tag{ 16 }$$

where the shape integral,

$$\begin{equation} F \equiv \int \eta ({{\bm x}},{{\bm z}}) \mathcal {F}_\ast ({{\bm z}}) g({{\bm x}},{{\bm z}};{{\bm \alpha}}) d{{\bm x}}d{{\bm z}}, \end{equation} \tag{ 17 }$$

weights the shape function over both the stellar distribution and the efficiency.

3.2.1. Maximizing Likelihood

Consider a survey that detects N_pl planets (some of which could be false positives) around N_* stars. We define the likelihood of the data as

$$\begin{equation} \tilde{L}= \left[ \prod _{i=1}^{N_{\rm pl}} df_i \right] \exp (-N_{\rm exp}) \,. \end{equation} \tag{ 18 }$$

Analogous to Equation (12), this likelihood represents the probability of each individual detection, labeled by i, times the exponential probability that no other planets were detected.

Applying Equations (13) and (14), we eliminate unnecessary constants, including the phase space volume. The likelihood function becomes

$$\begin{eqnarray} L &=& \left[C^{N_{\rm pl}} \prod _{i=1}^{N_{\rm pl}} g({{\bm x}}_i, {{\bm z}}_{j(i)};{{\bm \alpha}})\right] \exp (-N_{\rm exp}) \,, \end{eqnarray} \tag{ 19 }$$

where j(i) refers to the star that hosts planet detection i. Letting g_i represent the shape function evaluated at detection i, the log likelihood is

$$\begin{equation} \ln (L) = N_{\rm pl}\ln (C) + \sum _{i = 1}^{N_{\rm pl}} \ln (g_i) - N_{\rm exp}\,. \end{equation} \tag{ 20 }$$

The best-fit normalization, C, is found by maximizing the likelihood as the root of ∂ln L/∂C = 0:

$$\begin{equation} C = {N_{\rm pl}\over N_\ast F} \end{equation} \tag{ 21 }$$

using Equation (16).

Following TT02, we use this constraint to eliminate C from the likelihood:

$$\begin{eqnarray} \ln (L) &=& -N_{\rm pl}\ln (F) + \sum _i^{N_{\rm pl}} \ln (g_i) + \left\lbrace N_{\rm pl}\left[\ln \left(N_{\rm pl}\over N_\ast \right)-1 \right]\right\rbrace,\nonumber\\ \end{eqnarray} \tag{ 22 }$$

where the constant term in curly brackets can again be ignored.

To find the best-fit values of ${{\bm \alpha}}$, follow from the roots of $\partial \ln (L) /\partial {{\bm \alpha}}= 0$ as

$$\begin{equation} {\partial \ln F \over \partial {{\bm \alpha}}} = { 1 \over N_{\rm pl}} \sum _i^{N_{\rm pl}}{ \partial \ln (g_i) \over \partial {{\bm \alpha}}}\,. \end{equation} \tag{ 23 }$$

These N_α (the number of parameters in ${{\bm \alpha}}$) equations generally couple to each other and require numerical solution, but see the Appendix. In Equation (23), the detection efficiencies affect the left-hand side (via F), while detections determine the right-hand side.

Next, we summarize the basic steps in obtaining a best-fit solution. First, make a hypothesis for the functional form of the PLDF and express it in the form of Equation (14). Second, solve Equation (23) for the shape parameters. Third, using the best-fit ${{\bm \alpha}}$, calculate the normalization via Equations (21) and (17). The best-fit solution is now complete and can be used (e.g.) to estimate the NPPS as N_tot/N_* using Equation (15). Fourth, estimate the errors on the fit from the contours of the likelihood function in Equation (20), as explained in Section 3.2.3 below.

3.2.2. Interpretation

This formalism has a straightforward interpretation. The normalization C matches the numbers of detected and expected planets, N_pl = N_exp, as seen by comparing Equations (16) and (21).

The shape parameters similarly match the expected and detected averages of planet and host star observables. Consider the case of power-law distributions with a shape function

$$\begin{equation} g = \exp ({{\bm \alpha}}_{\rm pl}\cdot {{\bm x}}+ {{\bm \alpha}}_\ast \cdot {{\bm z}})\,. \end{equation} \tag{ 24 }$$

Here ${{\bm x}}$ and ${{\bm z}}$ are the logarithms of quantities being fit to a power law, and ${{\bm \alpha}}= \lbrace {{\bm \alpha}}_{\rm pl}, {{\bm \alpha}}_\ast \rbrace$ distinguishes the power-law exponents for planetary and stellar parameters. With this power-law shape function, Equation (23) gives

$$\begin{equation} \langle {{\bm x}}\rangle _F = \langle {{\bm x}}\rangle _{\rm obs}, \langle {{\bm z}}\rangle _F = \langle {{\bm z}}\rangle _{\rm obs}, \end{equation} \tag{ 25 }$$

where

$$\begin{equation} \langle {{\bm x}}\rangle _F \equiv {1\over F} \int {{\bm x}}\eta \mathcal {F}_\ast g d{{\bm x}}d{{\bm z}} \end{equation} \tag{ 26a }$$

$$\begin{equation} \langle {{\bm x}}\rangle _{\rm obs} \equiv {1 \over N_{\rm pl}} \sum _i^{N_{\rm pl}} {{\bm x}}_i \,, \end{equation} \tag{ 26b }$$

and similarly for the stellar parameters. Non-power-law shape functions correspond to a more complex averaging of the observables.

We now comment on the validity of treating planet occurrence as a Poisson process, i.e., as random, uncorrelated events. This approach may sound inconsistent with interactions in multi-planet systems, but it is not. In effect, the method investigates whether a random drawing from a PLDF (weighted by selection biases) provides a good fit to the observations. The possibility that planet interactions give rise to a complicated PLDF can be included and does not affect the validity of this procedure. Furthermore, the vast majority of planet pairs (and possible planet pairs as well) involve different stars which should be completely non-interacting. Thus, the ability of planet interactions in individual systems to affect the PLDF of all planets in all systems is quite weak.

This discussion emphasizes one of the main features of the maximum likelihood approach. For better and worse, the functional form of the PLDF must be prescribed. The appropriateness of the assumed family of functions must be evaluated by some external method of hypothesis testing.

3.2.3. Errors

To quantify the uncertainty in the shape parameters, we compare contours of the likelihood L to a Gaussian distribution. The maximum likelihood, L_max, is given by Equation (22) evaluated at best-fit ${{\bm \alpha}}$. The 1σ error ellipse covers ${{\bm \alpha}}$ values that satisfy ln (L) ⩾ ln (L_max) − 1/2. The general nσ contours follow ln (L) = ln (L_max) − n²/2. The uncertainty in the normalization C is only the range of values taken by C within the error ellipse of the ${{\bm \alpha}}$.

For one-dimensional errors on an individual shape parameter, we hold other parameters fixed at their best-fit values. With strong covariance, and skewed error ellipses, such one-dimensional errors would not be very meaningful. As we will show (in Figure 6), such covariance is mild in our analysis.

We emphasize that a maximum likelihood analysis does not evaluate the quality of a fit. The errors on ${{\bm \alpha}}$ measure the fit's precision. For instance, a large planet sample will precisely define an average power law even if the actual PLDF deviates strongly from a power law.

Our analysis ignores uncertainties in the planetary and stellar parameters. Here, we describe a simple way to include these errors, mainly to show that this is possible. A more rigorous approach would use Fisher information matrices, but is not developed here.

Uncertainties associated with an individual detection, i, affect the value of the shape function, g_i, that appears in the likelihood analysis. To include these uncertainties, g_i should be a weighted integral over the uncertainties,

$$\begin{equation} g_i \equiv \int \phi _i({{\bm x}}^{\prime } - {{\bm x}}_i,{{\bm z}}^{\prime } - {{\bm z}}_{j(i)}) g({{\bm x}}^{\prime },{{\bm z}}^{\prime };{{\bm \alpha}}) d{{\bm x}}^{\prime }d{{\bm z}}^{\prime }, \end{equation} \tag{ 27 }$$

where ϕ_i is an appropriately normalized Gaussian distribution (for instance) centered on the best-fit values of the planet detection and its host star. Ignoring uncertainties is equivalent to setting the error distributions, ϕ_i, to δ-functions.

Uncertainties in the target star properties affect the number of expected planets via the shape integral F. Errors can be included similarly to Equation (27) as

$$\begin{equation} \hspace{-0.7pc}F ={1 \over N_\ast }\sum _j^{N_\ast } \int \phi _j({{\bm z}}^{\prime } - {{\bm z}}_{j}) \eta ({{\bm x}},{{\bm z}}^{\prime })g({{\bm x}},{{\bm z}}^{\prime }) d{{\bm x}}d{{\bm z}}^{\prime } \end{equation} \tag{ 28a }$$

$$\begin{equation} \approx \int \phi ({{\bm z}}^{\prime } - {{\bm z}})\mathcal {F}_\ast ({{\bm z}}) \eta ({{\bm x}},{{\bm z}}^{\prime })g({{\bm x}},{{\bm z}}^{\prime }) d{{\bm x}}d{{\bm z}}^{\prime } d{{\bm z}}. \end{equation} \tag{ 28b }$$

The two forms above show how errors could be applied to each star j or (more approximately) how average errors could be assigned to the stellar distribution function. The integration over both ${{\bm z}}$ and ${{\bm z}}^{\prime }$ in Equation (28b) covers both the distribution of average stellar properties ${{\bm z}}$ (replacing the sum over individual stars in Equation (28b)) and the range of uncertainties ${{\bm z}}^{\prime }$ allowed by the error distributions ϕ. Estimates of the total number of planets N_tot of Equation (15) are modified by a similar inclusion of the error distribution ϕ.

Systematic uncertainties in the detection efficiencies are more difficult to quantify. Assigning detection efficiencies is the most dangerous step in the statistical analysis of survey data. It is safest to perform multiple analyses that cover a range of (hopefully reasonable and nearly complete) possibilities for detection efficiencies. That approach is the one taken here.

Any analysis must choose to study some subset of planetary and stellar parameters and to ignore others. The freedom to safely ignore a parameter in the PLDF only exists if detection efficiencies remain constant with respect to the ignored parameter. A detection efficiency that stays near 100% obviously qualifies as constant. As always, these conditions apply only to the region of parameter space being studied.

To understand the inherent danger in ignoring parameters, note that the PLDF is assumed to be averaged over all ignored parameters. However, an accurate prediction of the yield of a survey requires weighting the PLDF by the detection efficiencies before integrating over parameter space, as in Equation (16). Thus, ignoring a parameter in the PLDF makes it impossible to correct for selection biases related to that parameter.

A safer—though admittedly more difficult—way to ignore a parameter (that affects the detection efficiency) is to first include that parameter in the fit and then marginalize over it. With this caveat made explicit, we now show how the general formalism treats the case of a PLDF that depends only on planetary parameters, ${{\bm x}}$, or only on stellar parameters, ${{\bm z}}$, (most famously stellar metallicity).

3.3.1. Planetary Parameters Only

For a shape function, $g({{\bm x}};{{\bm \alpha}})$, that is independent of stellar properties, the shape integral simplifies to

$$\begin{equation} F = \int \eta _\ast ({{\bm x}}) g({{\bm x}};{{\bm \alpha}}) d{{\bm x}}\,, \end{equation} \tag{ 29 }$$

where the stellar-averaged efficiencies are

$$\begin{equation} \eta _\ast ({{\bm x}}) = \int \eta ({{\bm x}}, {{\bm z}}) \mathcal {F}_\ast ({{\bm z}}) d{{\bm z}}. \end{equation} \tag{ 30 }$$

The fit procedure is simplified in this case. The stellar averaging is done a single time and not for every choice of ${{\bm \alpha}}$ during optimization.

This case is the one most relevant to this work. Specifically in Section 4 we consider a shape function

$$\begin{equation} g = \exp (\alpha x + \beta y) = \left(R \over R_o\right)^{\alpha }\left(P \over P_o\right)^\beta \,. \end{equation} \tag{ 31 }$$

This power-law distribution in planetary radius and orbital period identifies the planetary parameters

$$\begin{equation} {{\bm x}}= \lbrace x,y\rbrace = \left\lbrace \ln (R/R_o), \ln (P/P_o)\right\rbrace \end{equation} \tag{ 32 }$$

as logarithms in radius and period, relative to reference values R_o and P_o, respectively. The shape parameters ${{\bm \alpha}}= \lbrace \alpha, \beta \rbrace$ are power laws. The Appendix gives an analytic solution for these best-fit power-law exponents for an ideal transit survey with unity discovery efficiency and no false positives.

3.3.2. Stellar Parameters Only

We now consider fitting a PLDF that depends on stellar parameters only. As cautioned above, this simplification would be a bad idea for transit surveys since detection efficiencies always vary with orbital period (due to the geometric transit probability) and often vary with planetary radius (for smaller planets that give shallow transits). We consider this case for other survey types and for completeness.

The primary conceptual difficulty is that the PLDF was defined in Equation (14) as a differential distribution in one or more planet parameters, ${{\bm x}}$. However, we can remove all references to ${{\bm x}}$ by integrating over the observed volume $V = \int d{{\bm x}}$ of parameter space. The non-differential PLDF,

$$\begin{equation} f({{\bm z}}) = C G({{\bm z}}; {{\bm \alpha}}_\ast) \,, \end{equation} \tag{ 33 }$$

gives the NPPS with characteristics ${{\bm z}}$. The integration of the shape function gives G ≡ gV.

The expected planet yield is still N_exp = N_*CF, but the shape integral simplifies to

$$\begin{equation} F = \int \eta ({{\bm z}}) \mathcal {F}_\ast ({{\bm z}}) G ({{\bm z}},\alpha _\ast)d{{\bm z}}\,. \end{equation} \tag{ 34 }$$

The best-fit solutions are still given by Equations (21) and (23), with ln g_i replaced by ln G_i.

This derivation demonstrates how the method could be applied to study the abundance of planets versus stellar metallicity and/or mass alone. We again caution that this approach assumes that the detection efficiencies do not vary with the properties of the detected planets.

Our most complete sample of planets covers 0.5 < P/day < 50 and 0.5 < R/R_⊕ < 20 (see Section 2.3 for details). Figure 3 plots planet counts of this full sample. The planet counts are binned by R (or P) in the left (right) panel. Raw counts are divided by the number of stars surveyed and the logarithmic bin size.⁵ Normalized this way, the data measure the PLDF as weighted by the detection efficiencies, η,

$$\begin{equation} \left\langle {\eta {d f \over d \ln R}}\right\rangle {\rm and} \left\langle {\eta {d f \over d \ln P}}\right\rangle \,. \end{equation} \tag{ 36 }$$

Brackets indicate marginalization over P and R, respectively. Aside from sampling noise, these values are independent of the survey size or choice of bin width.

Zoom In Zoom Out Reset image size

The maximum likelihood fits do not use this binned data. Rather, the joint power-law fits use the technique of Section 3.2 (specifically the special case of Section 3.3.1) to analyze the unbinned data. These fits use the detection efficiencies to estimate the true PLDF. For a visual comparison with the planet counts, we project the fits back into observational space in Figure 3. This projection requires weighting the PLDF by the detection efficiencies and marginalizing over P (or R) to give the quantities in Equation (36).

We now explain in detail how a power-law PLDF is modified by this projection. Though the best-fit period distribution in the right panel of Figure 3 has β = 1.3, the curve in observational space has β_obs ≃ 0.3, smaller by about one. A decrease of precisely 2/3 is due to the transit efficiency. The remainder (here ≃ −1/3) is due to the period dependence of η_disc. From Equation (4), this non-geometric correction could be as large as ≈ − 1, if all detections were S/N limited. Thus, β_obs could be lower than the underlying β by as much as −5/3. It is often assumed that the geometric correction of −2/3 adequately relate the underlying and observed power laws in P, but this approximation is insufficient.

Selection effects not only shift the slope of a power law, they can also introduce curvature. Since Figure 3 is log–log, any deviation of the power-law fits from a straight line is entirely due to selection effects. Though not easily visible, there is modest curvature to the projected period law in the right panel due to lower discovery efficiencies at longer periods.

The curvature in the fit to the size distribution is much more dramatic, as shown in the left panel of Figure 3. The sharp bend coincides remarkably well with the peak in planet counts below ∼3 R_⊕. The agreement is so good that we must emphasize that the break is controlled by the best estimates of η_disc, with no adjustable free parameters. Above ∼3 R_⊕, the Kepler discovery efficiency is unity for all periods considered. Thus, the slope of the projected size distribution at large radii matches the best-fit α. The break starts at R ≲ 3 R_⊕, because η_disc drops below 90% for these sizes, as shown in Figure 2. At small R, the power-law slope of the projected size distribution is α + 6.34, as fixed by the radial dependence of Equation (4).

Thus, the drop in planet counts toward small radii is very well modeled by selection effects. This agreement is strong evidence that the underlying planet counts continue to rise toward small sizes, down to the ∼0.5 R_⊕ limit of the observations.

The uncertainty of the fits are shown in Figure 3 by the thick red curve covering 1σ deviations from maximum likelihood. Figure 6 shows the corresponding error ellipse in the α and β parameters (the small black oval is for the full planet sample). The 1σ errors on α and β reported in Table 1 are one-dimensional errors holding the other parameters fixed. As explained in Section 3, the normalization C is not varied during the fit procedure, but follows from the values of α, β, and the planet-to-star ratio. The errors on C indicate the range of values taken inside the 1σ error ellipse of α and β.

Our best-fit size distribution has α = −1.73 ± 0.07. This value indicates that smaller planets are much more abundant, but that larger planets likely contain most of the mass. At constant planet density, α < −3 would give small planets most of the mass. However, larger planets tend to have lower density, ρ. Take ρ∝R⁻¹ as an example.⁶ Then α < −2 would give small planets more of the mass (over the range where the density law holds). Though the mass–radius relation remains uncertain, we tentatively conclude that larger planets dominate the mass distribution of Kepler planets.

The best-fit period law, β = 1.33 ± 0.03, indicates that planets are more closely packed further from the star, for logarithmic intervals in P or semimajor axis a. For β > 3/2, the planet density per linear interval in a would increase.

However, Figure 3 shows that a single power law is a qualitatively poor fit to the data due to a sharp drop in planet counts below P ∼ 3 days. After a broad maximum, planet counts drop slightly for periods above ∼10 days, but (as explained above) the drop is so shallow that the underlying period distribution continues to rise. The obvious break in the period distribution motivates our division of the planet sample below.

Since the period distribution deviates strongly from a simple power law, we divide the planets into short and long period samples. We use P = 7 days as the dividing line. The choice is only partly motivated by the proximity to the break in the period distribution, which is gradual and more prominent at shorter periods. We also want each sample to have significant numbers of planets and an adequate range of periods over which to measure a power-law slope. Planet counts for different data cuts are summarized in the N_pl column of Table 1.

Figure 4 shows independent power-law fits to the short and long period data, again projected against the planet counts. Two power laws vastly improve the qualitative fit to the period distribution (right panel). The large exponent at short periods β_fast ≈ 2.2 reflects the sharp drop in planet counts to short periods (despite the flattening of the planet counts for P > 3 days). Though planet counts decline to longer periods, the long period exponent β_slow ≈ 0.6 shows that the underlying distribution continues to rise. Longer period planets are needed to determine where the distribution turns over. It is quite exciting that we have not yet found where most of the planets lie.

Zoom In Zoom Out Reset image size

The size distributions in the left panel shows overlapping histograms for the short and long period samples. The longer period ("slow") sample has a steeper radius distribution than the shorter period ("fast") with α_slow − α_fast ≈ 0.5 ± 0.2. By eye, the difference in the distribution appears driven by the flatness of the fast sample between ∼3 and 15 R_⊕. To clarify this behavior, we now examine the differences between small and large radius planet samples.

We subdivide the planets into four samples by splitting both the shorter and longer period populations (still defined relative to 7 days) into "small" and "big" samples relative to 3 R_⊕. The results of the four independent power-law fits are plotted against planet counts in Figure 5. Figure 6 plots error ellipses to compare the slope changes and their statistical significance. The error ellipses of the subsamples are much larger than the full sample due to lower numbers and the reduced range in R and P (over which to define a slope).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The period distributions are only modestly affected by this split. At both short and long periods, the best-fit β of the small and big samples agree within the statistical uncertainties. Nevertheless, some qualitatively interesting features appear in the period distributions shown in the right panel of Figure 5. The small planet counts begin their decline below about 7 days (as noted by HKep11). The large planet sample remains flat toward shorter periods—with evidence of a peak at P ≈ 3 days. This peak is responsible for the flat appearance of the overall period distribution down to 3 days, as seen in the right panels of Figures 3 and 4.

The behavior of the size distributions is even more intriguing. For big planets, the best-fit radius law changes by α_{big, slow} − α_{big, fast} = −1.31 ± 0.47, with almost 3σ significance. Since the discovery efficiency of Kepler is near unity for these large planets, no obvious systematic uncertainties should be at work. For small planets, the change in the size distributions, α_{small, slow} − α_{small, fast} = 0.70 ± 0.47. While of more modest significance, this change is of the opposite sign.

Figure 7 plots the best-fit size distributions for each quadrant, with the full planet sample shown for comparison. At longer periods the size distribution flattens toward smaller radii, i.e., α is less negative for the small planets. For the shortest period planets the behavior is precisely opposite. The size distribution is quite shallow at large sizes and steepens for the small radius sample.

Zoom In Zoom Out Reset image size

The implications of this behavior are intriguing. The processes that form and/or migrate planets inside ∼7 days clearly disfavor planets of a few to several R_⊕. Atmospheric evaporation can at least qualitatively explain these trends. At short periods, it seems that only the most massive giants—true Jovians with R ≳ 11 R_⊕—can efficiently retain their atmospheres. Lesser giants are more readily stripped of their atmospheres, including ices that sublimate from the surface. The demotion of ice and gas giants to small rocky cores can explain the flattening of the large planet size distribution at short periods.

While the data are strongly suggestive, the evaporation hypothesis does have complications. First, even at periods below 7 days, lower gravity ice giants might retain their atmospheres (Ehrenreich & Désert 2011). Also, evaporation can only partly explain the sharp rise in small planets at short periods. Evaporation could provide some small cores, but not nearly enough. An additional mechanism for the preferential parking of smaller versus larger planets appears to be required. Nevertheless, it appears likely that atmospheric retention plays a key role in the observed size distributions, providing another pillar of support for the core–accretion hypothesis.

To investigate how systematic uncertainties in the discovery efficiency might affect our results, we experimented with different S/N thresholds. As discussed in Section 2.3, S/N thresholds of 10 and 20 should provide a reasonable bound on the uncertainties.

The results so far have described the more inclusive S/N >10 sample. Figure 3 includes dotted histograms that correspond to the higher S/N >20 sample. Planet counts are noticeably lower below 3 R_⊕ and across all periods.

Table 1 includes fits with both S/N thresholds. All of the qualitative trends discussed in this section persist in the higher S/N sample, with only mostly reduced significance. In particular the shifts to the size distributions of even the small planet sample—which is most affected by the adopted S/N—are qualitatively preserved.

Despite this robustness, η_disc is crucial for the statistical analysis of small planets. The final rows of Table 1 shows that setting η_disc = 1 (but still including the transit probability) gives very discrepant results. These discrepancies are expected due to the low discovery efficiency below 2 R_⊕, as shown in Figure 2.