Searching the Entirety of Kepler Data. II. Occurrence Rate Estimates for FGK Stars

We describe our input stellar and planet catalogs in Section 2. A full description of the process to create our planet catalog is the content of Kunimoto et al. (2020, hereafter Paper I). In Section 3, we describe our determination of both search and vetting completeness using injection/recovery tests. In Section 4, we give an overview of the ABC methodology, and we discuss our application of ABC to exoplanet occurrence rates in Section 5. We make our code available for public use on Github² under the BSD 3-Clause License (Kunimoto & Matthews 2020).

Our overall results are presented in Section 6, in which we discuss the dependence of exoplanet occurrence rates on stellar effective temperature (Section 6.2), planet radius (Section 6.3), and orbital period (Section 6.4). We also describe our incorporation of a simple catalog reliability model to assess the impact of a nonzero FP rate and better constrain our estimates (Section 6.5). In Section 7, we present both baseline and reliability-incorporated results over the potentially habitable, rocky exoplanet parameter space. Finally, in Section 8, we review the limitations of our methodology and give a final recommended ${\eta }_{\oplus }$ estimate.

We started with the 197,096 Kepler stars in the Q1–Q17 DR25 stellar catalog (Mathur et al. 2017), and calculated limb-darkening coefficients using T_eff, log g, and [Fe/H] from Claret & Bloeman (2011). With the arrival of stellar parallaxes in Gaia Data Release 2 (DR2), Berger et al. (2018) produced improved radii for 177,911 targets, yielding an average radius precision of less than 10% for most Kepler stars. Given that a fully updated set of stellar properties has not yet been released, we used these radii in tandem with the Mathur et al. (2017) catalog and only kept stars present in both catalogs.

We removed stars flagged in Berger et al. (2018) as likely binary stars (BIN flag = 1 or 3; 174,769 stars remained). We did not remove those flagged as binaries due to companions revealed with high-resolution imaging (BIN flag = 2), as these observations were only available for a subset of stars. Because the focus of our study is FGK dwarfs, we also removed stars with Evol flag >0 (116,637 stars remained), which indicate that they are unlikely to be on the main sequence.

To ensure each star's light curve had enough data to allow for the discovery of long-orbit planets, we required that the time length of the data (T_obs) is at least 2 yr, and the duty cycle (f_duty) was at least 0.6; in other words, at least 60% of the observations must be filled. After these cuts, 100,823 stars remained.

Lastly, we retained only FGK stars by using suggested T_eff limits from Pecaut & Mamajek (2013). This left 40,010 F- (6000 ≤ T_eff < 7300 K), 39,173 G- (5300 ≤ T_eff < 6000 K), and 17,097 K-type (3900 ≤ T_eff < 5300 K) stars, for a total of 96,280 stars in our sample.

Some stars in this sample may have been chosen as targets for reasons other than the Kepler exoplanet search program, such as for asteroseismology. These stars would be expected to exhibit different noise and variability properties than typical main-sequence stars, which could introduce a systematic bias in our results relative to studies that focus on only exoplanet search targets. We checked the investigation ID of each star in our sample using the Kepler Data Search & Retrieval form on MAST,³ and found that 290 did not have an "EX*" ID. In other words, 99.7% of our 96,280 FGK stars were selected for the exoplanet search program, and we do not expect a significant bias to be present.

Our full search and vetting pipeline is described in Paper I. In short, we obtained Q1–Q17 DR25 long-cadence PDC light curves from the Mikulski Archive for Space Telescopes (MAST).⁴ We detrended each light curve using the detrend5 routine from the Kepler Transit Model Codebase (Rowe 2016), where each observation was corrected by fitting a cubic polynomial to a segment (typically two days wide) centered on the time of measurement. We then 5σ-clipped the data, removing outliers only in the positive flux direction so as to leave deep transits untouched and removed data near data gaps. Then, we used a Kovács et al. (2002) box least-squares (BLS) algorithm to search for potential transits. After identifying an event in the light curve, we calculated its S/N by dividing the mean transit depth by the standard error of the mean, giving

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{\sqrt{N}}{\sigma }{T}_{\mathrm{dep}},\end{eqnarray} \tag{ 1 }$$

where T_dep is the mean transit depth, σ is the standard deviation of the observations, and N is the number of in-transit data points. This definition is comparable to the "effective" S/N described in Kovács et al. (2002). Following Rowe et al. (2014), we estimated σ using the standard deviation of all out-of-transit observations—defined as data outside of two transit durations of the center of the detected signal—and used the median absolute deviation (MAD) with σ = 1.48MAD (Hoaglin et al. 1983) to be more robust to outliers.

To define transit candidates (TCs), we followed the suggestion of Kovács et al. (2002) that the threshold for a significant detection with the BLS algorithm is S/N = 6. We also required at least three transits and the passing of an initial vetting stage to reject false alarms caused by instrumental and astrophysical systematics.

Each TC was passed through a vetting pipeline, involving both machine and manual triage. Automated candidacy tests were used to flag both noise false alarms and astrophysical FPs, while visual inspection was used as a "reality check" to confirm each surviving TC as a planet candidate (PC).

Around the 96,280 FGK stars considered, we identified 2623 PCs that matched with already known PCs as listed on the NASA Exoplanet Archive,⁵ defined as Kepler Objects of Interest (KOIs) with either a CONFIRMED or CANDIDATE disposition. Additionally, we introduce eight previously unknown candidates from our Paper I search, for a total of 2631 planets in our full catalog. By comparison, Kepler's Q1–Q17 DR25 pipeline identified 2829 PC KOIs corresponding to this stellar sample.

2.2.1. Confirmed and Candidate KOIs Missed

2.2.2. New PCs

2.2.3. Planet Properties

2.2.4. Dilution

The planet radius R_p can be determined from the fitted parameter R_p/R_s by multiplying by the known stellar radius. However, there may be one or more nearby stars that contribute light to the Kepler aperture, causing the measured transit depth to be diluted. In these cases, the planet radius can be underestimated. We make the assumption that the planet orbits the brighter primary star, in which case we apply a correction factor to the planet radius in the form of

$$\begin{eqnarray}&&{R}_{p,\mathrm{corr}}={R}_{p}\sqrt{1+{10}^{-0.4{\rm{\Delta }}m}},\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Delta }}m={m}_{\sec }-{m}_{\mathrm{pri}}$ is the Kepler magnitude difference between the primary and secondary star. For more than one companion per star, the previous equation becomes

$$\begin{eqnarray}&&{R}_{p,\mathrm{corr}}={R}_{p}\sqrt{1+\sum _{i=1}^{N}{10}^{-0.4{\rm{\Delta }}{m}_{i}}},\end{eqnarray} \tag{ 3 }$$

where the sum is for N companion stars with magnitude differences Δm_i.

We used the high-resolution imaging results from the Kepler Follow-Up Observation Program (Furlan et al. 2017) to correct the radii of planets around stars with a potential companion within 4'', the size of a Kepler pixel. Furlan et al. (2017) compiled observations for a total of 3557 KOIs, including those observed in the first three Robo-AO surveys (Law et al. 2013; Baranec et al. 2016; Ziegler et al. 2017), and provided a weighted average of correction factors across a variety of bands for 1891 KOIs with companions.

Ziegler et al. (2018) presented a fourth Robo-AO survey for 532 KOIs published after Furlan et al. (2017). Their results were provided as Δm in the LP600 band, which we approximate to be equal to the Kepler band for use in Equations (2) and (3). We also used our own adaptive optics imaging follow-up described in Paper I for three of our new FGK PCs (KIC-6126245 b, 6782399 b, and 7747788 b), none of which had a nearby stellar companion.

In total, 2578 of the 2631 PCs (98.0%) in our FGK sample had high-resolution imaging observations, and we applied correction factors to 679.

2.2.5. Final Planet Catalog

Our focus for this paper is on the occurrence rates of planets in a period–radius grid spanning orbital periods 0.78125 < P < 400 days and radii 0.5 < R_p < 16.0 R_⊕. Lower and upper limits on these properties were chosen so as to split the grid into logarithmically spaced bins comparable to bins used in previous grid-based works (e.g., Howard et al. 2012; Petigura et al. 2013; Mulders et al. 2015a). After applying the radius correction factors and including only candidates that fit these criteria, our final planet catalog involved 557 candidates around F-type stars, 1276 around G-type stars, and 700 around K-type stars, for a total of 2533 PCs. Table 3 summarizes the sizes of each star and planet sample, while Figure 1 shows the distribution of planets based on orbital period and radius.

Zoom In Zoom Out Reset image size

Table 3.  Number of Stars and Planets by Stellar Type

Type	${T}_{\min }$ (K)	${T}_{\max }$ (K)	${N}_{\mathrm{stars}}$	${N}_{\mathrm{planets}}$	${N}_{\mathrm{DR}25}$
FGK	3900	7300	96,280	2533	2700
F	6000	7300	40,010	557	639
G	5300	6000	39,173	1276	1338
K	3900	5300	17,097	700	723

Note. ${N}_{\mathrm{DR}25}$ gives the number of planets found in DR25 around the same sample of stars for comparison.

Download table as:  ASCIITypeset image

Bayesian inference is an increasingly popular approach of statistical inference on unknown parameters. In this framework, Bayes' theorem is used to estimate the posterior probability distribution $P({\boldsymbol{\theta }}| D)$ of a model with parameters ${\boldsymbol{\theta }}$ given the data D,

$$\begin{eqnarray}&&P(\theta | D)=\displaystyle \frac{P(D| \theta )P(\theta )}{P(D)},\end{eqnarray} \tag{ 11 }$$

where $P(D| {\boldsymbol{\theta }})$ is the likelihood function, indicating the compatibility of the data given the model; $P({\boldsymbol{\theta }})$ is the prior probability, representing initial beliefs toward the model; and P(D) is a normalization constant. The best-fit model parameters can be estimated from $P({\boldsymbol{\theta }}| D)$ such as by finding the posterior mode (most probable values of ${\boldsymbol{\theta }}$) or posterior median (50th percentile), with credible intervals representing our uncertainty about the model parameters.

For simple models, the likelihood function can typically be derived analytically. However, for more complex models, the likelihood may be unknown or too computationally expensive to evaluate. It is in these cases that the "likelihood-free" method of ABC steps in as an effective and rigorous way of performing an approximate Bayesian analysis.

ABC circumvents the need for a likelihood function by using our prior information along with an ability to simulate, or "forward model," the observed data under investigation. By simulating a large number of data sets and quantifying the "distance" between each data set and the observed data set, the distribution of model parameters that provides the best matches can be determined. This distribution serves as an approximation to the posterior probability distribution.

The specific form of ABC used here is the population MC (ABC-PMC) algorithm proposed by Beaumont et al. (2009), wherein multiple generations of simulated data are created and an adaptive importance sampling scheme is used to evolve the ABC posterior. We use the ABC-PMC algorithm implemented in cosmoabc, a Python ABC Sampler (Ishida et al. 2015), which is summarized here.

To initialize the ABC-PMC algorithm, we draw a set of M values from the prior distribution, called "particles," $\{{{\boldsymbol{\theta }}}^{i}\}$ with i ∈ [1, M]. M is chosen to be much larger than N, the number of samples needed to characterize the prior. For each particle, we generate a simulated data set ${D}_{S}^{i}$ and use a distance function ρ to calculate the distance between the simulated and real data set, ${\rho }^{i}=\rho (D,{D}_{S}^{i})$. From the whole set of M particles, we keep only the N particles with the smallest ρⁱ. These constitute the zeroth "generation" (S_t=0), and the 75% quantile of all ρ ∈ S_t=0 gives the distance threshold for the next iteration (_t=1). Each particle is assigned an equal weight, ${W}_{t=0}^{j}=1/N$, for $j\in [1,N]$.

An importance sampling technique is used to produce subsequent generations (t > 0). We draw a trial particle $({{\boldsymbol{\theta }}}_{\mathrm{try}})$ from the previous generation S_t − 1 with weights W_t − 1, and use it to simulate a catalog and find its associated distance, ρ_try. We store ${{\boldsymbol{\theta }}}_{\mathrm{try}}$ to the current generation S_t if ρ_try ≤ _t. This process is repeated until S_t is filled with N accepted particles. We then calculate the weights of each particle as

$$\begin{eqnarray}&&{W}_{t}^{j}=\displaystyle \frac{P({{\boldsymbol{\theta }}}_{t}^{j})}{{\sum }_{i=1}^{N}{W}_{t-1}^{i}N({{\boldsymbol{\theta }}}_{t}^{j};{{\boldsymbol{\theta }}}_{t-1}^{i},{C}_{t-1})},\end{eqnarray} \tag{ 12 }$$

where $P({{\boldsymbol{\theta }}}_{t}^{j})$ is the prior probability distribution calculated at ${{\boldsymbol{\theta }}}_{t}^{j}$, and $N({{\boldsymbol{\theta }}}_{t}^{j};{{\boldsymbol{\theta }}}_{t-1}^{i},{C}_{t-1})$ represents a Gaussian probability density function (PDF) centered at ${{\boldsymbol{\theta }}}_{t-1}^{i}$ with covariance matrix built from ${S}_{t-1}$ and calculated at ${{\boldsymbol{\theta }}}^{j}$.

Following the determination of the new weights, the algorithm repeatedly produces new generations until subsequent iterations no longer significantly change the ABC posterior. In cosmoabc, this convergence occurs when the number of draws necessary to construct a generation is much larger than N.

We assign independent uniform priors for each occurrence rate over $[0,{f}_{\max ,p,r})$. The upper limit for each bin is

$$\begin{eqnarray}&&{f}_{\max ,p,r}=C\times {\mathrm{log}}_{2}\left(\displaystyle \frac{{P}_{\max ,p}}{{P}_{\min ,p}}\right)\times {\mathrm{log}}_{2}\left(\displaystyle \frac{{R}_{p,\max ,r}}{{R}_{p,\min ,r}}\right),\end{eqnarray} \tag{ 13 }$$

with C = 2, small enough that proposals with more than three planets per factor of 2 in period are rare (Hsu et al. 2019). This is consistent with expectations based on long-term orbital stability.

It is within our exoplanet population simulator that many of the complexities that make a likelihood function infeasible to compute are able to be incorporated into the determination of occurrence rates.

One such complexity is the existence of selection effects. Youdin (2011) outlined three main selection effects to be accounted for as part of robust exoplanet population analysis. These are quantified as detection efficiencies, η, which give the ratio of detections to actual planets: (i) η_tr, the transit probability that the planet crosses our line of sight to the star; (ii) η_rec, the efficiency at which the detection pipeline recovers the planet; and (iii) ${\eta }_{\mathrm{fp}}=1/(1-{r}_{\mathrm{fp}})$, where r_fp is the rate of FP events that are detected as planets. The net detection efficiency of a given planet is found by multiplying all of the above efficiencies together. For our baseline results, we assumed the FP rate is low enough that it can be ignored for simplicity (${\eta }_{\mathrm{fp}}=1$). However, we discuss potential implications of this assumption in Section 6.5.

Importantly, these selection effects change on a per-star basis. For instance, our completeness model depends on both the physical properties of a star and the characteristics of its associated Kepler light curve. Our forward model allows us to take these into account and find a specific completeness for a planet around a specific star, with little sacrifice of computational efficiency. By comparison, studies that have used likelihood functions in occurrence rate statistics such as Burke et al. (2015) and Zink et al. (2019) have had to utilize star-averaged detection efficiencies that depend only on P and R_p, as incorporating information about individual stars would be too computationally expensive. In these cases, two planets with the same period and radius but host stars with vastly different properties would still be assigned the same completeness.

Furthermore, given that we are focused on specific period and radius bins, occurrence rates may be sensitive to the accuracy of a planet's membership in its correct bin. While orbital period is typically known to an accuracy of minutes or better, uncertainties in planet radius are significantly larger. First, measurement errors caused by fitting a transit model to a noisy light curve can cause the fitted ratio R_p/R_s to differ from its true value. Second, and more significantly, uncertainty in the star's radius used to derive R_p from R_p/R_s directly leads to uncertainty in the planet's radius, even if R_p/R_s is known exactly. As a result, a planet's "observed" radius bin may differ from its true radius bin, particularly if it is near the boundary between two bins. As in Hsu et al. (2019), our forward model is able to take into account these measurement uncertainties by simulating both true and observed stellar and planetary radii, while most other studies assume that a planet's properties are known exactly.

5.2.1. Step 1: Generate Planets

5.2.2. Step 2: Calculate Selection Effects

Transit Probability: Many planets will be undetected simply because they do not cross our line of sight to the star. We use the planet's semimajor axis $a={({{GM}}_{s}{P}^{2}/4{\pi }^{2})}^{1/3}$ and impact parameter i drawn previously to determine the planet's impact parameter,

$$\begin{eqnarray}&&b=\displaystyle \frac{a\cos i}{{R}_{s}},\end{eqnarray} \tag{ 14 }$$

requiring that b ≤ 1. In other words, the planet transits if the center of the planet passes inside the disk of the star. As in Hsu et al. (2019), we ignore the small number of transiting planets with b > 1, as large impact parameters are often associated with grazing eclipsing binaries and these planets are likely to be flagged as FPs. Thus, we set

$$\begin{eqnarray}{\eta }_{\mathrm{tr}}=\left\{\begin{array}{ll}1 & b\leqslant 1\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

Recovery Efficiency: We estimate the recoverability of each planet by taking into account pipeline search completeness, vetting completeness, and the probability that at least three transits occur in the Kepler window.

For search and vetting completeness, we use the combined search and vetting model as outlined in Section 3. We follow the same process outlined in Section 3 to estimate each simulated planet's transit S/N and N_tr to determine the corresponding P_det.

For the window probability, we use the binomial probability function described in Burke et al. (2015),

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{win},\geqslant 3} & = & 1-{\left(1-{f}_{\mathrm{duty}}\right)}^{M}-{{Mf}}_{\mathrm{duty}}{\left(1-{f}_{\mathrm{duty}}\right)}^{M-1}\\ & & -\displaystyle \frac{M(M-1)}{2}{f}_{\mathrm{duty}}^{2}{\left(1-{f}_{\mathrm{duty}}\right)}^{M-2},\end{array}\end{eqnarray} \tag{ 16 }$$

where M = T_obs/P. Thus, we find the total recovery efficiency for a given planet as

$$\begin{eqnarray}&&{\eta }_{\mathrm{rec}}={P}_{\det }{P}_{\mathrm{win},\geqslant 3}.\end{eqnarray} \tag{ 17 }$$

5.2.3. Step 3: Simulate Detected Exoplanet Population

We determine if a planet is detected by drawing from a Bernoulli distribution with probability

$$\begin{eqnarray}&&{\eta }_{\mathrm{tot}}={\eta }_{\mathrm{tr}}{\eta }_{\mathrm{rec}}.\end{eqnarray} \tag{ 18 }$$

At this point, we remove all planets flagged as undetected from the simulation and focus the remainder of our analysis on the recovered population.

5.2.4. Step 4: Incorporate Planet Radius Uncertainty

5.2.5. Step 5: Compare to Observed Population

We generate summary statistics for each bin in both observed and simulated catalogs,

$$\begin{eqnarray}&&{s}_{k}=\displaystyle \frac{{N}_{k}}{{N}_{s}},\end{eqnarray} \tag{ 19 }$$

where N_k is the number of planets in the kth bin. We use the fraction of planets per star rather than the absolute number of planets so as to allow for differences in the choice of N_s between catalogs. For instance, we could choose to run a quick inference by comparing our search results from all 96,280 FGK stars with a catalog that simulates planets around only 10,000 FGK stars.

It is at this point that we apply our distance function to quantify the distance between summary statistics and thus assess the agreement between the simulated and observed planet catalogs.

When modeling only a single period–radius bin at a time, such as in Hsu et al. (2018), the summary statistic for each catalog is scalar. The choice of distance function may be simply

$$\begin{eqnarray}&&\rho ({s}_{\mathrm{obs}},{s}_{\mathrm{sim}})={\left({s}_{\mathrm{obs}}-{s}_{\mathrm{sim}}\right)}^{2},\end{eqnarray} \tag{ 20 }$$

where $s=N/{N}_{s}$ is calculated for only the single bin of interest, and obs and sim refer to the observed and simulated catalogs respectively. However, when fitting multiple bins simultaneously, we use the distance suggestion of Hsu et al. (2019),

$$\begin{eqnarray}&&\rho ({s}_{\mathrm{obs},k},{s}_{\mathrm{sim},k})=\sum _{k}\displaystyle \frac{| {s}_{\mathrm{obs},k}-{s}_{\mathrm{sim},k}| }{\sqrt{{s}_{\mathrm{obs},k}+{s}_{\mathrm{sim},k}}},\end{eqnarray} \tag{ 21 }$$

inspired by the Canberra distance (Lance & Williams 1967). Hsu et al. (2019) found that this distance allowed ABC to converge more rapidly than other tested functions. This function also weights the absolute value of the differences in s_k by the square root of the sum, resulting in a similar fractional error in occurrence rates for all bins rather than a similar absolute error.

With our ABC framework set, we justified the number of bins to fit at once and verified that the algorithm was able to recover occurrence rates accurately with the appropriate choices.

Had we not incorporated planet radius uncertainty into the forward model, the placement of each simulated planet into a specific period–radius bin would be without ambiguity. The occurrence rates of each bin would not affect those of others, and thus fitting only one bin at a time would be an obvious choice due to computational efficiency.

Because our simulator takes into account measurement error, it may place a planet into a radius bin different from its true bin. If two neighboring bins have different occurrence rates, the number of planets exchanged across the radius boundary may be asymmetric. Furthermore, the edge bins being fit will display noticeable bias. As the simulator does not simulate planets with true radii above the upper limit of the top bin and below the lower limit of the bottom bin, the exchange of planets over these radii limits will be strictly one sided, and their occurrence rates will be overestimated.

These considerations necessitate the fitting of multiple bins simultaneously. However, fitting more parameters comes at the cost of the performance of the ABC-PMC algorithm, as it becomes less likely that the proposed values for all parameters will result in good agreement between the observed and simulated catalogs. Both the width of the ABC posterior and the computational time required for the algorithm to achieve convergence will increase significantly.

To explore these issues, we used our forward model to simulate 10 "true" planet catalogs in the $6.25\mbox{--}12.5$ day period range using the full 96,280 star sample. We used occurrence rates of $f=\{0.06,0.05,0.04,0.03,0.02,0.01,0.005\}$ for bins with boundaries ${R}_{p}=\{0.35,0.5,0.71,1,1.41,2,$ $2.83,4\}$ R_⊕. We fit a total of one, three, five, and seven bins simultaneously, centered on the $1.25\mbox{--}1.5{R}_{\oplus }$ bin (f = 0.03, or 3%), for each simulated catalog. As inputs to cosmoabc, we set the number of particles for the initial generation at 500 and all subsequent generations at 200. We considered the system converged when at least 2000 draws (10 times the size of each generation) were required to construct the next generation.

Figure 3 shows the final ABC posterior for the 6.25–12.5 day, $1\mbox{--}1.41{R}_{\oplus }$ occurrence rate after each run. As expected, the one-bin fit consistently overestimates the true occurrence rate due to the fact that simulated planets can only leak out of the bin. The average absolute difference between the ABC posterior median and the true 3% occurrence rate was 0.98%. The three-, five-, and seven-bin fits all show significant improvement, with ABC posteriors well clustered around the true occurrence rate. Average absolute differences were 0.24%, 0.15%, and 0.17%, respectively. We also observe the expected widening of the ABC posterior with more bins.

Zoom In Zoom Out Reset image size

When examining the results for the edge bins in each multibin-fit run, we confirmed that they tended to be overestimated compared to the interior bins. This was especially apparent for the bottom-edge bins, likely due to the fact that they were assigned the highest occurrence rates and thus had more outward leakage of planets than inward. These considerations prompted us to exclude the results for the two edge bins when performing multibin fits, and only report the results for the interior bins.

Overall, we agree with the conclusions of Hsu et al. (2019) that five to seven radius bins are the optimal choice, and that one should be careful when considering the results of edge bins.

We compare our occurrence rates with those reported by other studies, choosing F17, Mulders et al. (2015a), Petigura et al. (2013), and Fressin et al. (2013) as targets for comparison on the basis of having the most similar period and radius ranges. Given that occurrence rate studies calculate occurrence rates with different methods, account for completeness in different ways, use catalogs based on different amounts of available Kepler photometry, and more, direct comparison is difficult. Additionally, Petigura et al. (2013) inferred occurrence rates using only the first planet found in each system. Consequently, their occurrence rates are not estimates of the average number of planets per star as in the other works. Recognizing these challenges, F17 compared the ratios of occurrence rates between bins rather than the absolute occurrence rates of individual bins. We adopt the same approach here.

Table 5 shows our FGK results compared to Fressin et al. (2013; FGK), Mulders et al. (2015a; FGKM), and F17 (FGK) while Table 6 shows our G results compared to Petigura et al. (2013; GK) and Mulders et al. (2015a; G). We find that the main discrepancy is that the ratios of occurrence rates for 2–4 R_⊕ planets to 1–2 R_⊕ planets, across all period ranges and for both FGK- and G-type stars, are higher than all previous works. This is most noticeable when comparing to the older studies of Petigura et al. (2013) and Fressin et al. (2013). Petigura et al. (2013) found that 1–2 R_⊕ planets within 50 days are 1.2 times more common than those with radii 2–4 R_⊕. Fressin et al. (2013) found a similar ratio of 1.1. Meanwhile, Mulders et al. (2015a) found that planets with 1–2 R_⊕ are less common, with ratios of 0.9 (FGKM stars) and 0.6 (GK). Our results are even more favored toward the larger radius bin, with ratios of 0.6 (FGK) and 0.4 (G).

Table 5.  Rough Comparisons for FGK-type Stars

Rp (R⊕)	P (days)	This Work (%)	F17 (%)	M15 (%)	F13 (%)
1–2	<50	${16.2}_{-1.1}^{+1.2}$	⋯	16.3 ± 0.7	19.4 ± 2.0a
$2-2.8$	<50	${18.6}_{-1.2}^{+1.3}$	19.4 ± 1.4	12.7 ± 0.5	⋯
2–4	<50	${25.2}_{-1.5}^{+1.5}$	25.4 ± 1.6	18.6 ± 0.6	18.3 ± 1.3
4–8	<50	${1.6}_{-0.3}^{+0.4}$	⋯	3.1 ± 0.2	⋯
8–16	<50	${1.6}_{-0.3}^{+0.3}$	⋯	2.0 ± 0.2	⋯
1–2	<100	${18.1}_{-1.4}^{+1.4}$	⋯	16.3 ± 0.7b	23.0 ± 2.4a,b
$1.4-2.8$	<100	${36.5}_{-2.0}^{+2.1}$	43.1 ± 2.2	26.7 ± 0.8b	⋯
2–4	<100	${35.4}_{-2.1}^{+2.1}$	36.6 ± 2.2	23.0 ± 0.8b	23.5 ± 1.6b
4–8	<100	${3.1}_{-0.6}^{+6.7}$	⋯	4.4 ± 0.3b	⋯
8–16	<100	${2.6}_{-0.4}^{+0.5}$	⋯	2.6 ± 0.2b	⋯

Notes. Our results are FGK occurrence rates marginalized over periods down to 0.78125 days. F17 results are taken from Table 5 in Fulton et al. (2017). M15 results are taken from Table 6 in Mulders et al. (2015a), summing bins down to 0.68 days. F13 results are taken from Table 3 in Fressin et al. (2013), reported down to 0.8 days. For all results that involved summing multiple bins, we used error propagation to estimate the uncertainty.

^a $1.25\mbox{--}2$ R_⊕. ^b $P\lt 85$ days.

Download table as:  ASCIITypeset image

Table 6.  Rough Comparisons for G-type Stars

Rp (R⊕)	P (days)	This Work (%)	M15 (%)	P13 a (%)
1–2	<50	${14.1}_{-1.9}^{+2.0}$	12.9 ± 0.9	19.2 ± 1.7
2–4	<50	${33.2}_{-3.0}^{+3.0}$	20.8 ± 1.0	16.4 ± 1.2
4–8	<50	${1.8}_{-0.5}^{+0.6}$	3.3 ± 0.4	1.5 ± 0.3
8–16	<50	${1.8}_{-0.5}^{+0.5}$	1.8 ± 0.3	1.0 ± 0.4
1–2	<100	${16.0}_{-2.0}^{+2.2}$	16.0 ± 1.4b	25.0 ± 2.3
2–4	<100	${48.0}_{-3.8}^{+4.0}$	25.2 ± 1.2b	24.1 ± 1.8
4–8	<100	${4.1}_{-1.0}^{+1.2}$	4.8 ± 0.6b	2.8 ± 0.7
8–16	<100	${3.1}_{-1.8}^{+2.3}$	2.5 ± 0.4b	1.6 ± 0.5

Notes. Our results are G occurrence rates marginalized over periods down to 6.25 days (lower bound chosen to match the lower bound of the Petigura et al. 2013 results). M15 results are taken from Table 7 in Mulders et al. (2015a), summing bins down to 5.8 days. P13 results are taken from Figure 2 in Petigura et al. (2013), summing bins down to 6.25 days. For all results that involved summing multiple bins, we used error propagation to estimate uncertainty.

^aFraction of stars with planets instead of number of planets per star. ^b $P\lt 85$ days.

Download table as:  ASCIITypeset image

F17 also found a lower fraction of planets below 2 R_⊕ than older works. In particular, they found a P < 100 days, 1.41–2 R_⊕/2–2.83 R_⊕ ratio of 0.6, whereas Petigura et al. (2013) found 1.3. They explained this difference using their knowledge of a gap in the radius distribution between 1.5 and 2 R_⊕ and a peak near ∼2.5 R_⊕, which were both revealed in their study. Because these features were recovered in part due to their use of more precise stellar radii from spectroscopy, they suggested that the large (≈40%) radius uncertainties from photometry alone would scatter planets with true sizes between 2 and 2.83 R_⊕ to the 1.41–2 R_⊕ bin, both filling the gap and reducing the peak. Given that we used updated stellar radii from the Berger et al. (2018) catalog, which brought typical radius uncertainties down to ≈8% and found a similarly low ratio (0.4), the same explanation likely applies here. Thus, in combination with the results of F17, our results emphasize the sensitivity of planet occurrence rates to accurate stellar radii. We also argue that our results are more robust than previous works, given our use of updated stellar radii in combination with direct incorporation of planet radius uncertainties into our occurrence rates.

In the following sections, we discuss further comparisons to previous works in the context of interesting and informative features previously uncovered from exoplanet population analysis and occurrence rate estimates.

Marginalizing over the entire period–radius grid, we find occurrence rates of ${0.89}_{-0.16}^{+0.23}$ planets per F-type star, ${1.67}_{-0.16}^{+0.21}$ planets per G-type star, and ${2.56}_{-0.24}^{+0.29}$ per K-type star, compared to ${1.06}_{-0.07}^{+0.09}$ for the combined FGK sample.⁸ For planets within 200 days (i.e., omitting the 200–400 day bins which have low completeness and little to no planet detections below 2 R_⊕), we find occurrence rates of ${0.53}_{-0.05}^{+0.06}$ (F), ${1.17}_{-0.07}^{+0.08}$ (G), ${1.84}_{-0.13}^{+0.15}$ (K), and ${0.81}_{-0.04}^{+0.04}$ (FGK). Our results indicate a statistically significant trend of increasing planet occurrence toward cooler stars.

Figure 8 shows the entire radius distribution marginalized over P < 200 days for each stellar type. We indicate any estimates that involved marginalizing over a bin with no planet detections with a downward pointing arrow, representing an upper limit. For small planets, our results are in strong agreement with Howard et al. (2012) and Mulders et al. (2015a) that occurrence rates increase substantially toward cooler stars. Overall, for R_p = 1–2.83 R_⊕ and P < 200 days, we find occurrence rates of ${0.26}_{-0.02}^{+0.03}$ (F), ${0.67}_{-0.05}^{+0.05}$ (G), and ${1.20}_{-0.10}^{+0.11}$ (K).

Zoom In Zoom Out Reset image size

In other words, small planets around K-type stars are about twice as abundant than around G-type stars, and five times as abundant than around F-type stars. We also agree with Mulders et al. (2015b) that each distribution shows a clear drop-off in planets beyond 2.83 R_⊕. Note that some measurements on both sides of the transition are only upper limits, but this is due to a lack of planets in the 0.78–1.56 day period range, which have upper limits of less than 0.002 planets per star. The existence of the drop-off is not dependent on their contribution.

For planets beyond 2.83 R_⊕, we find that the trend becomes less clear. The 2.83–4 R_⊕ bin demonstrates no statistically significant difference between G and K occurrence rates, with a G – K difference of ${0.02}_{-0.05}^{+0.05}$ planets per star yet a K – F difference of ${0.10}_{-0.04}^{+0.05}$. Over the same radius bin, Mulders et al. (2015b) also found that G and K occurrence rates were indistinguishable while still significantly more common than around F-type stars. At larger radii where the distributions flatten out for all stellar types, we do not attempt to interpret trends, as the majority of occurrence rate estimates involve summing bins without planet detections and thus only represent upper limits.

Our FGK occurrence rate results marginalized over different period ranges are shown in Figure 9.

Zoom In Zoom Out Reset image size

First, we note that we recover the "Neptune desert," which is a dearth in Neptune- to sub-Jupiter-sized exoplanets in close-in orbits (P ≲ 3 days) that has been noted and studied in many previous works (e.g., Szabo & Kiss 2011; Mazeh et al. 2016; Owen & Lai 2018). On the lower-mass end of the desert, several studies have emphasized the role of photoevaporation on the mass loss of highly irradiated planets, while the dearth at the higher-mass end likely requires a different explanation (Ionov et al. 2018) such as tidal disruption that results from high-eccentricity migration (Owen & Lai 2018). In particular, we observe a significant decrease in occurrence rates for planets between 4 and 11.31 R_⊕ and the shortest orbital periods (0.78–3.13 days), but not larger planets, in line with expectations. Note that we did not detect any planets in the 5.66–11.31 R_⊕, 0.78–1.56 day range, nor in the 5.66–8.0 R_⊕, 1.56–3.13 day bins, so the dip is likely even lower than indicated in the plot. Meanwhile, at higher orbital periods (especially beyond 6.3 days), the distribution is flat over the same radii.

We now turn to smaller planets. An informative feature recently noted in exoplanet radius distributions is a gap between ∼1.5−2 R_⊕ for planets within P < 100 days, also known as the "radius valley" (F17; Van Eylen et al. 2018). The gap is accompanied by two peaks in radius, near ∼1.3 R_⊕ (super-Earths) and ∼2.4 R_⊕ (sub-Neptunes). This bimodal distribution had been predicted years earlier by numerical models involving the atmospheric erosion of highly irradiated low-mass planets (Lopez & Fortney 2013; Owen & Wu 2013), while its statistical significance was not established by completeness-corrected observations until F17. F17 attributed the revelation of the radius gap to their use of precise stellar radius measurements from the California-Kepler Survey (CKS; Petigura et al. 2017).

Our bin sizes were not small enough to be able to resolve this feature. This was a consequence of choosing logarithmically spaced bins large enough to be appropriate for the entire grid down to 0.5 R_⊕ and out to 400 days. Thus, we recomputed 0.78–100 day occurrence rates using the same P = {0.78, 1.56, 3.13, 6.25, 12.5, 25, 50, 100} days period bins, but much finer bins in radius space over the regime of interest: R_p = {1, 1.10, 1.21, 1.35, 1.49, 1.64, 1.81, 2, 2.21, 2.44, 2.69, 2.97, 3.28, 3.62, 4} R_⊕. The resulting distribution marginalized over the entire period range is compared with occurrence rates from Table 3 of F17, shown in the top panel of Figure 10.

Zoom In Zoom Out Reset image size

While we do find some evidence for a first peak near ∼1.4 R_⊕, a minimum at ∼1.7 R_⊕, and a second peak at ∼2.6 R_⊕, differences between our distribution and that of F17 are obvious. Planets smaller than 1.5 R_⊕ are considerably less abundant according to our sample, and we cannot confirm that our radius valley is statistically significant from these results alone. Meanwhile, our second peak is shifted to a higher radius than in F17, and it is twice as tall as the first peak rather than being comparable in height.

It should be noted that F17 only considered planets with hosts fainter than Kp = 14.2, given that the core sample of the CKS (Petigura et al. 2017) is magnitude limited, and the distribution of CKS planet radii above and below Kp = 14.2 is statistically different. In order to see if this could explain the differences between our results, we recalculated our occurrence rates using only our FGK stars with Kp < 14.2, which reduced our stellar sample from 96,280 to 30,688, and our P < 100 day planet population from 2377 to 833. This is shown in the bottom panel of Figure 10. While interpreting the results is difficult given that the occurrence rates are less well constrained by the data, we still do not recover the shape of their radius valley.

Potential explanations include our difference in occurrence rate methodology, where F17 used the IDEM and did not take into account uncertainty in planet radius. In particular, they estimated that the underlying radius distribution after removing the smear due to this uncertainty would cause the gap to become slightly deeper, but the sub-Neptune peak would be increased (see Figure 7 of F17). Furthermore, while the CKS sample represented a significant improvement in precise stellar radii over previous works, it did not yet incorporate Gaia DR2 parallaxes like the Berger et al. (2018) catalog used here. Berger et al. (2018) compared histograms of planet radii (uncorrected for occurrence rates) using their catalog and CKS-derived radii. Even under the same cuts utilized by F17, they found a higher number of sub-Neptunes in their sample (see Figure 8 of Berger et al. 2018).

We continue our analysis with these caveats and maintain our focus on overall planet occurrence rates for FGK stars (using the full 96,280 star sample, without the magnitude cut). We have the means to investigate the radius valley further, as a function of period, due to finding separate occurrence rates over specific period bins. We show these results in Figure 11. Owen & Wu (2017) revisited the radius valley following the results of F17 and developed an analytical model to demonstrate that photoevaporation should separate planets into bare cores (∼1.3 R_⊕) and those with double the core's radius (∼2.6 R_⊕). Starting with a primordial Kepler planet population and evolving the population under the effects of cooling contraction and mass loss by evaporation, Owen & Wu (2017) produced a prediction for the final radius distribution across different period ranges. Within P < 10 days, they found that super-Earths dominate, with only a small peak past 2 R_⊕. This is qualitatively similar to our P < 6.25 day occurrence rates. Not shown are our results for only P < 3.13 days, over which the sub-Neptune peak completely disappears. The 10–20 day bins in Owen & Wu (2017) demonstrated the most significant radius valley, with peaks at ∼1.3 R_⊕ and ∼2.6 R_⊕ exhibiting similar heights. Importantly, we clearly recover a similar bimodal distribution with strong peaks at ∼1.3 R_⊕ and ∼2.6 R_⊕ over roughly the same periods (6.25–25 days). Lastly, Owen & Wu (2017) found that the distribution became dominated by sub-Neptunes at larger orbital periods, with only a small super-Earth peak over 20–40 days and no such peak over 40–100 days. We find that the super-Earth peak disappears earlier, with no such peak over 25–100 days. Overall, our occurrence rates provide strong observational evidence in support of the Owen & Wu (2017) model and can inform future studies on theoretical explanations for the radius valley.

Zoom In Zoom Out Reset image size

Our FGK occurrence rate results marginalized over the different radius ranges are shown in Figure 12.

Zoom In Zoom Out Reset image size

The 1–2 R_⊕ and 2–4 R_⊕ distributions show clear increases in ${df}/d\mathrm{log}P$ with P up to a transition period P₀, followed by a decreasing trend for 1–2 R_⊕ and a flatter trend for 2–4 R_⊕. P₀ could indicate an important orbital distance down to which migration deposits planets, and differences between these distributions could indicate such mechanisms depend on planet size. We fit the function

$$\begin{eqnarray}&&F=\displaystyle \frac{{df}}{d\mathrm{log}P}={{CP}}^{\beta }\left(1-{e}^{-{\left(P/{P}_{0}\right)}^{\gamma }}\right)\end{eqnarray} \tag{ 22 }$$

from Howard et al. (2012) to these distributions using the extension of their maximum-likelihood method outlined in Petigura et al. (2018). For the ith bin, the log-likelihood of the model is

$$\begin{eqnarray}&&\mathrm{ln}{L}_{i}={n}_{\mathrm{pl},i}\mathrm{ln}F{\rm{\Delta }}{\boldsymbol{x}}+{n}_{\mathrm{nd},i}\mathrm{ln}(1-F{\rm{\Delta }}{\boldsymbol{x}}),\end{eqnarray} \tag{ 23 }$$

where ${\rm{\Delta }}{\boldsymbol{x}}={\rm{\Delta }}\mathrm{log}P\times {\rm{\Delta }}\mathrm{log}{R}_{p}$ is the size of the bin, ${n}_{\mathrm{pl},i}$ is the number of planets detected in the bin, and ${n}_{\mathrm{nd},i}\,={n}_{\mathrm{pl},i}/{f}_{i}-{n}_{\mathrm{pl},i}$ is the effective number of nondetections as estimated using the bin's occurrence rate f_i. The maximum-likelihood solution is obtained by maximizing the combined log-likelihood over all bins:

$$\begin{eqnarray}&&\mathrm{ln}L=\sum _{i=1}^{{n}_{\mathrm{bin}}}{L}_{i}.\end{eqnarray} \tag{ 24 }$$

We used MCMC sampling with emcee to explore the parameter space. The median and 68.3% credible interval for each distribution is shown in Equation (22), with associated parameters in Table 7.

Table 7.  Median and 68.3% Credible Interval Parameters for Equation (22), Describing the Shape of the Occurrence Rate Distribution with Orbital Period over Different Size Ranges

Rp (R⊕)	C	β	γ	P0 (days)
1–2	${0.36}_{-0.10}^{+0.07}$	$-{0.5}_{-0.1}^{+0.1}$	${2.42}_{-0.1}^{+0.1}$	${5.9}_{-0.5}^{+0.5}$
2–4	${0.33}_{-0.12}^{+0.08}$	$-{0.1}_{-0.1}^{+0.1}$	${2.3}_{-0.1}^{+0.1}$	${13.3}_{-1.5}^{+1.4}$

Download table as:  ASCIITypeset image

This function simplifies to two power laws far from P₀:

$$\begin{eqnarray}\displaystyle \frac{{df}}{d\mathrm{log}P}\propto \left\{\begin{array}{ll}{P}^{\alpha } & \mathrm{if}\ P\ll {P}_{0}\ ,\ \mathrm{where}\ \alpha =\beta +\gamma \\ {P}^{\beta } & \mathrm{if}\ P\gg {P}_{0}.\end{array}\right.\end{eqnarray} \tag{ 25 }$$

The 1–2 R_⊕ planet occurrence rate rises with P, with ${df}\propto {P}^{\alpha }d\mathrm{log}P$, where $\alpha ={1.9}_{-0.1}^{+0.1}$, up to a transition period ${P}_{0}={5.9}_{-0.5}^{+0.5}$ days. Beyond this, ${df}\propto {P}^{\beta }d\mathrm{log}P$, where $\beta \,=-{0.5}_{-0.1}^{+0.1}$. Comparatively, Petigura et al. (2018) looked at 1–1.7 R_⊕ bins and found a slightly higher initial increase ($\alpha ={2.4}_{-0.3}^{+0.4}$), a slightly higher transition period (${P}_{0}={6.5}_{-1.2}^{+1.6}$ days), and a slightly shallower decrease at longer orbital periods ($\beta =-{0.3}_{-0.2}^{+0.2}$), though the 68.3% credible intervals of the latter two parameters overlap with ours. Meanwhile, the long-period distribution found by Dong & Zhu (2013) was flat, with β = −0.10 ± 0.12.

The transition for 2–4 R_⊕ occurrence rates occurs farther out, at ${13.3}_{-1.5}^{+1.4}$ days, with a similar rapidly rising distribution at shorter orbital periods ($\alpha ={2.2}_{-0.1}^{+0.1}$) and a nearly flat distribution at longer orbital periods ($\beta =-{0.1}_{-0.1}^{+0.1}$). These are consistent with a ${P}_{0}={11.9}_{-1.5}^{+1.7}$ day transition period, $\alpha \,={2.3}_{-0.2}^{+0.2}$, and $\beta =-{0.1}_{-0.1}^{+0.1}$ for $1.7-4$ R_⊕ planets from Petigura et al. (2018). Dong & Zhu (2013) also found a nearly flat distribution for 2–4 R_⊕ with $\beta =0.11\pm 0.05$. However, while all three studies agree that ∼1–2 R_⊕ planets are more common than ∼2–4 R_⊕ planets before the small-planet transition, our results and those of Petigura et al. (2018) would indicate the opposite is true past the transition while Dong & Zhu (2013) found similar occurrence rates for both distributions. We believe this is due to our use of more up-to-date and precise stellar radii, causing many 1–2 R_⊕ planets to be pushed into the 2–4 R_⊕ bin.

All three studies indicate that the distributions of larger planets (here, 4–8 R_⊕ and 8–16 R_⊕) are inconsistent with this power-law cut-off model, with occurrence rates gradually increasing over the entire period range. Our 4–8 R_⊕ occurrence rates do jump suddenly at 3.1 days, though this is likely another look at the Neptune Desert described in Section 6.3.

For 8–16 R_⊕, we do no confirm the three-day "pile-up" of hot Jupiters present in radial velocity (RV) surveys (e.g., Cumming et al. 1999; Udry et al. 2003; Wright et al. 2009). This pile-up features strongly in various high-eccentricity migration scenarios (e.g., Fabrycky & Tremaine 2007; Wu et al. 2007; Wu & Lithwick 2011). However, other Kepler-based studies have called into question the pile-up (Howard et al. 2012; Fressin et al. 2013), and differences in overall hot Jupiter (P ≲ 10 days) occurrence rates between Kepler and RV surveys have been previously noted. In particular, Kepler hot Jupiter occurrence rates typically lie at around 0.4%–0.6% (e.g., Howard et al. 2012; Fressin et al. 2013; Mulders et al. 2015a; Petigura et al. 2018) with RV occurrence rates are at around 0.9%–1.2% (e.g., Marcy et al. 2005; Mayor et al. 2011). Our own hot Jupiter estimate should be intermediate between ${0.43}_{-0.09}^{+0.10}$% ($0.78-6.25$ days) and ${0.77}_{-0.14}^{+0.16}$% (0.78–12.5 days), consistent with other Kepler results. Dawson & Murray-Clay (2013) suggested that the pile-up was a feature of metal-rich stars ([Fe/H] ≥ 0) specifically, while the Kepler sample has systematically lower metallicity than RV samples. They largely recovered the pile-up in the Kepler sample when considering only stars with super-solar metallicity. Giant planet occurrence has also been shown to correlate strongly with host-star metallicity (Santos et al. 2003; Fischer & Valenti 2005; Petigura et al. 2018). An assessment of the presence of the pile-up in our planet catalog under these conditions will be left to a future paper focused on planet occurrence and its dependence on stellar metallicity.

While our baseline results incorporate catalog completeness, our planet sample is also not expected to be completely "reliable." Signals not caused by planet transits, whether transit-like noise or astrophysical FPs, may be erroneously classified as planets by the vetting pipeline. The corresponding concept of reliability refers to the fraction of planets in the catalog that are actually planets. We focus on our catalog's reliability against noise specifically, which is the largest concern for candidates near our S/N = 6, three-transit detection limit, including small, rocky planets in orbits with long periods.

The incorporation of reliability against noise into occurrence rate estimates remains an open question and was only first directly tackled in Bryson et al. (2019). Bryson et al. (2019) took a probabilistic approach to reliability models, fitting components of the reliability with functions over finely spaced bins in period–S/N space (where S/N is represented by the MES employed by the Kepler pipeline). Rather than assume a functional form of the reliability, we take a simplified approach using our reliability results from Paper I, summarized below.

We simulated noise using two data sets from the original 198,640 light curves searched. First, we recreated the inverted (INV) set described in Christiansen (2017) for their own reliability tests by inverting the light curves. Second, we recreated the Scrambled Group 1 (SCR1) set by reordering the Kepler quarters according to the first order described in Coughlin (2017). Each of these data sets allowed for the existence of realistic signals with noise properties similar to the real data while removing the possibility of planet transit detection. After searching and vetting these data, we estimated the fraction of noise FPs successfully classified as FPs (the "effectiveness" of the pipeline, E) to be 99.9% overall, having passed only 36 of 27,386 noise FPs as PCs. Using the definition of reliability from Thompson et al. (2018),

$$\begin{eqnarray}&&R=1-\displaystyle \frac{{N}_{\mathrm{FP}}}{{N}_{\mathrm{PC}}}\left(\displaystyle \frac{1-E}{E}\right),\end{eqnarray} \tag{ 26 }$$

where N_PC and N_FP are the number of observed PCs and FPs identified by the vetting pipeline, our pipeline has an overall reliability against noise FPs of 98.3%.

To set up the application of these results to our occurrence rate estimates, we determined reliability over a coarse grid in period–S/N space in an attempt to reduce the effect of small number statistics. We chose period bin edges equal to those used for our occurrence rates and make the assumption that the reliability of each bin is constant across that bin. While ideally we would find our reliability across only the FGK stars in our sample, a concern is that small number statistics would be more significant than for the full, ∼200,000-star sample. Thus, while we recognize that noise properties between the two samples should be different, we elected to use our results across all stars (top panel of Figure 13) to improve the S/N of the estimate. Our FGK-only results are included in the bottom panel of Figure 13 for comparison.

Zoom In Zoom Out Reset image size

Lastly, to incorporate reliability into the ABC methodology itself, we recall the discussion of selection effects from Youdin (2011) as outlined in Section 5.2. The FP-related selection effect was ${\eta }_{\mathrm{fp}}=1/(1-{r}_{\mathrm{fp}})$, where ${r}_{\mathrm{fp}}$ is the rate of FP events that are detected as planets. This is equal to 1 minus the reliability, giving ${\eta }_{\mathrm{fp}}=1/R$. Thus, we replace Equation (18) in our forward model with

$$\begin{eqnarray}&&{\eta }_{\mathrm{tot}}={\eta }_{\mathrm{tr}}{\eta }_{\mathrm{rec}}/R,\end{eqnarray} \tag{ 27 }$$

where R is found for a given planet according to its period and S/N.

We recalculated our FGK occurrence rates under these changes, with median and 68.3% credible interval results included in Table A1 alongside our baseline estimates. Reliability did not significantly impact our results, with posteriors for every cell overlapping significantly. This is not unexpected, given that our reliability is high ($\gt 90$%) outside only the long-period, low-S/N corner in period–S/N space. For the corresponding period–radius cells, it is difficult to assess the full impact of reliability given that the same areas have very low completeness and no planet detections, meaning estimates are already not as well constrained as other areas. However, we do tend to reduce our upper limits in the most affected areas. In particular, our 0.5–0.71 R_⊕, 200–400 day cell had its upper limit reduced from 14.4% to 13.1%, and our 0.71–1 R_⊕, 200–400 day cell had its upper limit reduced from 3.4% to 2.7%.

Given that we focused in our previous sections on high-reliability regimes (i.e., mainly R_p > 1 R_⊕, P < 200 days), we do not expect the lack of reliability incorporation to affect our previous analysis. However, for our upcoming terrestrial HZ planet frequency discussion, which will specifically depend on calculations over regions of low reliability, we will report both versions of occurrence rate results.

Our direct calculation over the 237–500 day bin represents a subset of the 237–860 day optimistic HZ. For the 1–1.5 R_⊕ range, we find an occurrence rate of ${0.05}_{-0.03}^{+0.04}$ planets per star. When incorporating reliability, we find an occurrence rate of ${0.05}_{-0.02}^{+0.03}$ planets per star. For the 0.75–1.5 R_⊕ range considered by Hsu et al. (2019), we find an occurrence rate of ${0.12}_{-0.06}^{+0.08}$ (${0.10}_{-0.04}^{+0.07}$ with reliability), which overlaps well with the ${0.16}_{-0.06}^{+0.11}$ estimate from Hsu et al. (2019). Lastly, we find a considerably increased estimate with substantial uncertainties for 0.5–1.5 R_⊕, at ${0.38}_{-0.19}^{+0.29}$ (${0.31}_{-0.15}^{+0.28}$), on account of the low completeness and poor constraints provided by the data. Note that neither of our studies had planet detections over any of these radius ranges, so these occurrence rates are best interpreted as upper limits.

Considering the entire 237–860 day HZ range requires extrapolating these results to longer periods. Interpreting results obtained via extrapolation should be done with added caution, considering we have demonstrated substantial uncertainties in sub-${\eta }_{\oplus }$ occurrence rate estimates and will necessarily have to make an assumption about the nature of planet distributions beyond the limit of our sensitivity. However, the ABC methodology requiring such extrapolation is not unique. All studies are limited by the Kepler mission duration to planets within ≈500 days, and small planets typically require more observed transits than larger planets in order to produce signals with sufficient S/N for detection. Other grid-based occurrence rates must make similar assumptions to our work. Meanwhile, studies that find a function to describe planet distributions with period and/or radius may integrate over a desired ${\eta }_{\oplus }$ range to produce an estimate, but the function itself will have been based on planets with larger sizes and/or shorter orbital periods. In these cases, the employed assumption is that the same model can also explain the ${\eta }_{\oplus }$ regime, which is not necessarily true; Mulders et al. (2018), for instance, found that their broken power-law model broke down outside of 400 days. Furthermore, we have shown that there are numerous period- and size-dependent small-scale variations (such as the clear radius valley for orbital periods within 25 days), indicating that a parametric occurrence rate model is not necessarily the best descriptor of the data even over the shorter orbital periods and larger planet sizes considered by these works.

With these caveats, we adopt the method of extrapolation used by Hsu et al. (2019) and assume that the differential occurrence rate derived over 237–500 days and a given radius range is the same over longer periods. Under this assumption, we estimate optimistic HZ occurrence rates of ${\eta }_{\oplus }={0.08}_{-0.04}^{+0.07}$ (${0.08}_{-0.04}^{+0.06}$) planets per star for 1–1.5 R_⊕, ${\eta }_{\oplus }={0.21}_{-0.10}^{+0.14}$ (${0.17}_{-0.08}^{+0.11}$) for 0.75–1.5 R_⊕, and ${\eta }_{\oplus }={0.66}_{-0.32}^{+0.51}$ (${0.53}_{-0.26}^{+0.48}$) for 0.5–1.5 R_⊕.

The 0.99–1.70 au conservative HZ from Kopparapu et al. (2013) corresponds to orbital periods of 360–809 days. Extrapolating our 237–500 day results over these periods, we find occurrence rates of ${\eta }_{\oplus }$=${0.05}_{-0.03}^{+0.04}$ (${0.05}_{-0.02}^{+0.04}$) planets per star for 1–1.5 R_⊕, ${\eta }_{\oplus }$=${0.13}_{-0.06}^{+0.09}$ (${0.11}_{-0.05}^{+0.07}$) for 0.75–1.5 R_⊕, and ${\eta }_{\oplus }$ = ${0.42}_{-0.20}^{+0.32}$ (${0.34}_{-0.16}^{+0.30}$) for 0.5–1.5 R_⊕.

The challenges involved with defining the bounds of ${\eta }_{\oplus }$ have motivated recent studies to instead report and compare ${{\rm{\Gamma }}}_{\oplus }$, the differential occurrence rate near the HZ (Youdin 2011; Foreman-Mackey et al. 2014; Burke et al. 2015). We follow Hsu et al. (2019) in defining ${{\rm{\Gamma }}}_{\oplus }$ using our 237–500 day, $0.75\,-1.5$ R_⊕ results, giving ${{\rm{\Gamma }}}_{\oplus }\equiv ({d}^{2}f)/[d(\mathrm{ln}P)d(\mathrm{ln}{R}_{p})]\,={0.23}_{-0.11}^{+0.16}({0.19}_{-0.08}^{+0.13})$.

Comparisons with other ${{\rm{\Gamma }}}_{\oplus }$ estimates in the literature are shown in Figure 14. The values from Pascucci et al. (2019) correspond to their Model #4 and #6 results given in their Table 2, which address the impact of photoevaporated cores by excluding planets within 12 and 25 days from their analysis, respectively. The ExoPAG SAG13 estimate, obtained via Kopparapu et al. (2018), is based on a meta-analysis of community-submitted G-type star occurrence rate studies, for which a broken power law was fit to a combined period–radius grid. The plotted central values were found by plugging in Earth's radius and period into their baseline power law, while the lower and upper uncertainties correspond to their pessimistic and optimistic power laws, respectively. The values from Burke et al. (2015) represent the allowable range from their sensitivity analysis. The Petigura et al. (2013) result is calculated by converting their 200–400 day, 1–2 R_⊕ extrapolated occurrence rate into a differential rate. The Dong & Zhu (2013) value is an extrapolation of their 1–2 R_⊕ best-fit function, evaluated at the orbital period of Earth and with error bars given by the propagation of uncertainty. All other values are the ${{\rm{\Gamma }}}_{\oplus }$ results explicitly reported in their respective papers.

Zoom In Zoom Out Reset image size

Our results compare most favorably to those from Pascucci et al. (2019), Hsu et al. (2019), Bryson et al. (2019), Kopparapu et al. (2018), and Dong & Zhu (2013), and are well within the allowable range of Burke et al. (2015). Meanwhile, the values from Youdin (2011), Garret et al. (2018), and Zink et al. (2019) all lie above our upper limits. Lack of consistency with the early Youdin (2011) result can readily be explained by the fact that it is based on a much older Kepler catalog, having included only planets with $P\lt 50$ days. Part of the disagreement with Zink et al. (2019) can be explained by their incorporation of transit multiplicity—in other words, they took into account a reduction in detection efficiency with planet detection order, which would result in an increased occurrence rate—though this is unlikely to explain the whole discrepancy (see Section 8). Nevertheless, we argue that our uncertainty estimates are more realistic than those of Garret et al. (2018) and Zink et al. (2019).

Turning to the impact of reliability, we again do not find significantly different occurrence rates. This is contrasted to the results of Bryson et al. (2019), who found that incorporating reliability dropped ${{\rm{\Gamma }}}_{\oplus }$ from ${0.205}_{-0.073}^{+0.106}$ down to ${0.083}_{-0.036}^{+0.050}$. We note three things: first, we reemphasize that given the total lack of planet detections in the relevant bins, our results are poorly constrained by the data and should be interpreted as upper limits. Had we been able to place better constraints on these occurrence rates, a clearer picture of the full impact of reliability would appear. Second, we only took into account reliability against noise and systematics, whereas Bryson et al. (2019) additionally considered reliability against astrophysical FPs. Their false-alarm-only reliability was ${{\rm{\Gamma }}}_{\oplus }={0.128}_{-0.049}^{+0.077}$. Lastly, our grid-based occurrence rate method means that the reliability of a given bin will only directly affect the occurrence rates of that bin, and only indirectly affect the occurrence rate of neighboring, simultaneously fit bins due to the leaking of planets between bins. Bryson et al. (2019) fit a joint power-law model across all planets with ${R}_{p}=0.75\mbox{--}2.5$ R_⊕ and $P=50\mbox{--}400$ days, meaning that the reliability of any planet would affect the fit results across the entire space considered. We would expect this to cause the reliability to have a greater affect on occurrence rates than in our work.

A subset of ${\eta }_{\oplus }$ not yet mentioned is the concept of ${\zeta }_{\oplus }$: the number of planets per star with radii and orbital periods within 20% of Earth's values, as introduced by Burke et al. (2015). Combining both an extrapolated analysis (assuming their 0.75–2.5 R_⊕, 50–300 day broken power law) and a direct analysis (recomputing a broken power law over 300–700 days), Burke et al. (2015) reported ${\zeta }_{\oplus }=0.10$ planets per star with an allowable range of 0.01–2. To compare to this value, we recalculated our occurrence rates with bin edges of R_p = {0.8, 1.0, 1.2} R_⊕ and P = 292–438 days. We recovered their central value almost exactly, finding ${\zeta }_{\oplus }={0.12}_{-0.06}^{+0.09}$ (${0.10}_{-0.06}^{+0.07}$).

For the definition of the HZ, Kopparapu et al. (2013) recommended the use of their conservative (0.99–1.70 au) boundaries. For the lower radius limit, we are inclined to consider potentially habitable rocky planets as those down to 0.75 R_⊕ rather than 0.5 R_⊕, as this is near the scientifically motivated 0.72 R_⊕ limit used by Zink & Hansen (2019). For the upper radius limit, the 1.5 R_⊕ radius considered throughout this section is already well motivated by both the characteristics of the radius valley and the findings of Rogers (2015). Meanwhile, incorporating reliability allows for better constraints on ${\eta }_{\oplus }$ upper limits.

For our reliability-incorporated, 0.75–1.5 R_⊕ conservative HZ, our {5, 15.9, 50, 84.1, 95}th percentiles are ${\eta }_{\oplus }=\{0.04,0.06,0.11,0.18,0.24\}$ planets per star. Thus, we suggest future exoplanet characterization and habitability-related missions to consider an upper limit (84.1th percentile) of <0.18 terrestrial HZ planets per Sun-like star.

Many of the limitations discussed in this section (reliability against astrophysical FPs, nonzero eccentricity, use of a Dirichlet prior, use of DR25 window functions) would lead to reduced occurrence rates should we apply them, so this upper limit should be robust to these concerns. However, transit multiplicity would indicate that our reported occurrence rates are underestimated. Hsu et al. (2019) estimated a less than ∼8% increase in their DR25 occurrence rates for long-period planets due to this effect. Our pipeline may be slightly more affected given its lower recovery efficiency for planets with few transits. Thus, we estimate that our ${\eta }_{\oplus }$ upper limit is robust to ∼10%.