THE RADIUS DISTRIBUTION OF PLANETS AROUND COOL STARS

In a perfectly idealized survey that is both 100% reliable and 100% complete, the occurrence rate of planets (in a survey, or in a specified bin) is simply

$$\begin{equation} {\rm NPPS} = \frac{N_p}{N_\star }, \end{equation} \tag{ 2 }$$

where N_p is the number of detected planets and N_⋆ is the number of stars surveyed. In practice, however, this must be corrected for both incompleteness and unreliability as follows:

$$\begin{equation} {\rm NPPS} = \frac{1}{N_{\star }} \sum _{i=1}^{N_p} w_i. \end{equation} \tag{ 3 }$$

Here the sum is over all detections and w_i is a weighting factor applied individually to account for the various necessary corrections. Generally, these weights can be thought of as

$$\begin{equation} w_i = \frac{(1-{\rm FPP}_i)}{\eta _i}, \end{equation} \tag{ 4 }$$

where FPP_i is the probability that signal i is a false positive and η_i is an individualized efficiency factor for the detection of planet i. In this work, we calculate the FPP_i according to the procedure in Morton (2012); see Section 3.4 for details.

We thus focus discussion here on the detection efficiency η_i, which is defined by the following thought experiment: If a very large number of planets identical to planet i were distributed randomly around all the stars in the survey, only a fraction η_i could have been detected. Conceptually, this can be factored (following Youdin 2011):

$$\begin{equation} \eta _i = \eta _{{\rm tr},i} \cdot \eta _{{\rm disk},i}, \end{equation} \tag{ 5 }$$

where η_tr is the geometric transit probability and η_disk is the "discovery efficiency": the fraction of planets in this thought experiment with transiting orbital geometries that could have been detected by the survey. Previous Kepler occurrence rate calculations (Howard et al. 2012; Swift et al. 2013; Dressing & Charbonneau 2013), have defined the discovery efficiency as

$$\begin{equation} \eta _{{\rm disk},i} = \frac{N_{\star,i}}{N_{\star }}, \end{equation} \tag{ 6 }$$

where N_{⋆, i} is the number of target stars around which planet i could have been detected. This number is determined by counting the stars around which a transit of planet i (at period P_i) would cause a photometric signal with signal-to-noise ratio (S/N) above some detection threshold (10 for Howard et al. 2012; 7.1 for Dressing & Charbonneau 2013 and Swift et al. 2013).

Though not spelled out explicitly, this way of calculating η_{disk, i} essentially boils down to two general ingredients: a hypothetical "ensemble-of-alternative-scenarios" S/N probability distribution ϕ_{S/N, i} for each planet in the survey, and the discovery efficiency as a function of S/N η_S/N for the transit detection pipeline. Given these two ingredients, Howard et al. (2012), Dressing & Charbonneau (2013), and Swift et al. (2013) all implicitly calculate the following:

$$\begin{equation} \eta _{{\rm disk},i} = \int _0^\infty \eta _{\rm S/N}(s) \cdot \phi _{{\rm S/N},i}(s) ds, \end{equation} \tag{ 7 }$$

where ϕ_{S/N, i} is a normalized probability distribution of S/N (s) for that planet determined by varying the stellar radius and noise properties according to each star in the survey, and η_S/N is assumed to be a step function at some detection threshold.

It is important to note, however, that a step function does not accurately characterize the true detection efficiency of the Kepler pipeline as a function of S/N. Fressin et al. (2013) found that the observed S/N distribution of the Batalha et al. (2013) catalog implied that the true behavior was more like an "S/N ramp," where the discovery efficiency was zero at S/N = 6 and unity at S/N = 16. Internal tests of the Kepler pipeline indicate that the true shape of the η_S/N function is similar to this (P. Tennenbaum 2013, private communication, based on poster at the 2nd Kepler Science Conference). An accurate estimate of this function is crucial to obtaining correct results in an occurrence rate calculation; we present the functional form we use in Section 3.

We also note that the traditional method of computing ϕ_{S/N, i}—simulating a planet of radius r_i transiting all of the survey stars at fixed period P_i—is designed to estimate the joint period–radius distribution of planets. This joint distribution is then summed up over period to find the period—marginalized radius distribution—this is how planet occurrence as a function of radius is determined in Howard et al. (2012), Dressing & Charbonneau (2013), E. Petigura & G. Marcy (2013, in preparation), and Petigura et al. (2013). However, as Figure 1 demonstrates, the smallest planets in the Cool KOI survey are likely only complete out to periods of a few tens of days, whereas many larger planets are detected on larger orbits. Thus, determining the period–marginalized radius function of planets using this traditional approach requires either restricting the analysis to radii above which the survey is substantially complete out to P = P_max, or only using a small fraction of the survey detections (out to periods of only a few tens of days) in order to probe below ∼1 R_⊕.

Zoom In Zoom Out Reset image size

We introduce here an alternative method of computing ϕ_{S/N, i} that allows for better reconstruction of the planet radius function at small radii without having to restrict analysis to a fraction of the available data. Instead of simulating the ensemble of alternative transit scenarios for planet i all at the single fixed period P_i, we suggest that these alternative scenarios could be assigned a distribution of periods according to a reasonable estimate of the true planet period distribution. This effectively amounts to a strategy of "pre-marginalization" that requires an assumption of the period distribution of planets but allows for every planet detection to be treated on an equal basis.

Making this adjustment to how η_{disk, i} is determined—by distributing hypothetical planets around all stars at all periods to calculate ϕ_{S/N, i}—also necessitates rethinking how transit probability is corrected for. That is, construction of ϕ_{S/N, i} must acknowledge that transit probability (as well as S/N) depends on planet period and host star radius. In other words, any process of building up a hypothetical S/N distribution from many instances of simulated transits must ensure that each instance is properly weighted by its actual likelihood of detection. The most straightforward way of ensuring this is to simulate an isotropic cloud of planets around each target star and include only the transiting configurations (e.g., impact parameter b < =1) in the S/N distribution. The transit probability factor η_{tr, i} then becomes the fraction of all these simulated planets that transit—this will typically not be the same as the individual geometric transit probability of planet i, since planet i has a single period P_i while the simulated population has a distribution of periods. In other words, calculating the radius function pre-marginalized over period in this way demands that only a single survey- and period-averaged transit probability be used when calculating the completeness correction for each planet—a function of only the assumed period distribution and the distribution of stellar radii of the target sample.

In all the Kepler planet occurrence calculations to date, the shape of the radius function has been explored only very coarsely, by calculating the occurrence rate in several different radius bins and either fitting a power law or qualitatively commenting on the shape. Howard et al. (2012) found a good fit to an R⁻² power law down to 2 R_⊕, and declined to comment for smaller planets. On the other hand, Fressin et al. (2013) and E. Petigura & G. Marcy (2013, in preparation) note that the occurrence rate of planets increases toward smaller radius but then appears to flatten out below about 2.8 R_⊕. Dressing & Charbonneau (2013) claim that the occurrence rate begins to decrease for planets smaller than 1–1.4 R_⊕, and Petigura et al. (2013) find that the planet occurrence rate decreases for planets smaller than 2 R_⊕.

Investigating the shape of the radius distribution in more detail requires a non-parametric approach, and also should avoid binning. Here we introduce the concept of using a weighted kernel density estimator (KDE) in order to accomplish this.

2.2.1. Weighted Kernel Density Estimation

A standard KDE attempts to estimate the true underlying probability distribution of a sample of data points using a function of the following form:

$$\begin{equation} \hat{\phi }(x) = \frac{1}{N} \sum _{i=1}^N k(x-x_i;\sigma _i), \end{equation} \tag{ 8 }$$

where N is the number of data points and k(x) is a zero-mean, normalized kernel function of arbitrary shape (commonly a Gaussian), with some bandwidth σ_i, that most generally can be different for each data point. This creates a smooth distribution out of a discrete data set, with the degree of smoothness controlled by the width parameter. The choice of width has tradeoffs in both directions: if the kernels are too narrow the estimator will be bumpy, but if they are too wide they can wash out real structure in the distribution. Often the width is selected to be the same for all points based on the number of data points, or sometimes a variable-width kernel is used, e.g., the distance to the nth nearest neighbor. The 1/N normalization factor assures that the integral of this density estimator over the whole parameter space is unity.

In order to use the KDE concept to properly reconstruct the radius function of planets detected in a transit survey, each data point has to be weighted appropriately, leading to a weighted KDE, or wKDE:

$$\begin{equation} \hat{\phi }^{P_{\rm max}}_r(r) = \frac{1}{N_\star } \sum _{i=1}^{N_p} w_i \cdot k(r-r_i;\sigma _i), \end{equation} \tag{ 9 }$$

where w_i = 1/η_i are the appropriately calculated individual weight factors that renormalize the kernels to correct for missing planets, as discussed in Section 2.1. The weights ensure that the shape of the radius function responds appropriately to the individual corrections, and the N_⋆ overall normalization ensures that the integral over all radii will return the NPPS, as desired in Equation (1). We note that because this function does not normalize to unity, it is not strictly a density function, but a "rate function," representing number of planets per unit radius rather than probability density.

Wang & Wang (2007) explore in detail several techniques for selecting the optimal kernel bandwidths, presenting in particular two methods of the form

$$\begin{equation} \sigma _{\rm opt} = 0.9 C n^{-1/5}, \end{equation} \tag{ 10 }$$

where C is a constant based on the sample mean, sample variance, or interquartile range of the data, and n is the number of data points. However, we show in Section 3 that for our particular case, using the actual observational uncertainties in planet radius to control the width of the individual smoothing kernels σ_i results in a smooth distribution, so we do not apply this technique. We do note, however, that if the planet radii were known more precisely, then a method such as this would be necessary in order to quantify an optimum smoothing width to avoid high-frequency features in the estimator.

2.2.2. De-biasing the wKDE and Calculating Variance

The wKDE $\hat{\phi }_r$ described in Equation (9) is an estimator of the true underlying radius distribution ϕ_r, and like any estimator, it has both bias and variance associated with it. These quantities must be determined in order to best understand what the data can tell us about the true distribution. Narsky & Porter (2013) explain how to calculate bias and variance for a standard KDE using a resampling technique called the "smoothed bootstrap"; we adopt this method for our purpose and describe it here.

While a traditional bootstrap technique involves resampling the observed data set with replacement to create new data sets, a "smoothed bootstrap" involves generating new data sets according to the estimated distribution. The density estimator is recalculated for each of these simulated data sets using the same procedure that generated the original estimator. The median offset of these resampled estimators relative to the original estimator (which can be directly observed) will then reflect the bias of the original estimator relative to the true underlying distribution, and the original estimator can thus be corrected accordingly. Similarly, the variance of the estimator can be determined by the scatter of these bootstrapped estimators.

Applying this method to a wKDE in the context of a transit survey is less straightforward than simple resampling according to the derived $\hat{\phi }_r$, but we borrow the same principle. Each smoothed bootstrap resampling data set is generated by simulating a whole new survey: drawing a set of planets according to the estimated distribution, assigning them isotropically distributed orbits around host stars drawn randomly from the target sample, calculating the S/N of each of the resulting transits, and then determining which of these planets would be detectable (using the appropriate η_S/N function). Each of these data sets is then converted into a new wKDE, and the bias may thus be assessed.

Figure 2 illustrates this procedure applied in a toy scenario. The true underlying planet radius distribution—a broken power law—is illustrated with the thick solid black line. The initial wKDE estimate of the radius distribution derived from one realization of "detected" planets is shown as a solid red line. Thin black lines illustrate 200 wKDEs resulting from resampled data sets generated according to the originally calculated wKDE distribution (red line). It is clear that for below 1.5 R_⊕ or so, the bootstrap resamplings underestimate the red line; this mimics the way the red line is biased with respect to the true distribution. So to correct for this, the difference between the solid red line and the median of the bootstrapped distributions is added to the initial wKDE estimate at each value of radius to obtain the de-biased wKDE (dashed red line), which matches the true distribution more closely. The error region around this de-biased estimate is then taken to be the same as the spread in the bootstrapped wKDEs about their median.

Zoom In Zoom Out Reset image size

In addition to illustrating how the bias and variance of these wKDEs may be calculated, this example also touches on two other important points. First, this procedure demonstrates that naive summing of weights to calculate NPPS (Equation (2)) is in fact most generally a biased estimator for the true NPPS—a non-intuitive but significant result. Second, the toy-model true radius distribution used here is a power law that continues to rise all the way down to 0.3 R_⊕. However, because of the detection sensitivity of the survey, only very rarely are any planets smaller than 0.5 R_⊕ detected. Thus, there is generically a turnover in the estimated planet radius distribution—the location of which depends on which planets happen to be observed. In fact, the estimated distribution sometimes even appears to turn over around ∼1 R_⊕—despite the fact that the true distribution continues to rise. Recognizing this is crucial to properly interpreting the radius function derived from Kepler data; we return to discussion of this point in Section 5.

In order to estimate the shape of the true period distribution of planets of all sizes, we make the simplifying assumption that the period distribution of planets is independent of their radii (see Section 5.1 for discussion regarding this assumption). We thus construct the distribution of log P from all the planet candidates in the sample, using a wKDE as described in Section 2.2. For the weights we use only the inverse transit probabilities, and enforce that the whole distribution is normalized to unity, creating the probability density function for log P. For the widths, we use σ = 0.15 (in log P) to create a smooth distribution. This is the period distribution function ϕ_P that we use in the following subsection, shown as the gray shaded region in Figure 3.

Zoom In Zoom Out Reset image size

The S/N of a transit signal is usually defined as follows:

$$\begin{equation} {\rm S/N} = \frac{\delta }{\sigma } \sqrt{N_{\rm tr} \cdot N_{\rm pts}}, \end{equation} \tag{ 11 }$$

where δ is the transit depth, σ is the one-point photometric uncertainty, N_tr is the number of transits observed, and N_pts is the number of photometric points per transit.

In order to construct ϕ_{S/N, i}, the distribution of S/Ns for every conceivable alternative scenario in which planet of radius r_i might have existed in this survey (around every other target star at any other potential period, under our pre-marginalization strategy; see Section 2.1), we use a Monte Carlo simulation. For each target star, 50,000 planets of radius r_i are generated with isotropic orbital inclinations and orbital periods according to the period distribution of Section 3.1. We then calculate the transit signal S/N for each planet in this simulation that has a transiting impact parameter, with δ being the depth of a Mandel & Agol (2002) transit model around that star averaged over the transit duration (using limb darkening parameters for each target star from Claret et al. 2012), σ being the published median 3 hr combined differential photometric precision (CDPP) for that star, N_pts being the total transit duration T_dur/(3 hr), and N_tr being the orbital period of the simulated planet divided by the total amount of time that star was observed by Kepler (number of quarters up until Q12, ×90 days, excepting Q1, which is 33 days). We note that while this formulation ignores details of the noise properties on the exact timescale of transit and the exact observing window function, this prescription for calculating S/N is exactly the use case of the CDPP values as recommended by the Kepler team (Christiansen et al. 2012), especially since most of the targets in this survey are faint and thus white-noise dominated on the timescale of transits. 2.6% of the planets simulated in this exercise—50,000 for each target star—have transiting geometries; this is the survey- and period-averaged transit probability discussed in Section 2.1.

Rather than repeat this exact procedure for every planet, it is sufficient to generate the distribution once for a fiducial planet radius r₀ (e.g., 1 R_⊕), and then create the distribution appropriate for any other planet with radius r by multiplying all the S/Ns in the fiducial distribution by (r/r₀)². Figure 4 illustrates these distributions for 1 R_⊕ and 2 R_⊕. Additionally, η_disk as a function of planet radius r may be determined by simply evaluating Equation (7) along a grid of S/N distributions corresponding to an evenly spaced grid in planet radius. Figure 5 illustrates this result: discovery completeness rises from zero just below r = 0.5 R_⊕ to nearly unity at r = 4 R_⊕.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before being able to apply the formalism discussed in Section 2.2 to derive an estimate of the planet radius function, we must first understand the radii of the detected planets. The radius of a transiting planet is extracted from its transit light curve, the depth of which reveals—approximately—the planet-star radius ratio r_p/R_⋆. However, there are subtle degeneracies between the radius ratio and other parameters of the fit (impact parameter in particular), especially when the duration of the transit begins to be comparable to the photometric integration time. The Kepler long-cadence integration time is 29.4 minutes, and since many of the stars in this particular study are relatively small, the durations of the transits are often only two to three times this. Consequently, the shapes of the transits are artificially smoothed, appearing more V-shaped than they would be with shorter-cadence data, and leading to even greater degeneracy between the planet/star radius ratio and the transit impact parameter.

Thus, in order to understand each planet radius, the NExScI archive catalog values (derived from the Kepler pipeline maximum-likelihood fits) are insufficient: a Markov Chain Monte Carlo (MCMC) approach is required. Our procedure for extracting, detrending and fitting the transit signals will be described in detail in a forthcoming publication (J. Swift et al., in preparation), but we briefly summarize the process here. Pre-search Data Conditioning Simple Aperture Photometry light curves are used from the 23rd data release for our analyses. Transit signals are selected from the light curves based on the catalog values for planet period and duration. Data within 1.5 times the duration of the transit center is preserved and the out of transit data is detrended using a linear model. We then use the emcee python module (Foreman-Mackey et al. 2012) to sample the posterior probability distributions for the parameters in our transit model based on Mandel & Agol (2002) using a quadratic stellar limb darkening with coefficients taken from Claret et al. (2012).

To convert the posterior PDFs for r_p/R_* of each transit signal to a posterior on r_p, we do a joint Monte Carlo sampling from the r_p/R_⋆ distribution and a distribution for R_⋆ according to the derived values from either spectroscopic studies (where available) or Dressing & Charbonneau (2013).

A transiting planet candidate may be an astrophysical false positive rather than a bona fide planet, as has been discussed extensively in literature regarding the Kepler survey. In this work, we calculate the false positive probabilities (FPPs) of each of the Cool KOI candidates using the procedure described in Morton (2012), supplemented by the dispositions on the NExScI archive (for dispositions of "FALSE POSITIVE" we assign FPP = 1). Calculating FPP for a transit candidate inevitably requires a prior assumption about planet occurrence; in Morton (2012), this is quantified by the "specific planet occurrence rate" f_{p, i}, which is an assumed occurrence rate of planets between 2/3 and 4/3 of the radius of planet i.

As the goal of the current study is determining exactly the quantity that is used for the planet prior in the FPP calculations (planet occurrence as a function of radius), and since the radius function study itself depends on FPP, there is an opportunity for iterative convergence. In other words, we may start with initial (not necessarily carefully calculated) values for f_{p, i}, determine the radius distribution $\hat{\phi }_r^{150}$ (as described in the following subsection) using the FPP values implied by this assumption, then recalculate f_{p, i} by integrating over the derived radius distribution. If these "post-calculated" f_{p, i} do not match the initial values used, then the FPPs can be recalculated using these new f_{p, i} values, and the process repeated until the output f_{p, i} match the input. The final result of this process is that of the 130 KOIs in this sample, 11 have FPP > 0.9 (and are thus essentially ignored in generating the radius distribution) and 95 have FPP < 0.10. FPPs for individual KOIs are listed in Table 1. We note that these FPP calculations are based on the Kepler pipeline-derived r/R_⋆ values, so for KOIs in this sample whose MCMC-derived radii are more than 30% discrepant from the Kepler pipeline-derived values, we manually set a lower limit to the FPP of 0.1; these KOIs (15) are indicated in Table 1 as well.

With the r_p posterior distribution p_{r, i}(r) for each KOI in hand, we then construct the wKDE to estimate the planet radius function. The discovery efficiency η_{disk, i} for each planet is given by integrating the curve in Figure 5 over p_{r, i}(r), and multiplied with the average transit probability of 0.026 gives the total efficiency factor η_i. We then calculate FPP values using the iterative method described in Section 3.4, which combined with η_i gives w_i, according to Equation (4), and subsequently the radius function $\hat{\phi }_r^{150}$. We choose the smoothing kernel for each planet to be p_{r, i}(r), without any additional smoothing, and de-bias the estimator and calculate its variance according to the procedure described in Section 2.2.2. Figure 6 illustrates the result, with the 1σ uncertainty region as the gray shaded region. Table 1 lists the planet radii with uncertainties, discovery efficiencies, and FPPs used to generate this radius distribution.

Zoom In Zoom Out Reset image size

Figure 1 illustrates very clearly why the detected population of small planets in the Cool KOI sample is indeed very likely incomplete, showing that where the detected period distribution of the smallest of the Cool KOIs drops off is right around the periods where the known short-period small KOIs would have become undetectable. This is the motivation behind the period redistribution procedure we use to calculate η_{disk, i} in Section 3—correcting for the undetectable longer-period small planets. However, the nature of this correction as applied in this work—using the implied all-planet period distribution for each planet—merits some discussion.

There are certainly both physical reasons and observational suggestions to believe that the planet period distribution is not completely independent of radius. In particular, Howard et al. (2012) finds (shown in their Figure 6) that the fraction of short-period planets that are large (4–8 R_⊕) is smaller than the fraction of longer-period planets that are large; in other words, the period distribution of larger planets decreases (heading toward shorter periods) sooner than does the distribution of smaller planets (2–4 R_⊕). Dong & Zhu (2012) present a similar finding. While there is not yet compelling evidence that this same effect has been detected for planets smaller than 2 R_⊕, simple physical considerations such as increasing stellar insolation (L. Weiss & G. Marcy 2013, in preparation) might reasonably contribute to a dearth of larger planets on short-period orbits. However, as there is no corresponding clear physical explanation for the absence of smaller planets in longer orbits, it is reasonable to assume that they do in fact exist, and that their period distribution might resemble the period distribution of the larger planets that are detected in such orbits.

In order to explore this in detail, we repeat the analysis assuming a modified period distribution, shown in the inset of Figure 7. Rather than using a wKDE built from the observed periods of the detected planets, we use a log-flat distribution for periods greater than 10 days, with an exponential cutoff short of 10 days. Qualitatively, this seems to approximately match the observed distribution, except that instead of tailing off toward longer periods it remains constant. The radius function that results from assuming this period distribution is plotted as the dashed line in the main panel of Figure 7, and differs only very slightly from the original analysis, within the noise inherent in the de-biasing procedure. The reason for this negligible difference is that the bulk of the correction for small planets happens when it is assumed that the planets observed with ≲1 day periods can also exist at periods of ∼10s of days—in other words, a 0.7 R_⊕ planet is just as undetectable at a period of 40 days as at 100 days. We thus conclude that our results are not very sensitive to the details of the longer-period–radius distribution, but rather depend only on the assumption that small planets do indeed exist at periods beyond which they are detectable in a manner roughly similar to that of larger planets.

Zoom In Zoom Out Reset image size

The method we have used to calculate the period–marginalized radius function contains several significant differences from previous planet occurrence rate studies. In particular, we use a survey- and period-averaged approach to calculating both the discovery efficiency of a planet and the transit probability correction factor, as opposed to the more traditionally used method of correcting each planet individually—that is, simulating the planet around other stars at only its observed period, and using each planet's individual transit probability. We summarize this difference as calculating the radius function while "pre-marginalizing" over period.

Figure 8 illustrates the difference in discovery efficiency calculated using the method used in this work (η_{disk, i} as described in Section 3) and that which would be calculated by only simulating planets around other stars at the single detected period: $\eta _{{\rm disk},i}^{\rm simple}$. In this figure, the size of the points is proportional to planet radius, with the annuli representing the ±1σ uncertainties from the radius posteriors. Predictably, planets discovered at short periods—predominantly small—get smaller η_{disk, i} (larger correction factor) than the simple "post-marginalized" calculation; this is because the majority of the period-redistributed simulated planets will be at longer periods than the original, and thus have lower S/N. Planets discovered at longer periods tend to have the opposite effect.

Zoom In Zoom Out Reset image size

We also note that to get η_{disk, i} in our calculations, we integrate the posterior distribution of planet radius over the efficiency curve in Figure 5; for $\eta _{{\rm disk},i}^{\rm simple}$, we simulate the alternative scenarios as having the planet radius fixed to be the median of its posterior—this explains the scatter of this difference around a monotonic relationship to period. Additionally, the efficiency curve we use is not just a simple S/N threshold, as has been used in some previous studies.

In addition, in order to correct for the planets missed by Kepler due to random orbital orientations, we do not correct each planet individually for its observed transit probability; rather, we use the same single survey- and period-averaged transit probability for every planet. The justification for this is that only way to properly determine η_i is through the full Monte Carlo procedure discussed in Section 3, where planets are simulated isotropically, the S/N distribution of the ones with transiting orientations are used to calculate η_{disk, i}, and the fraction of the whole simulation that transits is used as η_{tr, i}. Alternatively, one could imagine ignoring the factorization entirely and just assigning S/N = 0 (with detectability = 0) to the non-transiting planets in order to directly calculated η_i, and the result would be the same (though less intuitive).

Figure 9 illustrates how the results of this work would differ if we used the same data, but did our calculations according to the more widely used methods of completeness correction. The dotted line illustrates a calculation in the style of Dressing & Charbonneau (2013) or Howard et al. (2012). Each planet is corrected individually for its own transit probability rather than using a global average, and the detection efficiency curve is taken to be a step function at S/N = 7.1. We see that the occurrence rate of small planets is significantly underestimated. The dashed line is an improved version of this calculation, using the detection efficiency as a function of S/N from Figure 5, but still correcting for transit geometry using only individual transit probabilities and without period redistribution in the η_disk calculation (this is analogous to the method used by Petigura et al. 2013). We do not de-bias either the dotted or dashed distributions. The overall normalization is lower than our calculations (solid line) and the qualitative shape of the radius distribution has shifted, showing a peak around 1.2 R_⊕ and a clear decline below, rather than the continued rise to below 1 R_⊕ that we find.

Zoom In Zoom Out Reset image size

This qualitative reason for this discrepancy is quite simple: the dashed and dotted lines in Figure 9 are plotted down to below the planet radius where the survey is able to reliably detect planets on 150 day orbits (in this case, about 1 R_⊕). Whereas the method used in this paper "pre-marginalizes" over period, enabling extrapolation of the smaller-radius population to longer periods, the "simple" individual method does not do this. The simple method implies "post-marginalization," which can only be valid if the survey is actually sensitive enough to planets in the longest-period and smallest-radius corner of the period–radius space under consideration. So the simple method is not wrong if properly applied, but it cannot be used to estimate the radius function down to as low a radius as the method advocated in this paper. Additionally, it is easy to see from Figure 9 how the shape of the simple-method estimator is much more sensitive to random fluctuations in the data (due to the large correction factors for long-period planets) than the pre-marginalized estimator; this can also cause spurious peaks in the distribution that are strongly dependent on the particular realization of the data that we happen to observe.

While we have demonstrated and discussed differences between the planet radius function derived using two different techniques of calculating individual weights, the question still remains whether or how well either can actually accurately reconstruct the true radius distribution.

To explore this, we run a Monte Carlo experiment where we first generate 100 different transit survey data sets with planet frequency and radii drawn from a known distribution—a broken power law normalized to three planets/star on the interval [0.3,4.0] R_⊕, proportional to $r_p^{0.5}$ for r_p < 1.5 R_⊕ and to $r_p^{-3}$ for r_p > 1.5 R_⊕ (the same distribution used in the de-biasing demonstration in Figure 2). The periods of these simulated planets are assigned according to the log-flat/exponential distribution illustrated in Figure 7. For each of these mock data sets, we derive and de-bias the wKDE estimator $\hat{\phi }_r^{150}$ two different ways: first calculating weights w_i according to the procedure described in Section 3, and second calculating weights according to the "simple" prescription, without period redistribution and using individual transit probabilities.

Figure 10 illustrates the results of this experiment. Each thin red line is a de-biased wKDE derived using the methods presented in this paper to calculate detection efficiencies; each thin blue line is a de-biased wKDE using the individual period/individual transit probability (simple) method. The thicker lines and colored bands indicate the respective median and ±1σ ranges from this ensemble. The white dashed line is the true underlying distribution from which the planet radii are drawn.

Zoom In Zoom Out Reset image size

There are several important points this experiment demonstrates. First, both methods correctly recover the underlying distribution nearly perfectly for radii larger than 1.5 R_⊕, validating the accuracy of the technique. Below 1.5 R_⊕, the two methods begin to diverge, with ours sticking closely to the true distribution until r_p ≲ 0.5 R_⊕, but the simple method beginning to significantly underestimate planet occurrence by r_p ≲ 1 R_⊕. Indeed, the typical result of the simple method is an estimate of a turnover in the radius function at ∼1 R_⊕—a strong qualitative discrepancy from the true distribution.

We also note that at the very low end of the distribution (r_p ≲ 0.5 R_⊕), neither method correctly reconstructs the continued rise, though the period-redistribution method does still keep the true distribution within the 1σ percentile range while the simple method does not. This may be understood by realizing that if a very small planet does happen to be detected (which will only occur in rare cases since the detection efficiency is only ∼5% at r_p = 0.5 R_⊕), it will necessarily be at a very short orbital period, so the weight factor it receives without period correction will be strongly underestimated. This confirms our qualitative understanding of the difference between the two methods—the post-marginalized simple strategy can only be valid down to a radius where the survey is decently complete at the maximum allowed period (this is apparently somewhere between 1 and 1.5 R_⊕), whereas the pre-marginalized method enables confident extrapolation to significantly below this point.

Perhaps most dramatically, Figure 10 illustrates that even if the true distribution continues to rise down to arbitrarily small radii, the estimator for its shape will typically turn over significantly above the actual detection limit of the survey, and begin to flatten at even larger radius. In other words, our derived radius distribution from actual Kepler data that we show in Figure 6 could very plausibly reflect a true underlying distribution that keeps rising continuously down to below 0.5 R_⊕—despite the fact that a log-binned histogram of the estimated distribution looks like it turns over below 1.5 R_⊕.

Finally, we note that the analysis of Petigura et al. (2013) has found, in apparent contradiction with our results, that planets in the 2.0–2.8 R_⊕ radius bin are more common than planets in the 1.4–2.0 R_⊕ bin, which are in turn more common than in the 1.0–1.4 R_⊕ bin. There are several possible explanations for this discrepancy. First, the bins in question are logarithmically spaced, which immediately exaggerates the presence of any turnover, though this alone is not sufficient to explain the difference in results. Second, while the bins used by Petigura et al. (2013) are chosen to all have significant numbers of detected planets, the completeness does vary from ∼70% to ∼10% within the smallest, longest-period bin used for the radius function reconstruction, potentially leaving room for some of the effects that distinguish the blue curve from the red in Figure 10. Another contributing explanation may also be that the only planets counted in Petigura et al. (2013) were the most-detectable planets in each system, which will lead to preferential underestimation of occurrence in the smallest radius bins. Or most simply, this difference may just be a function of the stellar target sample considered—Petigura et al. (2013) used only solar-like (G/K) stars in their target sample, whereas we use only stars with T_eff < 4000 K. If this is the reason, further study of Kepler results should distinguish a radius function that changes shape with increasing host-star temperature.