Occurrence Rates of Planets orbiting FGK Stars: Combining Kepler DR25, Gaia DR2 and Bayesian Inference

TL;DR

This study estimates the occurrence rates of planets around FGK stars using Kepler and Gaia data, employing Bayesian methods to provide robust upper limits and inform future mission designs for rocky planets in habitable zones.

Contribution

It combines Kepler DR25, Gaia DR2, and Bayesian inference to accurately estimate planet occurrence rates across a wide range of sizes and periods, including near the habitable zone.

Findings

Estimated occurrence rate for 0.75-1.5 R⊕ planets with 237-500 day periods is less than 0.27 per star.

Upper limits on occurrence rates are robust within approximately 10%.

Provides recommended ranges for future mission planning targeting rocky planets.

Abstract

We characterize the occurrence rate of planets, ranging in size from 0.5-16 R$_\oplus$, orbiting FGK stars with orbital periods from 0.5-500 days. Our analysis is based on results from the "DR25" catalog of planet candidates produced by NASA's Kepler mission and stellar radii from Gaia "DR2". We incorporate additional Kepler data products to accurately characterize the efficiency of planets being recognized as a "threshold crossing events" (TCE) by Kepler's Transiting Planet Search pipeline and labeled as a planet candidate by the robovetter. Using a hierarchical Bayesian model, we derive planet occurrence rates for a wide range of planet sizes and orbital periods. For planets with sizes $0.75-1.5$ R$_\oplus$ and orbital periods of 237-500 days, we find a rate of planets per FGK star of $<0.27$ ($84.13$th percentile). While the true rate of such planets could be lower by a factor of $\sim~2$ (primarily due to potential contamination of planet candidates by false alarms), the upper limits on the occurrence rate of such planets are robust to $\sim~10\%$. We recommend that mission concepts aiming to characterize potentially rocky planets in or near the habitable zone of sun-like stars prepare compelling science programs that would be robust for a true rate in the range $f_{R,P} = $ $0.03-0.40$ for $0.75-1.5$ R$_\oplus$ planets with orbital periods in 237-500 days, or a differential rate of $\Gamma_\oplus \equiv (d^2 f)/[d(\ln P)~d(\ln R_{p})] = $ $0.06-0.76$.

Figures

Tables

Table 1. Table 1: Parameters for p det & vet subscript 𝑝 det vet p_{\mathrm{det\&vet}}

$N_{tr}$	$\alpha$	$\beta$	$c$
3	33.3884	0.264472	0.699093
4	32.8860	0.269577	0.768366
5	31.5196	0.282741	0.833673
6	30.9919	0.286979	0.859865
7-9	30.1906	0.294688	0.875042
10-18	31.6342	0.279425	0.886144
19-36	32.6448	0.268898	0.889724
$\ge 37$	27.8185	0.32432	0.945075

Table 2. Table 2: Inferred Planet Occurrence Rates

Period	Radius	Combined Detection &	No Vetting
(days)	(R⊕)	Vetting Efficiency	Efficiency
$0.50-1.00$	$0.50-0.75$	${2.5}_{-1.7}^{+2.6}\times {10}^{-3}$	${1.71}_{-0.58}^{+1.34}\times {10}^{-3}$
$0.50-1.00$	$0.75-1.00$	${2.70}_{-0.92}^{+0.80}\times {10}^{-3}$	${2.5}_{-1.0}^{+1.2}\times {10}^{-3}$
$0.50-1.00$	$1.00-1.25$	${1.7}_{-1.1}^{+1.2}\times {10}^{-3}$	${1.93}_{-0.74}^{+1.15}\times {10}^{-3}$
$0.50-1.00$	$1.25-1.50$	${1.57}_{-0.69}^{+0.74}\times {10}^{-3}$	${1.73}_{-1.07}^{+0.87}\times {10}^{-3}$
$0.50-1.00$	$1.50-1.75$	${8.3}_{-5.3}^{+6.0}\times {10}^{-4}$	${7.3}_{-4.0}^{+5.3}\times {10}^{-4}$
$0.50-1.00$	$1.75-2.00$	${2.7}_{-1.9}^{+2.6}\times {10}^{-4}$	${4.1}_{-2.1}^{+4.3}\times {10}^{-4}$
$0.50-1.00$	$2.00-2.50$
$0.50-1.00$	$2.50-3.00$
$0.50-1.00$	$3.00-4.00$
$0.50-1.00$	$4.00-6.00$	${2.9}_{-2.0}^{+1.4}\times {10}^{-4}$	${3.2}_{-1.3}^{+1.4}\times {10}^{-4}$
$0.50-1.00$	$6.00-8.00$
$0.50-1.00$	$8.00-12.00$
$0.50-1.00$	$12.00-16.00$
$1.00-2.00$	$0.50-0.75$	${2.6}_{-1.7}^{+1.4}\times {10}^{-3}$	${2.7}_{-1.5}^{+2.2}\times {10}^{-3}$
$1.00-2.00$	$0.75-1.00$	${2.9}_{-1.1}^{+1.6}\times {10}^{-3}$	${2.49}_{-1.77}^{+0.93}\times {10}^{-3}$
$1.00-2.00$	$1.00-1.25$	${1.03}_{-0.68}^{+0.96}\times {10}^{-3}$	${1.02}_{-0.65}^{+0.65}\times {10}^{-3}$
$1.00-2.00$	$1.25-1.50$	${2.2}_{-1.3}^{+2.0}\times {10}^{-3}$	${2.7}_{-1.2}^{+3.1}\times {10}^{-3}$
$1.00-2.00$	$1.50-1.75$	${9.7}_{-8.3}^{+6.7}\times {10}^{-4}$	${6.4}_{-3.4}^{+8.4}\times {10}^{-4}$
$1.00-2.00$	$1.75-2.00$	${6.4}_{-3.3}^{+4.9}\times {10}^{-4}$	${6.1}_{-3.7}^{+3.3}\times {10}^{-4}$
$1.00-2.00$	$2.00-2.50$	${2.3}_{-1.5}^{+2.4}\times {10}^{-4}$	${2.9}_{-1.4}^{+2.5}\times {10}^{-4}$
$1.00-2.00$	$2.50-3.00$	${2.7}_{-2.1}^{+3.6}\times {10}^{-4}$	${2.2}_{-2.1}^{+4.1}\times {10}^{-4}$
$1.00-2.00$	$3.00-4.00$
$1.00-2.00$	$4.00-6.00$
$1.00-2.00$	$6.00-8.00$
$1.00-2.00$	$8.00-12.00$	${2.5}_{-1.0}^{+2.7}\times {10}^{-4}$	${2.9}_{-1.8}^{+1.5}\times {10}^{-4}$
$1.00-2.00$	$12.00-16.00$	${4.1}_{-1.2}^{+6.9}\times {10}^{-4}$	${4.4}_{-3.1}^{+2.6}\times {10}^{-4}$
$2.00-4.00$	$0.50-0.75$	${2.65}_{-0.77}^{+0.86}\times {10}^{-2}$	${2.49}_{-0.75}^{+0.62}\times {10}^{-2}$
$2.00-4.00$	$0.75-1.00$	${1.20}_{-0.45}^{+0.15}\times {10}^{-2}$	${10.0}_{-1.7}^{+3.7}\times {10}^{-3}$
$2.00-4.00$	$1.00-1.25$	${2.5}_{-1.3}^{+1.9}\times {10}^{-3}$	${2.3}_{-1.4}^{+1.0}\times {10}^{-3}$
$2.00-4.00$	$1.25-1.50$	${9.2}_{-2.3}^{+2.8}\times {10}^{-3}$	${9.4}_{-2.3}^{+2.8}\times {10}^{-3}$
$2.00-4.00$	$1.50-1.75$	${5.6}_{-2.8}^{+2.5}\times {10}^{-3}$	${5.4}_{-2.4}^{+3.1}\times {10}^{-3}$
$2.00-4.00$	$1.75-2.00$	${2.9}_{-1.2}^{+1.9}\times {10}^{-3}$	${2.9}_{-1.5}^{+5.8}\times {10}^{-3}$
$2.00-4.00$	$2.00-2.50$	${2.46}_{-1.19}^{+0.74}\times {10}^{-3}$	${2.5}_{-1.6}^{+1.1}\times {10}^{-3}$
$2.00-4.00$	$2.50-3.00$	${7.5}_{-4.7}^{+5.1}\times {10}^{-4}$	${9.8}_{-5.8}^{+8.7}\times {10}^{-4}$
$2.00-4.00$	$3.00-4.00$	${1.28}_{-0.63}^{+0.53}\times {10}^{-3}$	${1.21}_{-0.61}^{+1.04}\times {10}^{-3}$
$2.00-4.00$	$4.00-6.00$	${6.1}_{-4.1}^{+4.7}\times {10}^{-4}$	${7.1}_{-4.1}^{+4.6}\times {10}^{-4}$
$2.00-4.00$	$6.00-8.00$
$2.00-4.00$	$8.00-12.00$	${9.0}_{-3.6}^{+8.2}\times {10}^{-4}$	${8.4}_{-5.1}^{+7.6}\times {10}^{-4}$
$2.00-4.00$	$12.00-16.00$	${1.62}_{-0.54}^{+0.87}\times {10}^{-3}$	${1.84}_{-0.82}^{+0.65}\times {10}^{-3}$
$4.00-8.00$	$0.50-0.75$	${9.2}_{-2.9}^{+2.2}\times {10}^{-2}$	${7.6}_{-2.6}^{+3.1}\times {10}^{-2}$
$4.00-8.00$	$0.75-1.00$	${3.18}_{-0.73}^{+1.32}\times {10}^{-2}$	${2.90}_{-0.97}^{+0.77}\times {10}^{-2}$
$4.00-8.00$	$1.00-1.25$	${2.40}_{-0.45}^{+0.70}\times {10}^{-2}$	${2.26}_{-0.52}^{+0.66}\times {10}^{-2}$
$4.00-8.00$	$1.25-1.50$	${1.91}_{-0.50}^{+0.43}\times {10}^{-2}$	${1.93}_{-0.48}^{+0.70}\times {10}^{-2}$
$4.00-8.00$	$1.50-1.75$	${1.40}_{-0.40}^{+0.44}\times {10}^{-2}$	${1.38}_{-0.41}^{+0.40}\times {10}^{-2}$
$4.00-8.00$	$1.75-2.00$	${6.9}_{-3.3}^{+4.7}\times {10}^{-3}$	${7.0}_{-3.2}^{+5.2}\times {10}^{-3}$
$4.00-8.00$	$2.00-2.50$	${1.03}_{-0.28}^{+0.39}\times {10}^{-2}$	${1.07}_{-0.32}^{+0.30}\times {10}^{-2}$
$4.00-8.00$	$2.50-3.00$	${1.10}_{-0.30}^{+0.28}\times {10}^{-2}$	${1.10}_{-0.30}^{+0.27}\times {10}^{-2}$
$4.00-8.00$	$3.00-4.00$	${3.2}_{-1.8}^{+1.9}\times {10}^{-3}$	${3.6}_{-1.8}^{+2.2}\times {10}^{-3}$
$4.00-8.00$	$4.00-6.00$	${1.54}_{-1.02}^{+0.80}\times {10}^{-3}$	${1.49}_{-0.74}^{+0.76}\times {10}^{-3}$
$4.00-8.00$	$6.00-8.00$	${6.1}_{-3.7}^{+4.0}\times {10}^{-4}$	${5.7}_{-4.2}^{+5.6}\times {10}^{-4}$
$4.00-8.00$	$8.00-12.00$	${1.17}_{-0.54}^{+0.87}\times {10}^{-3}$	${1.19}_{-0.60}^{+0.84}\times {10}^{-3}$
$4.00-8.00$	$12.00-16.00$	${1.65}_{-0.59}^{+1.05}\times {10}^{-3}$	${1.75}_{-0.71}^{+0.84}\times {10}^{-3}$
$8.00-16.00$	$0.50-0.75$	${6.6}_{-3.3}^{+5.9}\times {10}^{-2}$	${5.4}_{-2.9}^{+3.5}\times {10}^{-2}$
$8.00-16.00$	$0.75-1.00$	${5.2}_{-1.3}^{+1.4}\times {10}^{-2}$	${4.81}_{-0.77}^{+1.33}\times {10}^{-2}$
$8.00-16.00$	$1.00-1.25$	${3.36}_{-0.94}^{+1.18}\times {10}^{-2}$	${2.90}_{-0.48}^{+0.79}\times {10}^{-2}$
$8.00-16.00$	$1.25-1.50$	${3.56}_{-0.74}^{+0.84}\times {10}^{-2}$	${3.25}_{-0.92}^{+0.78}\times {10}^{-2}$
$8.00-16.00$	$1.50-1.75$	${1.71}_{-0.62}^{+0.91}\times {10}^{-2}$	${1.64}_{-0.55}^{+0.78}\times {10}^{-2}$
$8.00-16.00$	$1.75-2.00$	${1.29}_{-0.62}^{+0.53}\times {10}^{-2}$	${1.29}_{-0.59}^{+0.72}\times {10}^{-2}$
$8.00-16.00$	$2.00-2.50$	${3.09}_{-0.80}^{+0.82}\times {10}^{-2}$	${3.18}_{-0.56}^{+0.53}\times {10}^{-2}$
$8.00-16.00$	$2.50-3.00$	${2.65}_{-0.45}^{+0.43}\times {10}^{-2}$	${2.63}_{-0.53}^{+0.48}\times {10}^{-2}$
$8.00-16.00$	$3.00-4.00$	${7.0}_{-3.2}^{+3.8}\times {10}^{-3}$	${7.2}_{-3.8}^{+3.5}\times {10}^{-3}$
$8.00-16.00$	$4.00-6.00$	${3.5}_{-1.7}^{+1.5}\times {10}^{-3}$	${3.9}_{-1.7}^{+2.5}\times {10}^{-3}$
$8.00-16.00$	$6.00-8.00$	${1.07}_{-0.70}^{+1.12}\times {10}^{-3}$	${1.096}_{-0.450}^{+0.095}\times {10}^{-3}$
$8.00-16.00$	$8.00-12.00$	${1.58}_{-0.73}^{+0.74}\times {10}^{-3}$	${1.45}_{-0.60}^{+1.59}\times {10}^{-3}$
$8.00-16.00$	$12.00-16.00$	${1.02}_{-0.66}^{+0.72}\times {10}^{-3}$	${9.7}_{-6.2}^{+6.4}\times {10}^{-4}$
$16.00-32.00$	$0.50-0.75$
$16.00-32.00$	$0.75-1.00$	${8.7}_{-2.9}^{+2.5}\times {10}^{-2}$	${7.2}_{-1.4}^{+2.0}\times {10}^{-2}$
$16.00-32.00$	$1.00-1.25$	${3.4}_{-1.6}^{+1.5}\times {10}^{-2}$	${3.0}_{-1.4}^{+1.7}\times {10}^{-2}$
$16.00-32.00$	$1.25-1.50$	${2.47}_{-0.83}^{+2.57}\times {10}^{-2}$	${2.52}_{-1.00}^{+0.74}\times {10}^{-2}$
$16.00-32.00$	$1.50-1.75$	${2.40}_{-0.68}^{+0.82}\times {10}^{-2}$	${2.28}_{-0.26}^{+0.30}\times {10}^{-2}$
$16.00-32.00$	$1.75-2.00$	${5.5}_{-3.3}^{+5.3}\times {10}^{-3}$	${5.1}_{-3.2}^{+5.4}\times {10}^{-3}$
$16.00-32.00$	$2.00-2.50$	${5.46}_{-0.89}^{+1.11}\times {10}^{-2}$	${5.4}_{-1.2}^{+1.0}\times {10}^{-2}$
$16.00-32.00$	$2.50-3.00$	${3.88}_{-0.68}^{+0.79}\times {10}^{-2}$	${3.80}_{-0.70}^{+1.06}\times {10}^{-2}$
$16.00-32.00$	$3.00-4.00$	${2.07}_{-0.39}^{+0.59}\times {10}^{-2}$	${2.14}_{-0.62}^{+0.60}\times {10}^{-2}$
$16.00-32.00$	$4.00-6.00$	${4.6}_{-2.8}^{+4.4}\times {10}^{-3}$	${3.4}_{-1.6}^{+2.8}\times {10}^{-3}$
$16.00-32.00$	$6.00-8.00$	${2.8}_{-1.1}^{+1.7}\times {10}^{-3}$	${2.6}_{-1.8}^{+1.7}\times {10}^{-3}$
$16.00-32.00$	$8.00-12.00$	${2.1}_{-1.1}^{+1.5}\times {10}^{-3}$	${1.9}_{-1.0}^{+1.4}\times {10}^{-3}$
$16.00-32.00$	$12.00-16.00$	${2.58}_{-1.23}^{+0.67}\times {10}^{-3}$	${2.6}_{-1.6}^{+2.1}\times {10}^{-3}$
$32.00-64.00$	$0.50-0.75$
$32.00-64.00$	$0.75-1.00$	${4.1}_{-2.3}^{+4.0}\times {10}^{-2}$	${3.7}_{-2.5}^{+4.9}\times {10}^{-2}$
$32.00-64.00$	$1.00-1.25$	${6.0}_{-2.5}^{+2.7}\times {10}^{-2}$	${5.1}_{-1.9}^{+2.6}\times {10}^{-2}$
$32.00-64.00$	$1.25-1.50$	${9.5}_{-6.3}^{+11.3}\times {10}^{-3}$	${8.4}_{-5.5}^{+9.6}\times {10}^{-3}$
$32.00-64.00$	$1.50-1.75$	${1.44}_{-0.70}^{+1.11}\times {10}^{-2}$	${1.33}_{-0.90}^{+0.85}\times {10}^{-2}$
$32.00-64.00$	$1.75-2.00$	${2.3}_{-1.2}^{+1.1}\times {10}^{-2}$	${2.1}_{-1.2}^{+1.1}\times {10}^{-2}$
$32.00-64.00$	$2.00-2.50$	${5.0}_{-1.2}^{+1.5}\times {10}^{-2}$	${4.8}_{-1.2}^{+1.5}\times {10}^{-2}$
$32.00-64.00$	$2.50-3.00$	${6.6}_{-1.3}^{+1.4}\times {10}^{-2}$	${6.2}_{-1.3}^{+1.5}\times {10}^{-2}$
$32.00-64.00$	$3.00-4.00$	${2.07}_{-0.76}^{+0.90}\times {10}^{-2}$	${2.09}_{-0.79}^{+1.01}\times {10}^{-2}$
$32.00-64.00$	$4.00-6.00$	${4.6}_{-2.6}^{+3.1}\times {10}^{-3}$	${4.0}_{-2.4}^{+2.8}\times {10}^{-3}$
$32.00-64.00$	$6.00-8.00$	${5.0}_{-2.5}^{+3.0}\times {10}^{-3}$	${5.0}_{-2.3}^{+2.8}\times {10}^{-3}$
$32.00-64.00$	$8.00-12.00$	${2.8}_{-1.7}^{+2.7}\times {10}^{-3}$	${2.7}_{-1.9}^{+2.4}\times {10}^{-3}$
$32.00-64.00$	$12.00-16.00$	${5.7}_{-1.9}^{+2.8}\times {10}^{-3}$	${5.2}_{-2.0}^{+3.1}\times {10}^{-3}$
$64.00-128.00$	$0.50-0.75$
$64.00-128.00$	$0.75-1.00$
$64.00-128.00$	$1.00-1.25$	${9.8}_{-4.9}^{+6.3}\times {10}^{-2}$	${7.6}_{-3.5}^{+5.6}\times {10}^{-2}$
$64.00-128.00$	$1.25-1.50$	${1.8}_{-1.3}^{+2.0}\times {10}^{-2}$	${1.24}_{-0.95}^{+1.75}\times {10}^{-2}$
$64.00-128.00$	$1.50-1.75$	${2.8}_{-1.8}^{+2.3}\times {10}^{-2}$	${2.5}_{-1.4}^{+2.1}\times {10}^{-2}$
$64.00-128.00$	$1.75-2.00$	${2.0}_{-1.2}^{+1.8}\times {10}^{-2}$	${1.8}_{-1.1}^{+1.4}\times {10}^{-2}$
$64.00-128.00$	$2.00-2.50$	${6.2}_{-1.9}^{+2.1}\times {10}^{-2}$	${5.5}_{-1.6}^{+1.9}\times {10}^{-2}$
$64.00-128.00$	$2.50-3.00$	${3.9}_{-1.7}^{+1.4}\times {10}^{-2}$	${3.7}_{-1.4}^{+2.1}\times {10}^{-2}$
$64.00-128.00$	$3.00-4.00$	${3.6}_{-1.1}^{+1.3}\times {10}^{-2}$	${3.35}_{-0.89}^{+1.21}\times {10}^{-2}$
$64.00-128.00$	$4.00-6.00$	${1.48}_{-0.66}^{+0.62}\times {10}^{-2}$	${1.40}_{-0.52}^{+0.68}\times {10}^{-2}$
$64.00-128.00$	$6.00-8.00$	${4.8}_{-2.9}^{+3.8}\times {10}^{-3}$	${4.8}_{-2.5}^{+4.7}\times {10}^{-3}$
$64.00-128.00$	$8.00-12.00$	${8.6}_{-4.3}^{+4.3}\times {10}^{-3}$	${8.6}_{-3.5}^{+4.3}\times {10}^{-3}$
$64.00-128.00$	$12.00-16.00$	${4.9}_{-3.3}^{+4.3}\times {10}^{-3}$	${4.2}_{-2.3}^{+4.3}\times {10}^{-3}$
$128.00-256.00$	$0.50-0.75$
$128.00-256.00$	$0.75-1.00$
$128.00-256.00$	$1.00-1.25$
$128.00-256.00$	$1.25-1.50$
$128.00-256.00$	$1.50-1.75$	${3.9}_{-2.4}^{+2.0}\times {10}^{-2}$	${3.2}_{-2.0}^{+2.4}\times {10}^{-2}$
$128.00-256.00$	$1.75-2.00$	${2.5}_{-2.0}^{+2.4}\times {10}^{-2}$	${1.9}_{-1.3}^{+2.1}\times {10}^{-2}$
$128.00-256.00$	$2.00-2.50$	${6.2}_{-2.9}^{+2.6}\times {10}^{-2}$	${5.5}_{-2.0}^{+2.8}\times {10}^{-2}$
$128.00-256.00$	$2.50-3.00$	${7.5}_{-2.1}^{+2.4}\times {10}^{-2}$	${6.5}_{-2.3}^{+2.1}\times {10}^{-2}$
$128.00-256.00$	$3.00-4.00$	${3.0}_{-1.7}^{+1.9}\times {10}^{-2}$	${2.8}_{-1.1}^{+1.6}\times {10}^{-2}$
$128.00-256.00$	$4.00-6.00$	${1.05}_{-0.54}^{+0.74}\times {10}^{-2}$	${9.6}_{-5.5}^{+8.1}\times {10}^{-3}$
$128.00-256.00$	$6.00-8.00$	${9.9}_{-3.9}^{+9.3}\times {10}^{-3}$	${9.7}_{-5.2}^{+5.5}\times {10}^{-3}$
$128.00-256.00$	$8.00-12.00$	${2.09}_{-0.80}^{+0.96}\times {10}^{-2}$	${1.90}_{-0.65}^{+0.62}\times {10}^{-2}$
$128.00-256.00$	$12.00-16.00$
$256.00-500.00$	$0.50-0.75$
$256.00-500.00$	$0.75-1.00$
$256.00-500.00$	$1.00-1.25$
$256.00-500.00$	$1.25-1.50$
$256.00-500.00$	$1.50-1.75$	${1.66}_{-0.96}^{+1.19}\times {10}^{-1}$	${1.28}_{-0.75}^{+0.89}\times {10}^{-1}$
$256.00-500.00$	$1.75-2.00$	${1.59}_{-0.83}^{+1.02}\times {10}^{-1}$	${1.12}_{-0.66}^{+0.82}\times {10}^{-1}$
$256.00-500.00$	$2.00-2.50$	${5.9}_{-3.6}^{+5.4}\times {10}^{-2}$	${4.1}_{-2.7}^{+3.8}\times {10}^{-2}$
$256.00-500.00$	$2.50-3.00$	${10.0}_{-3.8}^{+5.9}\times {10}^{-2}$	${7.6}_{-3.9}^{+3.3}\times {10}^{-2}$
$256.00-500.00$	$3.00-4.00$	${3.4}_{-2.0}^{+2.4}\times {10}^{-2}$	${2.7}_{-1.7}^{+2.4}\times {10}^{-2}$
$256.00-500.00$	$4.00-6.00$	${3.2}_{-1.5}^{+1.5}\times {10}^{-2}$	${2.6}_{-1.5}^{+2.1}\times {10}^{-2}$
$256.00-500.00$	$6.00-8.00$
$256.00-500.00$	$8.00-12.00$	${2.1}_{-1.1}^{+1.5}\times {10}^{-2}$	${1.62}_{-0.70}^{+1.07}\times {10}^{-2}$
$256.00-500.00$	$12.00-16.00$

Table 3. Table 3: HZ Occurrence Rate Estimates a a The true rate of planets in this period range could be lower by as much as a factor of ∼ 2 similar-to absent 2 \sim~{}2 , due to uncertainty in the reliability of planet candidates (see § 4.3 & Thompson et al. , 2018 ; Beaumont et al. , 2009 ) . Upper limits are believed robust to ∼ 10 % similar-to absent percent 10 \sim~{}10\% (see § 4.3 ).

	Prior
Radius	Dirichlet		Uniform
(R⊕)	$f$	$\Gamma$	$f$	$\Gamma$
$0.75-1.00$
$1.00-1.25$
$1.25-1.50$
$1.50-1.75$	${0.077}_{-0.040}^{+0.061}$	${0.67}_{-0.35}^{+0.53}$	${0.12}_{-0.06}^{+0.07}$	${1.06}_{-0.53}^{+0.63}$
$1.75-2.00$	${0.12}_{-0.05}^{+0.06}$	${1.2}_{-0.5}^{+0.6}$	${0.14}_{-0.07}^{+0.09}$	${1.39}_{-0.72}^{+0.90}$
$0.75-1.50$
$1.00-1.75$	${0.17}_{-0.07}^{+0.09}$	${0.40}_{-0.07}^{+0.10}$	${0.35}_{-0.14}^{+0.21}$	${0.83}_{-0.34}^{+0.50}$

Table 4. Table 4: ExoPaG Occurrence Rates

Period	Radius	Type F	Type G	Type K
(days)	($R_{\oplus }$)
$10.00-20.00$	$0.67-1.00$	${1.32}_{-0.65}^{+0.98}\times {10}^{-1}$	${3.9}_{-1.7}^{+2.1}\times {10}^{-2}$	${1.27}_{-0.49}^{+0.59}\times {10}^{-1}$
$10.00-20.00$	$1.00-1.50$	${1.23}_{-0.34}^{+0.37}\times {10}^{-1}$	${7.1}_{-1.6}^{+1.7}\times {10}^{-2}$	${7.2}_{-1.7}^{+1.9}\times {10}^{-2}$
$10.00-20.00$	$1.50-2.25$	${3.4}_{-1.4}^{+2.0}\times {10}^{-2}$	${5.14}_{-0.94}^{+1.16}\times {10}^{-2}$	${7.3}_{-1.3}^{+1.5}\times {10}^{-2}$
$10.00-20.00$	$2.25-3.38$	${2.4}_{-1.0}^{+1.1}\times {10}^{-2}$	${6.14}_{-0.86}^{+0.99}\times {10}^{-2}$	${8.1}_{-1.1}^{+1.2}\times {10}^{-2}$
$10.00-20.00$	$3.38-5.06$		${8.8}_{-4.0}^{+4.8}\times {10}^{-3}$	${6.5}_{-3.7}^{+5.7}\times {10}^{-3}$
$10.00-20.00$	$5.06-7.59$		${3.1}_{-1.6}^{+2.2}\times {10}^{-3}$	${5.7}_{-3.4}^{+2.3}\times {10}^{-3}$
$10.00-20.00$	$7.59-11.39$		${1.46}_{-0.86}^{+1.43}\times {10}^{-3}$	${4.0}_{-2.2}^{+3.9}\times {10}^{-3}$
$10.00-20.00$	$11.39-17.09$		${1.7}_{-1.0}^{+1.2}\times {10}^{-3}$
$20.00-40.00$	$0.67-1.00$	${2.3}_{-1.2}^{+1.9}\times {10}^{-1}$	${6.7}_{-3.8}^{+5.9}\times {10}^{-2}$	${1.51}_{-0.70}^{+0.91}\times {10}^{-1}$
$20.00-40.00$	$1.00-1.50$	${6.3}_{-2.7}^{+3.8}\times {10}^{-2}$	${4.7}_{-1.6}^{+2.1}\times {10}^{-2}$	${7.3}_{-2.7}^{+3.1}\times {10}^{-2}$
$20.00-40.00$	$1.50-2.25$	${6.3}_{-2.5}^{+2.7}\times {10}^{-2}$	${4.3}_{-1.2}^{+1.1}\times {10}^{-2}$	${8.3}_{-1.6}^{+1.9}\times {10}^{-2}$
$20.00-40.00$	$2.25-3.38$	${6.1}_{-2.1}^{+2.3}\times {10}^{-2}$	${9.7}_{-1.2}^{+1.3}\times {10}^{-2}$	${6.7}_{-1.2}^{+1.4}\times {10}^{-2}$
$20.00-40.00$	$3.38-5.06$	${1.15}_{-0.75}^{+1.59}\times {10}^{-2}$	${9.3}_{-4.4}^{+6.2}\times {10}^{-3}$	${1.30}_{-0.67}^{+0.91}\times {10}^{-2}$
$20.00-40.00$	$5.06-7.59$	${1.21}_{-0.63}^{+1.03}\times {10}^{-2}$	${2.6}_{-1.5}^{+2.1}\times {10}^{-3}$	${6.2}_{-3.8}^{+7.6}\times {10}^{-3}$
$20.00-40.00$	$7.59-11.39$			${5.8}_{-2.9}^{+4.2}\times {10}^{-3}$
$20.00-40.00$	$11.39-17.09$		${3.6}_{-1.9}^{+2.9}\times {10}^{-3}$
$40.00-80.00$	$0.67-1.00$
$40.00-80.00$	$1.00-1.50$		${8.2}_{-3.0}^{+4.0}\times {10}^{-2}$	${9.3}_{-4.3}^{+5.2}\times {10}^{-2}$
$40.00-80.00$	$1.50-2.25$	${4.9}_{-2.3}^{+3.8}\times {10}^{-2}$	${4.7}_{-1.6}^{+2.1}\times {10}^{-2}$	${9.8}_{-2.3}^{+3.1}\times {10}^{-2}$
$40.00-80.00$	$2.25-3.38$	${3.8}_{-1.6}^{+2.9}\times {10}^{-2}$	${1.34}_{-0.19}^{+0.21}\times {10}^{-1}$	${7.1}_{-1.9}^{+2.8}\times {10}^{-2}$
$40.00-80.00$	$3.38-5.06$		${5.9}_{-4.0}^{+6.6}\times {10}^{-3}$	${2.6}_{-1.2}^{+1.5}\times {10}^{-2}$
$40.00-80.00$	$5.06-7.59$		${1.03}_{-0.42}^{+0.91}\times {10}^{-2}$	${6.9}_{-4.2}^{+6.4}\times {10}^{-3}$
$40.00-80.00$	$7.59-11.39$		${5.0}_{-3.3}^{+4.5}\times {10}^{-3}$	${7.5}_{-3.7}^{+7.5}\times {10}^{-3}$
$40.00-80.00$	$11.39-17.09$		${5.7}_{-2.9}^{+4.4}\times {10}^{-3}$	${8.3}_{-5.1}^{+7.0}\times {10}^{-3}$
$80.00-160.00$	$0.67-1.00$
$80.00-160.00$	$1.00-1.50$			${1.67}_{-0.89}^{+1.03}\times {10}^{-1}$
$80.00-160.00$	$1.50-2.25$	${9.7}_{-6.2}^{+5.4}\times {10}^{-2}$	${7.7}_{-2.8}^{+3.7}\times {10}^{-2}$	${1.01}_{-0.29}^{+0.46}\times {10}^{-1}$
$80.00-160.00$	$2.25-3.38$	${6.8}_{-4.0}^{+5.1}\times {10}^{-2}$	${1.14}_{-0.19}^{+0.29}\times {10}^{-1}$	${1.22}_{-0.28}^{+0.34}\times {10}^{-1}$
$80.00-160.00$	$3.38-5.06$		${2.8}_{-1.1}^{+1.4}\times {10}^{-2}$	${1.5}_{-1.0}^{+1.4}\times {10}^{-2}$
$80.00-160.00$	$5.06-7.59$	${3.1}_{-1.8}^{+3.9}\times {10}^{-2}$	${4.9}_{-3.1}^{+5.6}\times {10}^{-3}$	${1.97}_{-0.94}^{+1.32}\times {10}^{-2}$
$80.00-160.00$	$7.59-11.39$		${2.13}_{-0.87}^{+0.89}\times {10}^{-2}$
$80.00-160.00$	$11.39-17.09$		${1.14}_{-0.57}^{+1.06}\times {10}^{-2}$
$160.00-320.00$	$0.67-1.00$
$160.00-320.00$	$1.00-1.50$
$160.00-320.00$	$1.50-2.25$		${1.26}_{-0.49}^{+0.81}\times {10}^{-1}$	${1.27}_{-0.61}^{+0.69}\times {10}^{-1}$
$160.00-320.00$	$2.25-3.38$	${1.57}_{-0.64}^{+1.17}\times {10}^{-1}$	${9.6}_{-2.9}^{+4.7}\times {10}^{-2}$	${1.17}_{-0.39}^{+0.57}\times {10}^{-1}$
$160.00-320.00$	$3.38-5.06$		${2.0}_{-1.2}^{+1.7}\times {10}^{-2}$	${4.5}_{-2.5}^{+2.8}\times {10}^{-2}$
$160.00-320.00$	$5.06-7.59$		${2.7}_{-1.4}^{+1.3}\times {10}^{-2}$
$160.00-320.00$	$7.59-11.39$	${4.4}_{-2.6}^{+3.9}\times {10}^{-2}$	${2.10}_{-0.96}^{+1.58}\times {10}^{-2}$	${3.0}_{-1.7}^{+4.7}\times {10}^{-2}$
$160.00-320.00$	$11.39-17.09$
$320.00-640.00$	$0.67-1.00$
$320.00-640.00$	$1.00-1.50$
$320.00-640.00$	$1.50-2.25$	${6.7}_{-3.7}^{+3.5}\times {10}^{-1}$	${6.3}_{-2.4}^{+2.6}\times {10}^{-1}$	${3.9}_{-2.0}^{+2.8}\times {10}^{-1}$
$320.00-640.00$	$2.25-3.38$	${3.8}_{-2.2}^{+3.1}\times {10}^{-1}$	${1.22}_{-0.63}^{+1.01}\times {10}^{-1}$	${2.64}_{-0.96}^{+1.45}\times {10}^{-1}$
$320.00-640.00$	$3.38-5.06$		${3.5}_{-2.5}^{+3.3}\times {10}^{-2}$	${1.32}_{-0.73}^{+1.12}\times {10}^{-1}$
$320.00-640.00$	$5.06-7.59$		${3.0}_{-1.5}^{+3.3}\times {10}^{-2}$
$320.00-640.00$	$7.59-11.39$
$320.00-640.00$	$11.39-17.09$

Equations

$$p_{\det ,DR25}(SNR)=c\times \gamma (\alpha ,\beta \times SNR)/\Gamma (\alpha ),$$

$$p_{\mathrm{vet,Mulders}}(R_{p},P)=cR_{p}^{a_{R}}\left\{\begin{array}[]{ll}\left(P/P_{\mathrm{break}}\right)^{a_{P}},~{}\mathrm{if~{}}P<P_{\mathrm{break}}\\ \left(P/P_{\mathrm{break}}\right)^{b_{P}},~{}\mathrm{otherwise}\end{array}\right.,$$

$$p_{det\&vet}(SNR,N_{tr})=c_{N_{tr}}\times \gamma ({\alpha }_{N_{tr}},{\beta }_{N_{tr}}\times SNR)/\Gamma ({\alpha }_{N_{tr}}),$$

$$SNR=d\ast \sqrt{N_{tr}\ast D}/\mathrm{CDPP}$$

$$\rho (s_{\mathrm{obs},k},s_k^{\ast })=\sum _k\frac{|s_{\mathrm{obs},k}-s_k^{\ast }|}{\sqrt{s_{\mathrm{obs},k}+s_k^{\ast }}},$$

Full text

methods: data analysis — methods: statistical — catalogs — planetary systems — stars: statistics

††software: ABC (Ford et al., 2018), CORBITS (Brakensiek & Ragozzine, 2016), ExoplanetsSysSim (Ford et al., 2019), ExoJulia SysSimData (Ford, 2019), Matplotlib (Hunter, 2007)

1 Introduction

In 2009, NASA’s Kepler mission began its mission to characterize the abundance and diversity of exoplanets. Four years of precise photometry of 189,516 target stars (as part of its exoplanet survey) allowed the Kepler project to identify 4,034 strong planet candidates via transit, over half of which have been confirmed by another method or validated based on detailed analyses of light curves and statistical considerations (Thompson et al., 2018). The Kepler mission was enormously successful, with an impressive array of discoveries, both anticipated and surprising. However, one of the primary goals of the mission, characterizing the frequency of Earth-size planets in or near the habitable zone of solar-like stars has remained elusive. A combination of factors, from intrinsic stellar variability to failure of Kepler’s reaction wheels, resulted in the Kepler data being sensitive to such planets around a small fraction of the stars surveyed and the detection efficiency of such planets depending sensitively on their host star properties. As a result, considerable care is necessary when inferring the rate of small planets near the habitable zone from Kepler data.

Several previous studies have estimated the rate of planets as a function of planet size and orbital period based on Kepler data. Early studies had less data, but just as importantly, the accuracy of early studies were constrained by limited knowledge of the Kepler pipeline’s detection efficiency and the host star properties. Further complicating matters, Hsu et al. (2018) showed that one of the most common algorithms for estimating the planet occurrence rate was systematically biased for planets near Kepler’s detection threshold. While this bias is modest for much of the parameter space explored by Kepler, it can be of order unity for small planets near the habitable zone. Only a few recent studies analyzed the planet occurrence rates using a hierarchical Bayesian model, so as to avoid this bias (Foreman-Mackey et al., 2014; Burke et al., 2015; Hsu et al., 2018). Foreman-Mackey et al. (2014) inferred planet occurrence rates using a hierarchical Bayesian model with non-parametric model and a prior that favored the planet occurrence rate changing smoothly. This study analyzed a list of planet candidates from a non-standard pipeline that was restricted to finding a single planet around each star. Burke et al. (2015) inferred planet occurrence rates using a parametric model (i.e., power-law in planet size and orbital period) and a hierarchical Bayesian model. Hsu et al. (2018) inferred planet occurrence rates using a non-parametric model, a hierarchical Bayesian model and Approximate Bayesian Computing (ABC), and performed extensive tests on simulated catalogs to validate the algorithm. However, none of these studies made use of Kepler’s final planet catalog (known as DR25; Thompson et al., 2018), additional DR25 data products characterizing the detection pipeline’s efficiency (Christiansen, 2017; Burke & Catanzarite, 2017a, b, c; Coughlin, 2017), or recent improvements in the knowledge of host star properties from the ESA’s Gaia mission (Gaia Collaboration et al., 2018). Recent studies indicate that the updates in stellar properties thanks to Gaia can have a significant effect on the inferred planet occurrence rate (Shabram et al., 2019).

In this manuscript, we report occurrence rates inferred using a hierarchical Bayesian model and Approximate Bayesian Computing (ABC), and incorporating several improvements relative to previous studies, including the DR25 catalog of planet candidates, updates to the host star properties from Gaia, and an improved model for the Kepler efficiency for detecting and vetting of planet candidates. We review the statistical methodology in §2.1 and describe improvements to our model since Hsu et al. (2018) in §2.2. In §3.1, we describe the selection of target stars, and present the resulting planet occurrence rates in §3.2. We conclude with a summary of key findings, implications for the rate of Earth-size planets in or near the habitable zone of FGK stars, and future of occurrence rate studies in §4.

2 Methodology

2.1 Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) is a tool for performing Bayesian inference when the likelihood is not available. For example, ABC can be applied to perform inference using a hierarchical Bayesian model to infer the properties of a population based on observations from a survey. Real-world surveys can have many complexities (e.g., measurement uncertainties, target selection, instrumental effects, data reduction pipeline, detection probability that depends on unobserved properties) that make it impractical to write down the correct likelihood. Instead, ABC relies on: 1) the ability to generate samples from a forward-model for the intrinsic population and survey, and 2) a physically-motivated distance function that quantifies the distance between the actual survey results and the results of a simulated survey. For a sufficiently large number of simulated catalogs and a well-chosen distance function, the ABC posterior estimate converges around the true posterior. Thus, ABC is ideally suited to inferring planet occurrence rates based on transit surveys such as Kepler.

Here we provide a qualitative description of naive ABC and the sequential importance sampling method applied in this study. The occurrence rate of planets within a given range of sizes and orbital periods is characterized by model and a set of hyperparameters. We apply a model that is a piece-wise constant rate per logarithmic bin in orbital period and planet size. The hyperparameters are simply the occurrence rates for each bin in the 2-d grid. For this study, we chose a period-radius grid with similar period and radius bins: $P=$ {0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 500} days & $R_{p}=$ {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 6, 8, 12, 16} R*⊕*. We repeatedly draw planets from this model, simulate observations by the Kepler mission, compute a set of summary statistics (in this case the number of detected planets in each bin) and compare the number of planets (within a certain range of planet sizes and/or orbital periods) that would have been detected in the simulated dataset to the actual number of such planets detected by the Kepler mission. A naive application of ABC would involve drawing hyperparameters from their prior distribution and keeping parameters that resulted in matching the observed number of planets in each bin exactly. While straightforward, this approach is computationally unfeasible for our application.

The implementation of ABC used in this study is Population Monte Carlo, wherein multiple generations of simulated data are created. In each generation, a population of model parameters are used to generate and evaluate multiple simulated catalogs. Simulated catalogs with a distance less than some threshold (“$\epsilon$”) are kept and the values of their hyperparameters inform the values of hyperparameters proposed for the following generation. Each set of hyperparameter values is assigned a weight, and sequential importance sampling is used to ensure that the weights properly account for the prior probability distribution at each generation. This sequential importance sampling continues until one of several stopping criteria are met which indicate sufficient convergence has been achieved. For a full description of ABC and the particulars of the ABC-PMC algorithm we use, see Hsu et al. (2018).

While we generally apply the ABC procedure detailed in Hsu et al. (2018), there were a few modifications to account for this study. In this study, each generation of ABC-PMC consists of 100 to 500 “particles” (increased from 50) where each particle is a simulated catalog produced from a draw of occurrence rates from importance sampling. We found an increased number of particles to be helpful when applying the algorithm to infer several occurrence rates simultaneously, particularly when some of the rates were weakly constrained by the data (see §2.4). Instead of always using the same number of target stars and value of $\tau$ (which controls the width of the importance sampler’s proposal distribution), we split the ABC simulation into three steps. For the first 15 generations of each ABC simulation, we create simulated catalogs with only $N_{*}=1000$ targets and set $\tau=1.5$. This accelerates the ABC-PMC algorithm as it “zooms in” on the relevant region of parameter space. For future generations, we set $\tau=2$ and generate larger simulated catalogs of $N_{*}=10,000$ targets, so as to provide a good importance sampling density for future generations. The ABC-PMC algorithm continues for more generations until satisfying the stopping criteria described in Hsu et al. (2018). In the third step, if the simulation was using less than 200 particles, we run one more generation of ABC using 200 particles drawn using the distance goal ($\epsilon$) of the last generation to provide a larger posterior sample for plotting and computing percentiles. To prevent to selection of an outlier as the reported result, we perform five ABC-PMC simulations for each set of parameters and report the median value for each quantile.

2.2 Model Improvements

The majority of the physical model in the Exoplanet System Simulator (known as “SysSim”) used in this study has remained the same as in our previous application of ABC to estimate the occurrence rates of planet candidates around FGK stars based on the Q1-Q16 Kepler catalogs in Hsu et al. (2018) (see Section 3). SysSim (Ford et al., 2019) is available publicly on Github. The version used in this study which implements the older v0.6 release of the programming language Julia (Bezanson et al., 2017) is available at https://github.com/dch216/ExoplanetsSysSim.jl/tree/hsu_etal_2019-v1.0 while a version updated to be compatible with the current release of Julia is available at https://github.com/ExoJulia/ExoplanetsSysSim.jl/releases/tag/v1.0.0, with both repositories under the MIT Expat License. Below we summarize several improvements made for this study.

2.2.1 Stellar Properties

The second Gaia data release (DR2) (Gaia Collaboration et al., 2018) provides significantly improved stellar parameters for the overwhelming majority of the Kepler target stars. Improvements in stellar radii determinations result in more accurate planet radii (and more accurate assigning of planet candidates to their appropriate radii bins). For this study, we use the stellar radius from Gaia DR2 (Andrae et al., 2018). For comparing observations to simulations, we use the median transit depths computed from the MCMC posterior chains provided with the Kepler DR25 catalog (Hoffman & Rowe, 2017) (see §2.2.2: MCMC Transit Depths). Note that this differs from the values listed in the DR25 catalog which are based on a maximum likelihood fit. For assigning each planet to a radius bin, we also use planet-star radius ratios which we derived from the MCMC posterior chains provided along with the DR 25, accounting for limb-darkening, plus the transit depth and impact parameter of each state in the Markov chain. An initial cross-match was performed using https://github.com/megbedell/gaia-kepler.fun. We filter this target list further to remove poor planet search targets, as described in §3.1.

Note that we do not report results which use the stellar radii derived by Berger et al. (2018), as that catalog makes use of additional observations (e.g., spectroscopy, asteroseismology) that are not available for all stars. Comparing the stellar radii values for the same targets reveals that the Berger catalog generally reports larger values than Gaia DR2 for stars with radii $R_{*}<0.75R_{\odot}$ and vice versa at larger radii. Additionally, the Berger catalog reports a much narrower spread of radii uncertainties than the Gaia catalog, although the Gaia catalog has smaller uncertainties on average. Following-up on a suggestion from the anonymous referee, we repeated our ABC simulations making use of stellar properties from the Berger catalog instead of Gaia DR2. Typically, we found very similar results, but there were some modest differences in regions of the period-radius grid which have the most planet candidates and thus best measurement precision. While incorporating additional observational constraints (as done by the Berger catalog) likely improves the precision for individual stars, using such an inhomogeneous set of stellar parameters risks biasing our results since high-resolution spectroscopic observations are available for nearly all stars with known planet candidates, but only a small fraction of stars without any planet candidates.

Additionally, Gaia’s precise parallax measurements allow for more accurate identification of main-sequence stars based on their position in the color-luminosity diagram. We use this information (in place of estimates based on $\log g$) to define a clean sample of FGK main-sequence target stars. Finally, we make use of Gaia’s astrometric information to identify targets likely to be multiple star systems with stars of similar masses/luminosities. The details of our target selection are described in §3.1. This results in a total of $79,935$ Kepler targets for which the Kepler DR25 pipeline and robovetter identified $2,525$ planet candidates with $P=0.5-500$d and $R_{p}=0.5-16R_{\oplus}$.

2.2.2 Incorporating Kepler DR 25 data products

Planet Catalog: This study makes use of the Kepler DR 25 planet catalog and its associated data products (Thompson et al., 2018). This represents a substantial improvement over the Q1-16 catalog used in Hsu et al. (2018). The DR25 catalog makes use of the light curves from the full Kepler prime mission and incorporates several improvements to the data reduction, planet search, and vetting algorithms. Even more importantly, DR25 was the first Kepler catalog to be generated using a uniform data reduction, planet search and vetting process. For planet occurrence studies, uniformity of processing is particularly important, as it enables statistical modeling of the detection process.

MCMC Transit Depths: We replace the catalog maximum-likelihood estimator value of transit depth with the median value of the transit depth posterior, which was determined by taking the DR25 MCMC posterior chains (Hoffman & Rowe, 2017) for each planet candidate and applying the quadratic limb-darkening law from Mandel & Agol (2002). The median depth value is then used to determine the planet-star radius ratio and therefore planet radius in a consistent manner with the simulated planets.

Planet Detection Efficiency: When computing occurrence rates, it is important to use an accurate detection efficiency model for the probability that a planet with known properties is included in the result catalog of planet candidates. Hsu et al. (2018) used a detection efficiency model from Christiansen et al. (2015) that was calibrated to transit injection tests performed using a previous pipeline. In this study, we update the planet detection efficiency model based on transit injection tests based on the DR 25 pipeline Christiansen (2017). Initially, we adopted the detection probability for a transiting planet as a function of the expected effective signal-to-noise ratio (SNR),

[TABLE]

where $\gamma$ is the incomplete gamma function, $\Gamma$ is the gamma function, $\alpha=30.87$, $\beta=0.271$, and $c=0.94$, based on §5 of Christiansen (2017). This expression models the probability that a transiting planet is detected by the pipeline as a threshold crossing event, but does not attempt to model the probability that a threshold crossing event is promoted to a Kepler Object of Interest or a planet candidate. The Kepler pipeline identifies threshold crossing events for further investigation based on a “multiple event statistic” or MES exceeding a threshold of 7.1. The MES is defined to be what the pipeline returns and can not be computed exactly without applying the Kepler pipeline to a lightcurve. As described later in this section, our definition of expected SNR represents the best available approximation of the expected MES in the absence of applying the Kepler pipeline to a light curve. Christiansen (2017) suggest that occurrence rate studies fit custom detection efficiency models that are optimized for the specific stellar sample and range of orbital periods being studied. While most previously published occurrence rate studies have ignored this advice, we fit a custom detection efficiency model to the Christiansen (2017) simulations for this study’s baseline model (described below). We explicitly give Eqn. 1, because we also report results using Eqn. 1 so as to facilitate comparisons to previous studies and our baseline model.

Most previous studies have simply assumed that all true planets that are identified as threshold crossing events survive the vetting process. Mulders et al. (2018) presented a first attempt to estimate the impact of the vetting process. In this study, we consider three models for the combined planet detection process (i.e., being identified as a threshold crossing event by the pipeline and labeled as a planet candidate by the robovetter). First, we adopt $p_{\mathrm{det,DR25}}$ from Eqn. 1 and assume perfect vetting (i.e., $p_{\mathrm{vet}}=1$). Second, we adopt the product of $p_{\mathrm{det,DR25}}$ and the Mulders et al. (2018) model for vetting efficiency,

[TABLE]

where $R_{p}$ is planet radius (in R⊕), $P$ is the orbital period (in days), $c=0.93,a_{R}=-0.03,P_{\mathrm{break}}=205,a_{P}=0.00$ and $b_{P}=-0.24$111These parameter values differ from the values reported in Mulders et al. (2018) which were based on fitting the same functional form to planet candidates with a vetting score of $\geq~{}0.9$ and using the DR25 stellar parameters. This study uses updated parameters provided by G. Mulders that are updated to use Gaia DR2 stellar parameters and make use of all Kepler DR 25 KOIs that are labeled as a planet candidate by the robovetter, regardless of the vetting score.. Note that this model implies a vetting efficiency that does not explicitly depend on the effective transit SNR, but rather has only an indirect dependence via the planet size and orbital period. Since the power of the vetting process depends primarily on the ability of the robovetter to features in the shape of the light curve, it is likely that the vetting efficiency indeed depends more directly on the effective measurement noise for the target star. Another concern with this model is that it assumes the probability of a planet being detected by the pipeline and a planet passing vetting are uncorrelated. Since both the pipeline detection efficiency and vetting efficiency depend on the transit light curve (and thus key properties such as the transit signal to noise and the number of transits observed), assuming that these are independent seems unwise. These (as well as results described in §2.5) motivated us to create a new model of the planet detection and vetting process.

For the third model, we derive a combined model for the probability of a transiting planet being detected and labeled as a planet candidate by the robovetter using the pixel-level transit injection tests (Christiansen, 2017) and the associated robovetter results (Coughlin, 2017).

[TABLE]

where $N_{tr}$ is the number of “valid” transits observed by Kepler and values for $\alpha$, $\beta$, and $c$ are given in Table 1. (The Kepler photometric reduction pipeline assigned weights to each flux measurement, so as to deweight observations that may be spurious for a variety of reasons (e.g., measurement near a data gap or anomaly). Following Christiansen (2017), a transit is labeled as “valid” if the central flux measurement during the transit has a weight greater than 0.5.)

We follow the general procedure recommended in Christiansen (2017) for fitting a probability as a function of the effective SNR. However, we only consider an injected transit to have been detected if it also is labeled as a planet candidate by the robovetter. The choice of $N_{tr}$ as the second model parameter was motivated by the observation that the vetting efficiency decreases for planets with longer orbital periods and thus smaller $N_{tr}$. Theoretically, we expect that the dispersion in the measured SNR will depend directly on the number of transits that are contributing to the candidate and only indirectly on orbital period. For example, if two planets have a similar orbital period, but fewer transits of one planet were observed (e.g., due to falling on the CCD module that died early in the Kepler mission), then the distribution of effective SNR would be broader for the planet with fewer observed transits. We adopt this as our baseline model and compare the results using these three models for the combined planet detection process in §4.1.

$1-\sigma$* Depth Function:* Each of the above planet detection models depends on computing the expected effective SNR. In Hsu et al. (2018), we used:

[TABLE]

where the numerator represents the “signal”, $N_{tr}$ is the number of transits, $D$ is the transit duration, and $d$ is the fractional transit depth. For “noise”, we took the mission average of the 4.5 hour duration combined differential photometric precision (CDPP) for each target star and assumed that the fractional noise scaled with one over the square root of the number of measurements.

In this study, we make use of a new data product provided by DR25, the “1-$\sigma$ depth function” (OSDF), which describes the true transit depth that would be expected to result in a MES of unity for a given target star, orbital period, and transit duration, after averaging over the epoch of transit (Burke & Catanzarite, 2017a). Thus, the expected MES and effective SNR become simply the transit depth divided by the OSDF for the target star interpolated to the orbital period and appropriate transit duration. The CDPP is effectively an averaged version of the OSDF that assumes a specific transit duration and averages over orbital period and phase. Thus, we expect the OSDF to be a more accurate measure of the effective noise as it accounts for how the photometric noise deviates from white noise. For planets with a modest number of transits, the SNR can be affected by whether they occur at times with more or less noise than typical. This is particularly significant for planets with orbital periods greater than 90 days, since there is a smaller number of transits and the number of transits observed with one CCD module may differ significantly from the number of transits observed with another CCD module.

The original Kepler DR25 OSDFs are defined for every star, all 15 pulse durations that were searched by the pipeline, and for each of a a large ($\sim 10^{4}$) number of orbital periods (Burke & Catanzarite, 2017a). In order to reduce the memory requirements, we downsample these through linear interpolation to a common period grid for all stars of 1000 logarithmically-spaced periods from 0.49 to 700 days, which captures the vast majority of the information in these OSDFs.222The OSDFs for each star are available from the Kepler mission at https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_completeness_reliability.html. Our downsampled versions are in a more convenient format and available on GitHub at https://github.com/ExoJulia/SysSimData. Though the amount of memory required for SysSim to use OSDFs is much higher than before, it’s run-time performance did not increase significantly.

When evaluating the noise of a particular planet, we typically use bi-linear interpolation between the period and duration. However, OSDFs are not defined for all periods and durations because the Kepler Transit Search pipeline did not search every combination of period and duration. When the simulated planet’s duration was shorter than the shortest duration searched (and available in OSDF), we assume that the SNR was diluted by the ratio of the searched duration to the actual duration since the search pipeline effectively adds noise to such short duration transits.

When the simulation planet’s duration was longer than the longest duration searched at that period, we adopt the noise of the longest duration searched. The SNR is reduced due to the signal only accumulating for the duration searched rather than the true duration.

The “timeout” corrections for which stars and durations were searched in DR25 have not been incorporated to our model. This affects only a small fraction of our target stars, as discussed in §4.3.

Window Function: For the Kepler pipeline to detect a planet candidate, there must be good data for at least three transits. This becomes non-trivial for long period planets and/or stars not observed for the full mission duration. Following (Burke et al., 2015), Hsu et al. (2018) adopted a binomial model for the probability that at least three transits were observed, based the expected number of transits if there were no data gaps and the target-specific duty cycle.

In this study, we make use of the DR25 Window Function data products which tabulate the probability of detecting at least three transits as a function of orbital period (Burke & Catanzarite, 2017a). As with OSDFs, the provided Kepler DR25 Window Functions are a function of star, duration, and period. These were downsampled to 1000 linearly-spaced periods from 0.5 to 700 days using linear interpolation.333The original versions are also at https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_completeness_reliability.html and our downsampled versions are available on GitHub at https://github.com/ExoJulia/SysSimData. For a simulated planet, SysSim uses bilinear interpolation on these window functions to determine the value of the window function for each simulated planet. As in Hsu et al. (2018), the SNR-detection probability is multiplied by this value to return the final probability of detecting a planet.

Since window functions are identical for targets that were observed in the same sequence of quarters, we identify the 100 most common window functions for targets that were observed for at least 4 quarters as part of the Exoplanet target list (determined by requiring "EX" to be in the Investigation ID as recorded on MAST (archive.stsci.edu/kepler). These 100 window functions provide an exact match for practically all of our targets. Any targets for which the exact window function was not available were assigned a randomly drawn window function.

Minor improvements affecting detection efficiency: We also incorporate several minor model improvements when calculating the effective signal to noise. For example, Hsu et al. (2018) ignored limb-darkening. In this study, we account for limb-darkening when computing the transit depth, using the limb-darkening parameters from the Kepler DR 25 stellar catalog, the same values as used for the Kepler pipeline run that produced the DR25 planet candidate catalog. (Note: These limb-darkening parameters are based off of Kepler DR25 stellar parameters, so as to be consistent with DR25 data products such as the OSDF.) Additionally, we account for dilution of transit depth due to contamination as tabulated in the Kepler Input Catalog. Finally, we have improved the calculation of the transit duration, by using Eqn. 15 of Kipping (2010), so as to avoid making the small angle approximation.

2.3 Sky-averaged detection probabilities using CORBITS

In order to improve the wall-time efficiency of our ABC calculation, we have implemented a probabilistic simulated catalog approach. Previously, when simulating observations to determine the detection probability of each simulated planet, we assigned a geometric transit probability of either zero or unity depending on whether the planet would appear to transit for a single observer orientation. In this study, we assign each potentially detectable planet a weight proportional to its transit detection probability marginalized over all possible observer orientations. The geometric transit probability for each planet is multiplied by a transit detection efficiency that is marginalized over all impact parameters (that result in a transit) to give the total detection probability. For a single planet, the sky-averaged geometric transit probability for each planet is $R_{\star}/[a(1-e^{2})]$, where $R_{\star}$ is the stellar radius, $a$ is the semi-major axis and $e$ is the orbital eccentricity.

For stars with multiple planets, SysSim calculates the sky-averaged geometric transit probability for each planet, as well as each pair of planets and each number of transiting planets in the system. The CORBITS package (Brakensiek & Ragozzine, 2016) makes this computationally practical by using semi-analytic methods to determine the transit probability of each combination of planets associated with a given planetary system. 444CORBITS requires that orbits not “cross”, i.e., if $a_{1}<a_{2}$, then it requires that $a_{1}\times(1+e_{1})<a_{2}\times(1-e_{2})$. Physically, such configurations are unlikely due to long-term orbital stability. Therefore, when multiple planets are assigned to a single star, we reject configurations that result in crossing orbits and redraw until a configuration without crossing orbits is created (up to a maximum of 20 times).

Using the sky-averaged detection probability instead of the detection probability for a single viewing geometry reduces the number of simulated planetary systems needed to accurately infer planet occurrence rates. As part of our model and algorithm verification, we determined that the number of planetary systems drawn could be reduced by a factor of $\simeq~{}8$ relative to the number of target stars without significantly affecting the accuracy of the inferred occurrence rates. This reduction offsets the overhead from CORBITS, which increases the calculation time from single observer by a factor of $\simeq~{}6$. Therefore, the overall reduction is about one-third off the previous calculation time. While this provides only a modest performance improvement for this study, we anticipate that the sky-averaged observing geometry mode will be particularly valuable for future studies that investigate the occurrence rate of systems with multiple transiting planets.

In Hsu et al. (2018), the distance function used by ABC was simply the absolute value of the difference in the ratios of detected planets to stars surveyed for the simulated and observed planet catalogs. Since this study averages over viewing geometries, simply summing the detection probabilities would result in the expected number of planet detections which is not directly comparable with the observed number of planets of the actual Kepler data. The expected number of planet detections has a smaller variance than the actual number of planet detections. Therefore, after calculating the detection probabilities for each planet in a simulated catalog, we label each planet as either detected or not-detected based on a Bernoulli draw with success probability equal to each planet’s total detection probability.

2.4 Inferring multiple occurrence rates simultaneously

In Hsu et al. (2018), we inferred the occurrence rate for each “bin” in period and radius independently of other bins. Orbital periods are measured so precisely that there is effectively no ambiguity about which period bin a planet should be assigned to. In principle, three effects could cause a planet to be assigned to a different bin than it would be if its properties were known perfectly. First, measurement errors cause the observed planet-star radius ratio to differ from it’s true value. Fortunately, this happens for only a small fraction of planets. Indeed, Hsu et al. (2018) found that this effect did not result in significant correlation between neighboring planet radii bins. Second, uncertainties in the stellar properties cause the measured planet radius to differ from the true planet radius, even if the planet-star radius ratio were known precisely. Nearly all previously published studies, including Hsu et al. (2018), have neglected this effect. Shabram et al. (2019) showed that uncertainties in host star radii could have a significant effect on planet occurrence rates. In this study, we infer multiple occurrence rates simultaneously, so as to account for this effect. A third source of uncertainty (unrecognized contamination) is discussed in 4.3, but is not modeled in this study.

Once Gaia provided accurate parallaxes for most Kepler target stars (Gaia Collaboration et al., 2018), accurate stellar radii and uncertainties were derived for those targets (Andrae et al., 2018). This makes it feasible to accurately model the effects of uncertainty in stellar radii on planet occurrence rates. Since planet radii are proportional to the host star radius, uncertainty in host star radius directly translate into uncertainty in the planet radius. Even with Gaia DR2 parallaxes, the uncertainty in stellar radius can result in a planet being assigned to the wrong radius bin. For each planet we simulate, the true transit depth is based on the true planet radius and the stellar radius from Gaia (Andrae et al., 2018). We draw an assumed stellar radius from an equal mixture of two half-normal distributions with median equal to the Gaia best-fit radius and use Gaia’s upper and lower uncertainties for the standard deviations of the half-normals (Andrae et al., 2018). If the uncertainty in stellar radius (in either direction) is less than $6\%$ the value of the radius, then we increase that uncertainty to $6\%$ the stellar radius (based on the distribution of stellar radius uncertainties in Berger et al. (2018)). Then, we draw the observed transit depth centered on the true transit depth with a width based on its SNR and the diagonal noise model of Price & Rogers (2014) that accounts for finite integration time. Next, we compute the “observed” planet radius from the observed plant-star radius ratio and the assumed stellar radius. In most cases, the observed planet radius results in the planet being assigned to the same bin as it would be if its radius were known precisely. However, sometimes, the observed planet radius results in the planet being assigned to a neighboring bin containing slightly larger or smaller planets. If one were to include the uncertainty in stellar properties, but only modeled planets drawn from a single bin, then the inferred planet occurrence rate would be biased to be larger than the true rate, since simulated planets can “leak” beyond of the boundaries of the radius bin. The magnitude of the effect depends on the size of the bin in radius relative to the size of uncertainties in stellar properties and on the differences between occurrence rates in neighboring bins. In tests, we found that this effect could lead to occurrence rates being inflated by $\sim~{}20\%$ of the measured rate (for bins focusing on small planets with a width of 0.25 R⊕), if each bin were analyzed individually, rather than in a group.

2.4.1 Model Parameterization

Fortunately, our ABC approach naturally accounts for this effect, as long as we simultaneously model the occurrence rate of planets in each size bin. When performing inference with occurrence rates for multiple radius bins, we performed calculations using two sets of priors, each with their own model parameterization.

First, we assign each bin in period and radius for the entire grid its own occurrence rate, with a uniform prior. We adopt a uniform prior for $f_{i,j}$ over $[0,f_{\mathrm{\max,i,j}})$. The upper limit was set to $f_{\mathrm{\max,i,j}}=C\times\log(P_{\max,j}/P_{\min,j})/\log(2)\times\log(R_{\max,i}/R_{\min,i})/\log(2)$, with $C=2$. This upper limit is small enough that proposals with more than 3 planets per factor of 2 in period are extremely rare, consistent with expectations based on long-term orbital stability. This has the advantage that the prior for any individual bin is easy to understand. However, it has the disadvantage that the total number of planets within a range of orbital periods could be larger than plausible if one considers long-term orbital stability.

When deriving our best estimate of planet occurrence rates in the $\eta_{\oplus}$ regime, we use a different parameterization. For each bin in orbital period (indexed by $j$), we infer: 1) $f_{\mathrm{tot},j}$, the total occurrence rate for all planets within the range of orbital periods and the full range of planet sizes being considered, and 2) a vector $f_{\mathrm{rel},i,j}$: the fraction of planets in the $j$th orbital period bin that have sizes falling within the range of the $i$th radius bin. Thus, the occurrence rate for planets in the $i$th radius bin and $j$th period bin is given by $f_{i,j}=f_{\mathrm{tot,j}}f_{\mathrm{rel},i,j}$. We adopt a uniform prior for $f_{\mathrm{tot,j}}$ over $[0,f_{\mathrm{\max,tot,j}})$. The upper limit $f_{\mathrm{\max,tot,j}}=3\times\log(P_{\max,j}/P_{\min,j})/\log(2)$ is motivated by long-term orbital stability, as it is very rare for more than 3 planets to have orbital periods within a factor of 2 of each other (though a few cases are known to exist for very small planets). For each $f_{\mathrm{rel},i,j}$ (a vector with $j$ fixed), we adopt a Dirichlet prior with concentration parameters proportional to $\log(R_{\max,i}/R_{\min,i})$. The smallest concentration parameter is set to unity. The Dirichlet prior ensures that $\sum_{i}f_{\mathrm{rel},i,j}=1$ to within numerical precision. This prior has the advantage that the rate of planets summed over all sizes within one period bin is constrained. Since the Dirichlet prior over the fractions in each radius bin makes it more challenging to visualize the prior, we report results and show posterior distributions for both sets of priors in §2.4.

For most periods and sizes, the choice of prior did not have a significant impact on the posterior for occurrence rates. The differences were more noticeable for small, long-period planets. Therefore, we will report results for both choices of prior in the $\eta_{\oplus}$ regime.

2.4.2 Distance Function

Now that we allow for uncertainty in stellar properties, we revisit the choice of distance function. We performed tests on simulated data to verify our algorithms and to compare the performance of multiple distance functions. Simply summing the absolute value of the difference in the ratios of detected planets to stars surveyed for the simulated and observed planet catalogs often resulted in an ABC posterior with highly disparate widths for different bins. Since bins with a greater number of observed planets have a greater contribution to the total distance, ABC would result in precise occurrence rates for well-populated bins, but much less precise estimates of occurrence rates for bins with low rates of planet candidates. For many scientific purposes, one would prefer a similar fractional error in occurrence rates for all bins, rather than a similar absolute error in occurrence rates for all bins. Therefore, using the Canberra distance (Lance & Williams, 1967) as a foundation we choose a new distance function that weights the absolute value of the difference in the ratio of detected planets to targets for each bin by the square root of the sum of the ratio of detected planets to targets for the simulated and observed planet catalogs for that bin. Our distance is given by

[TABLE]

where, for each $k$th period-radius bin, $s_{\mathrm{obs},k}$ is the ratio of number of planet candidates detected by Kepler to the number of target stars searched and $s^{*}_{k}$ is ratio of the number of planets detected to target stars in the simulated catalog. We found that this distance function allows ABC to converge more rapidly when inferring occurrence rates for multiple bins simultaneously. We tested that this algorithm accurately simulates both the number of planets detected and the dispersion in the number of planets detected for the simulated catalogs. We also tested that this distance function accurately estimates the true planet occurrence rate and its uncertainty (by comparing to simulations that do not make use of sky-averaged detection probabilities; see §2.5 and Figure 1).

2.4.3 Importance Sampler Proposal Distribution

Initially, we attempted to apply the sequential importance sampler in the ABC-PMC algorithm using a multivariate Gaussian for the proposal distribution, as described in Beaumont et al. (2009) and Hsu et al. (2018). However, this proved computationally prohibitive. For example, a proposal of a negative value for any one bin would result in rejecting the entire proposal; when considering several bins that have small occurrence rates, such proposals are common. In order to make the ABC-PMC algorithm practical, we experimented with a variety of replacement proposal distributions. Eventually, we settled on using a Beta distribution for the proposal for both parameterizations. When using independent priors for each $f_{i,j}$, the ABC-PMC algorithm uses a Beta proposal distribution for each $f_{i,j}$. When using a Dirichlet prior for the relative occurrence rates within a given period bin, the ABC-PMC algorithm using a Beta proposal for $f_{\mathrm{tot},j}$ and each $f_{\mathrm{rel},i,j}/f_{\mathrm{\max,tot,j}}$. After each proposal, the initially proposed values for the relative rates were rescaled according to $f_{\mathrm{rel},i,j}\leftarrow f_{\mathrm{rel},i,j}/\sum_{i}f_{\mathrm{rel},i,j}$, so that the rescaled rates sum to unity. We validated the algorithm on simulated data sets and found that this dramatically improved the computational efficiency.

2.5 Verification and Validation

Before applying our improved model and ABC algorithm to actual Kepler catalogs, we performed extensive tests to validate and and verify the model, the algorithms, and our implementations.

Figure 1 shows a comparison of results using either the sky-averaging mode (described in §2.3) or the single observing geometry mode for one particular bin (orbital periods of 8-16 days and planet sizes of 1.25-1.5R⊕) where the ”true“ catalog is set to $f=3\%$. The blue (orange) points show the median planet occurrence rate inferred for this bin using a single observing geometry (sky-averaged viewing geometry) and the vertical error bars indicate the 68.3% credible region. The horizontal axis indicates the number of occurrence rates inferred simultaneously. For any single simulated catalog, the number of observed planets is likely greater (or less) than the expected number of observed planets, leading to the posterior median rate being larger (or smaller) than the true rate ($f=3\%$). We performed several simulations to verify that the median posterior rate from multiple runs was distributed about the true value. Figure 1 shows the results of multiple analyses for a single catalog demonstrating that we recover a posterior median close to the true value.

First, we will focus on results for three bins. Note that SysSim yields very similar results, regardless of whether using a single observing geometry (blue) or when using averaging over all viewing geometries (orange). As expected, including stellar radii uncertainties (green) results in an increased width of posterior distribution relative to an analysis which assumes stellar radii are known perfectly (orange).

Next, we compare results as a function of the number of occurrence rates inferred simultaneously. Each method (single observing geometry or sky-averaged viewing geometries, excluding or including uncertainties in stellar properties) gives statistically indistinguishable results for the median occurrence rate, as long as the number of model parameters being inferred simultaneously does not exceed seven. When analyzing three or five radius bins simultaneously, the posterior widths are very similar. In either case, the posterior width is observed to increase modestly when the model accounts for uncertainties in stellar radii. When analyzing nine radius bins simultaneously, the posterior width has increased noticeably. This behavior is typical of the sequential importance sampler used by the ABC-PMC algorithm. When performing inference on nine parameters simultaneously, it becomes less likely that the proposed values for all nine parameters will result in a precise match for the number of observed planets in all nine bins. As a result, the computational time required for ABC-PMC to achieve a given distance increases significantly. In practice, the final distance achieved by the ABC-PMC algorithm begins to increase substantially as the number of dimensions increases beyond seven, resulting in an increased width of the ABC-posterior. Thus, we must balance the desire to minimize bias with the desire to minimize unnecessary increase in posterior width due to computational limitations. Minimizing bias leads us to perform inference on at least three bins simultaneously. Avoiding an unnecessary increase in posterior width limits us to performing inference on $\leq~{}7$ bins simultaneously.

Next, we performed several tests to further fine tune the choice of how many occurrence rates to infer simultaneously. When inferring occurrence rates for multiple radius bins simultaneously, we find that only the top and bottom bins display noticeable bias. This is expected since the model effectively assumes that there are no detectable planets with radius greater than the upper limit of the top bin or less than the lower limit of the bottom bin. Since the bias is proportional to the difference between the true and assumed rate of such planets, it is maximized when the assumed rate is zero. When inferring occurrence rates for just three bins simultaneously, one could worry that the estimated rates for both edge bins will be overestimated and this could have an indirect effect on the estimated rate for the central bin. In practice, the bias for the central bin is a second-order effect and significantly smaller than the other uncertainties. A related effect occurs when two bins (neighboring in radii) have a large difference in occurrence rate. If the “leaking” of planets out-of and in-to the radius bin in question is not symmetric, then the impact of stellar uncertainties is amplified compared to other sets of bins where the occurrence rate is similar above and below the bin in question. This demonstrates that it is important to model the uncertainties in stellar radii and to simultaneously infer the planet occurrence rates when estimating planet occurrence rates for a wide range of planet radii and orbital periods.

Based on the above tests, we settled on reporting results using simulations that simultaneously infer planet occurrence rates for 5 (or 7) radius bins, reporting the results for the interior 3 (or 5) radius bins. (One exception is that we also report the occurrence rate for our largest bin, 12-16 R⊕, as explained below.) In order to span the full range of planet radii and to avoid significantly overestimating the ABC-posterior width, we perform multiple simulations, each using 5 to 7 radius bins. The specific bin boundaries chosen for the subsets were: $R_{p}=$ {0.25, 0.5, 0.75, 1, 1.25, 1.5}, {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5}, {1, 1.25, 1.5, 1.75, 2, 2.5}, {1, 2, 2.5, 3, 4, 6}, {3, 4, 6, 8, 12, 16} R*⊕* for each period range.

Formally, one should be cautious in using results from an “edge bin”. Therefore, for each subset we throw out the two edge bins and take estimates for only interior bins from each subset. In the case of bins which are interior to multiple subsets, we combine the final populations produced by the separate runs on each subset to produce the occurrence rate estimates. The one edge bin which we do not throw out is the 12-16 R⊕ bin. Even though the 12-16 R⊕ bin is an edge, the rate of planets larger than 16 R⊕ is so small that there is no practical impact on the measured occurrence rates for 12-16 R⊕ planets. While we include occurrence rates for $0.25-0.5$ R⊕ planets in our model calculations, we do so primarily so as to ensure that our estimates for 0.5-0.75 R⊕ planets are accurate. Thus, we do not include the inferred rates for $0.25-0.5$ R⊕ planets in Figure 2 out of an abundance of caution.

3 Application to DR25

3.1 Catalog Selection

The second Gaia data release (DR2) (Gaia Collaboration et al., 2018) provides significantly improved stellar parameters for the overwhelming majority of the Kepler target stars. Improvements in stellar radii determinations result in more accurate planet radii and more accurate assigning of planet candidates to their appropriate radii bins. Additionally, Gaia’s precise parallax measurements allow for more accurate identification of main-sequence stars based on their position in the color-luminosity diagram. We use this information (in place of estimates of $\log g$) to define a clean sample of FGK main-sequence FGK target stars. Finally, we make use of Gaia’s astrometric information to identify targets likely to be multiple star systems with stars of similar masses/luminosities.

We construct a catalog of Kepler target stars starting from the Kepler DR25 stellar properties catalog and augmented with data from Gaia DR2 (Gaia Collaboration et al., 2018). The selection criteria used to select Kepler targets likely to be main-sequence FGK stars are as follows (applied in order):

Require that the Kepler magnitude and the Gaia G magnitude are consistent. We compute the median and standard deviation of the difference in the Kepler magnitude and Gaia G magnitude based on the initial cross-matched catalog based on position. If a target’s Kepler magnitude minus Gaia G magnitude differs from the median by more than 1.5 times the standard deviation, then it is excluded from our cleaned catalog. This filter makes sure that a Kepler target is not erroneously matched to a Gaia target corresponding to a background star. 2. 2.

Require Gaia GOF_AL $\leq 20$ and Gaia astrometric excess noise $\leq 5$ cut. These criteria indicate a poor astrometric fit. This is often caused by an unresolved astrometric binary star (Evans, 2018). 3. 3.

Require Gaia Priam processing flags that indicate the parallax value is strictly positive and both colors are close to the standard locus for main-sequence stars. This rejects Kepler targets which are unlikely FGK main sequence stars or are so distant that Gaia does not have a good parallax measurement and the stellar radius will be highly uncertain. 4. 4.

Require Gaia parallax error is less than $10\%$ the parallax value. This excludes some faint/distant Kepler targets for which accurate stellar radii would not be available. Of all our target selection criteria based on data quality this results in the largest decrease in the number of target stars and planet candidates. In principle, future observations (e.g., improved Gaia data releases) could increase the number of stars with accurate stellar properties, perhaps allowing for more precise planet occurrence rate measurements. However, for FGK main-sequence stars, this criterion only excludes distant and hence faint targets for which Kepler does not have significant sensitivity to small planets in the habitable zone. Therefore, we expect that even future Gaia data releases and other follow-up observations to improve stellar characterization will not lead to significant improvements in the constraints on the occurrence rate of small planets in or near the habitable zone of main-sequence FGK stars. 5. 5.

Require that Kepler DR 25 provide a valid Kepler stellar mass, data span, duty cycle, and limb-darkening coefficients, and that Gaia’s Apsis module provide a valid stellar radius based on Gaia data (Andrae et al., 2018). 6. 6.

Require that the Kepler target was observed for $>4$ quarters and must have been on the Exoplanet target list for at least one quarter. These criteria exclude targets that were observed by Kepler for purposes other than the exoplanet search, as these are unlikely to be FGK main sequence stars. A small number of stars that were part of the exoplanet search will also be excluded if they were only observed for $\leq~{}4$ quarters. Kepler does not have significant sensitivity to small planets in the habitable zone of such FGK stars. 7. 7.

Require the target to have a color $0.5\leq B_{p}-R_{p}\leq 1.7$ from Gaia photometry. This color cut results in selecting FGK stars and is more precise than using the temperature from the Kepler Input Catalog and more uniform than using temperatures from the DR25 stellar catalog. 8. 8.

Require the target star to have a luminosity, $L\leq 1.75L_{MS}(B_{p}-R_{p})$, where $L_{MS}(x)=10^{2.62-3.74x+0.962x^{2}}$ represents the luminosity of a main-sequence star for a given $B_{p}-R_{p}$ color. To compute $L_{MS}(B_{p}-R_{p})$ we fit a quadratic function to the observed values of $\log_{10}L$ as a function of $B_{p}-R_{p}$ for our Kepler targets. We reject those targets whose observed $L$ deviates from $L_{MS}$ by more than 75% and refit the model to the remaining targets. We iterate this process six times before converging on our final model for the main-sequence luminosity $L_{MS}$. This rejects targets that are significantly above the main sequence, either due to stellar evolution (primarily for F stars) or due to a multiple star system where stars other than the primary contribute >75% of the total flux. While this final criteria likely removes some viable targets for planet hunting (e.g., F stars starting to evolve off the main sequence), it is a small fraction of the total stars searched and Kepler does not have significant sensitivity to small planets near their habitable zone.

The total number of FGK Kepler target stars remaining after these cuts is $79,935$. While the number of targets is significantly less than the total number observed by Kepler, a clean sample of target stars is preferable for performing planet occurrence rate studies. In principle, a less restrictive set of cuts might provide a larger stellar sample and provide more precise estimates for occurrence rates of large planets. However, such a strategy is unlikely to be useful for measuring the occurrence rate of small planets near the habitable zone, since the vast majority of stars excluded were not FGK main-sequence stars for which Kepler had significant sensitivity to small planets in or near the habitable zone. Our selection criterion intentionally exclude M stars, which will be the subject of a future study.

The Kepler DR25 pipeline and robovetter identified $2,524$ planet candidates associated with these targets with $P=0.5-500$d and $R_{p}=0.5-16R_{\oplus}$. One noticeable change from our previous study is the relatively few planet candidates in the $\eta_{\oplus}$ regime. The primary cause for this change is the usage of updated stellar radii from Gaia, which are often larger than the previously determined Kepler stellar radii. This trend tends to boost the inferred planet radii. As a result, the majority of long period, small radii planet candidates shifted to larger radii bins within the period-radius grid once we incorporated Gaia stellar radii. For example, the inferred radius KOI 7016.01(Jenkins et al., 2015; Mullally et al., 2018, also known as Kepler-452 b) reported in DR25 is $1.06^{+0.2}_{-0.1}$ R⊕, but incorporating the updated stellar radii from Gaia DR results in the estimated radius increasing to 1.51 R⊕ (with uncertainties increasing proportionally). While the best-estimate for the planetary radius no longer falls within the $1-1.25$ R⊕ bin, the uncertainty in the radius of such planets results in it to contributing to the estimated occurrence rate.

In order to explore which $\eta_{\oplus}$ regime planet candidates were removed due to our stellar cuts (as opposed to updated stellar radii), we created a separate target list where we relax the two stellar cuts which most significantly reduced the number of stars in our sample: the FGK luminosity cut and the cuts designed to remove suspected binaries. This catalog has a total number of $139,232$ target stars, with $3,170$ planet candidates within the limits of the period-radius grid. Most significantly, we recover two long-period, small radius planet candidates (associated with targets KIC 5097856 and 5098334) in the $P=256-500$d, $R_{p}=1.25-1.5R_{\oplus}$ bin. After investigating the properties of these two planet candidates, we determined that the candidates were not in our final catalog because their associated target stars had poor astrometric GOF which suggests that their host star is likely part of an unresolved binary. If true, then the unmodeled flux from the binary companion would be diluting the transit depth, causing the true planet radius to be larger than currently estimated. We conclude that our process for identifying a clean sample of target stars worked as intended.

3.2 Baseline Planet Occurrence Rates

The values in each bin of Figure 2 state the median of the ABC-posterior for the occurrence rates over the full period-radius grid. Results are tabulated in Table 2. The uncertainties shown are the differences between the median and either the 15.87th or 84.13th percentile (whichever has the larger absolute difference) of the ABC-posterior for each occurrence rate. These rates make use of the combined detection and vetting efficiency (Eqn. 3). For period-radius bins with zero or one observed planet candidates, we report the 84.13th percentile as an upper limit instead of a median rate and uncertainty. For these bins, the posterior distribution is asymmetrical and summarizing the posterior with a point estimate runs a greater risk of misinterpretation. Figure 5 depicts the ABC posterior for the occurrence rate for several long period, small radius bins to demonstrate this asymmetry.

The reported upper limits are conservative in the sense that a full calculation (i.e., inferring rate for all radius bins simultaneously within a period range, rather than just 5 at a time) would likely cause a modest reduction in the upper limit. This is because we set a maximum total rate of 3 per factor of 2 in period (as a stability criterion), so a full calculation which adds well-constrained bins at larger radii would reduce the remaining total rate that can be allotted to the bins with only upper limits. Given the occurrence rates inferred, the effect is negligible for most of the parameter space considered. While a few effects (e.g., unmodeled contaminating flux, detection efficiency decreasing for systems with multiple transiting planets, pipeline timeouts) could cause our model to overestimate the detection efficiency for individual planetary systems, each of these affect only a small fraction of stars (see §4.3). Therefore, we expect that any revisions to our upper limits due to these effects would be limited to a few percent. However, for bins with so little constraining data the estimated posterior may be significantly affected by the choice of prior (see §4.2).

In Figure 3, we show the occurrence rate as a function of planet size, marginalized over two separate period ranges, 0.5-8d and 16-128d. The planet radius valley Fulton et al. (e.g., 2017); Van Eylen et al. (e.g., 2018) is apparent for periods of 16-128d, but not for periods less than 8 days. Our use of a cleaned target star sample is designed to reduce effects due to contamination from binary stars, similar to Teske et al. (2018). The depth and width of the valley are expected to depend on orbital period on theoretical grounds (Lopez & Fortney, 2014; Owen & Wu, 2017; Gupta & Schlichting, 2019). A narrow and deep valley whose location depends on period would appear to be broader and shallower when one marginalizes over a broad range of periods, as in this figure. Therefore, the depth of the valley shown in Figure 3 should be interpreted as a minimum depth when comparing to predictions of model predictions over a narrower range of periods. Comparing the rate of 1.25-1.5 R⊕ planets as a function of orbital periods, we find no evidence that the occurrence rate is decreasing beyond 64 days, as predicted by some models.

\startlongtable

4 Discussion

Using Kepler DR25, Gaia DR2, SysSim and ABC model, we provide accurate estimates of the exoplanet occurrence rate around FGK main sequence stars as a function of planet size and orbital period. Our results span a wide range of orbital periods and planet radii, including the $\eta_{\oplus}$ regime, without the need to extrapolate. Our Bayesian approach provides direct and meaningful constraints on the $\eta_{\oplus}$ regime, despite the limited number of such planet candidates (see §4.2). We find planet occurrence rates somewhat lower than some previous studies (e.g., Figure 17 of Burke et al., 2015). Just as importantly, we derive statistically valid uncertainties in planet occurrence rates, while accounting for many effects that have often been neglected in previous studies. The slightly lower rate and larger uncertainty may have significant implications for scientists planning future planet surveys or characterization missions.

We performed extensive tests of our model and examined the relative importance of various upgrades to the model. One of the most significant factors resulting in differences with the Q1-Q16 rates presented in Hsu et al. (2018) is the usage of Gaia DR2 stellar parameters. In addition to improving the accuracy of planet radii, the Gaia data enabled the selection of a clean sample of FGK main sequence target stars for occurrence rate studies. The number of small planets at long orbital periods was significantly reduced after focusing on a clean FGK target star sample and incorporating improved stellar radii. This is one factor that contributes to our inferred occurrence rates being lower than some previous studies.

We find that it is important to account for the vetting efficiency when characterizing the occurrence rate of small planets or planets with orbital periods greater than a month (see §4.1). We provide an improved model for the combined detection efficiency (i.e., including both being identified by the Kepler planet search pipeline and labeled as a planet candidate by the robovetter) in Eqn. 3. We conclude with a discussion of the areas for further improvement in future occurrence rate studies (§4.3) and implications for mission planning (§4.3.1).

4.1 Implications of Vetting Efficiency

Previous studies of planet occurrence rates have assumed that all true planets detected by Kepler’s Transiting Planet Search module are labeled as planet candidates in the Kepler planet catalog. With Kepler’s DR25 data release, vetting of threshold crossing events was performed by an automated process, the “robovetter”. This opened the door to explicitly modeling the vetting efficiency, i.e. the probability with which the Kepler robovetter identifies a true exoplanet as a planet candidate.

In order to explore the importance of vetting efficiency, we calculated the occurrence rate using three different models:

• •

Using the detection efficiency model in Eqn. 1 and setting the vetting efficiency to unity. This is analogous to previous exoplanet occurrence rate studies.

• •

Using the detection efficiency model in Eqn. 1 and adopting the vetting efficiency in Eqn. 2, a function of orbital period and planet radius based on the model of Mulders et al. (2018).

• •

Using a single model for the combined detection and vetting efficiency in Eqn. 3.

The occurrence rates produced by these three different models are compared for two bins in different detection regimes in Fig. 4. When comparing results using a vetting efficiency of unity to results using the Mulders vetting efficiency model, we find that the Mulders model can significantly affect the inferred planet occurrence rate, particularly for small, long-period planets (bottom panel). As the occurrence rate of such planets is of particular interest, this motivated us to develop an improved model for the vetting efficiency. We make use of the pixel-level transit injections that were processed by both the Kepler planet search pipeline and the robovetter to empirically fit a model for the probability of a true transiting planet being both detected by the pipeline and labeled as a planet candidate.

Our results demonstrate the importance of both accounting for the vetting efficiency and of using an accurate vetting efficiency model. When comparing the results with our model for the combined detection and vetting efficiency, we find significant changes in the inferred planet occurrence rates relative to previous calculations that neglected vetting efficiency, only for bins with $P>32$d or $R_{p}<2R_{\oplus}$. For the remainder of the bins, the difference in the median of the ABC-posteriors is typically considerably less than the width of the ABC-posterior.

We compared the results using our combined detection and vetting efficiency to results using the Mulders model in two regimes: short orbital periods with well-characterized occurrence rates with $P=8-16$d and $R_{p}=1.75-3R_{\oplus}$ (using 3 bins), and long orbital periods with weaker constraints on the occurrence rates ($P=237-500$d and $R_{p}=0.75-2R_{\oplus}$ (using 5 bins). We find that in both of these cases, the inferred occurrence rates are consistent within measurement uncertainties. The Mulders vetting model yields a mean rate only slightly higher (no more than $\sim~{}12\%$) in the first case. In the long-period regime, the Mulders vetting model finds a lower rate. While the fractional difference in the mean occurrence rates is up to $\sim~{}33\%$, this is still significantly smaller than the uncertainty due to the Kepler sample size. We conclude that our inferences for planet occurrence rates are robust to the choice of vetting model.

4.2 Occurrence Rates of Planets near the Habitable Zone

Given the widespread interest in the occurrence rate of small planets in or near the habitable zone ($\eta_{\oplus}$), we performed additional calculations for the occurrence rate of planets with radii $\{0.75-1$, $1-1.25$, $1.25-1.5$, $1.5-1.75$, $1.75-2\}$ R⊕ and orbital periods 237-500 days using different choices of window functions and vetting efficiencies.

4.2.1 Importance of Individual Model Improvements

Our results using an independent uniform prior for each bin’s occurrence rate can be seen in Fig. 4 (lower panel) where we present occurrence rate estimates assuming:

1. no loss of planets due to the vetting process and the simple binomial window function (comparable to the methodology of Hsu et al. (2018)),
2. no loss of planets due to the vetting process and DR25 target-by-target window functions,
3. combined detection and vetting efficiency and DR25 target-by-target window functions, and
4. Mulders et al. (2018) vetting efficiency and DR25 target-by-target window functions.

The inclusion of a target-by-target set of window functions can result in diminishing the inferred occurrence rates in the $\eta_{\oplus}$ regime due to better accounting for detection probability when there are few transits. While this effect is small compared to measurement uncertainties for the smallest planets (i.e., $<1.25R_{\oplus}$), it can be significant for larger planets.

Accounting for the vetting efficiency can significantly affect the inferred occurrence rate compared to previous estimates that have assumed all planets pass the vetting process. The occurrence rates based on the combined detection and vetting efficiency model generally increases relative to assuming perfect vetting. This demonstrates the importance of accounting for the vetting process when computing planet occurrence rates.

For orbital periods beyond 128 days, Kepler observations place interesting upper limits on the occurrence rate of Earth-size planets, even if there are no detections in a given period-radius bin. Table 2 and Figure 2 provide upper limits based on adopting independent uniform priors for each bin’s occurrence rate. While appealing in its simplicity, adopting independent uniform priors for each bin is equivalent to assuming a total rate of planets in that period range that is peaked at large planet occurrence rates. This can be understood by considering the sum of $n$ multiple independent random variables, each uniformly distributed between zero and $A$. The probability distribution for the sum is an Irwin-Hall distribution which has mean of $nA/2$ and a standard deviation of $\sqrt{nA/12}$. Assuming independent uniform distributions for the rate of planets in each of several radius bins is equivalent to assuming a prior distribution for the total rate that is more concentrated than assuming a uniform prior for the total rate. The location of the prior mode is at half of the sum of the upper limits of the individual priors. Thus, the shape of the posterior at low rates is sensitive the choice of prior, when there are no (or few) detections. Even in the case of no detections, the posterior still drops off rapidly towards high rates, since the lack of detections strongly excludes high rates. These effects can be seen clearly by comparing the dotted lines in panels f and l of Figure 5.

While the choice of prior did not have a significant effect for most of the bins in our period-radius grid, it does affect the rates derived in the $\eta_{\oplus}$ regime. Given the relatively weak observational constraints, the use of a peaked prior on the total rate of small planets in the habitable zone seems unwise.

Therefore, for the results presented in §4.2.2, we performed additional simulations specifically focusing on this regime using a Dirichlet prior over radius bins with boundaries subset $R_{p}=$ {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5}. This allows us to compute the posterior distribution for the rate of small, long-period planets, while assuming a uniform prior on the rate of planets of all sizes between 0.5 and 2.5 R⊕. We do not include larger planets, so as to improve the computational efficiency of our ABC algorithm. The rate of planets in each bin larger than 2.5 R⊕ is well constrained and the sum is less than $\simeq~{}40\%$. In principle, one could attempt to use the rate of larger planets to motivate a narrower prior on the rate of planets between 0.5 and 2.5 R⊕. Based on the posteriors computed, we confirm that including larger planets would have a negligible effect on the ABC posterior. Therefore, we do not reduce the upper limit for the total rate in $P=237-500d$ and $R_{p}=0.5-2.5$ R⊕ planets, as any corrections would be modest and the reliability of the larger planet candidates may not be well known.

4.2.2 $\eta_{\oplus}$ Rate Estimates

Our ABC posterior distributions for $\eta_{\oplus}$ (or $\Gamma_{\oplus}$) summarize the current state of our knowledge about the rate of small planets in the habitable zone of FGK stars (see panels a, c, e, g, i and k of Figure 5). We summarize our results by reporting upper limits (84.13th percentile) and median (50th percentile) occurrence rate of small planets over the period range $P=237-500$d using our combined detection and vetting efficiency and two choices of prior. Note that the our recommended rates in this regime are based on the Dirichlet prior for relative rates, and so these rates differ from the rates reported in Table 2 and Figure 2 based on independent uniform priors for each bin. Comparing the results with the two priors can be useful as a test to establish how sensitive these occurrence rates are to the choice of prior. Summing over the three interior bins, the {5, 15.87, 50, 84.13, 95}th percentiles for the summed occurrence rate assuming the Dirichlet prior over the range $P=237-500$d, $R_{p}=$ $0.75-1.5$ R*⊕* are $f_{R,P}=$ {$0.06$, $0.10$, $0.16$, $0.27$, $0.36$}. The corresponding percentiles for the differential rates are $\Gamma_{\oplus}\equiv(d^{2}f)/[d(\ln P)~{}d(\ln R_{p})]=$ { $0.12$, $0.19$, $0.32$, $0.51$, $0.69$} for the same range. We consider the above rates to represent the current best estimates.

For the sake of examining the sensitivity of the results to the choice of prior, we performed a similar analysis, adopting an independent uniform prior for each radius bin. The results can be seen in panels b, d, f, h, j and l of Figure 5. Using this prior, the {5, 15.87, 50, 84.13, 95}th percentiles for the summed occurrence rate over the range $P=237-500$d, $R_{p}=$ $0.75-1.5$ R*⊕* are $f_{R,P}=$ {$0.27$, $0.40$, $0.67$, $1.1$, $1.3$} while the differential rates are $\Gamma_{\oplus}=$ { $0.51$, $0.78$, $1.3$, $2.0$, $2.5$}. While we expect use use of independent uniform priors to result in an overestimate of summed occurrence rates, we include these values to quantify the sensitivity of these results to the choice of prior. While these two priors are quite different (e.g., compare dotted lines for panels k & l), the ABC posteriors are qualitatively similar for large rates (where data provide strong constraints), but differ significantly for small rates (where the posterior shape is primarily based on information coming from the prior).

4.2.3 Comparison to previous studies

The median of our posterior for the occurrence rate of small planets in or near the habitable zone of sun-like stars (see Fig. 5, panel k) is similar to estimates based on extrapolation from recent studies and well within the range of possibilities suggested by sensitivity analyses (e.g., Burke et al., 2015; Mulders et al., 2018; Bryson et al., 2019). Fig. 6 shows the ABC posterior distribution for $\Gamma_{\oplus}$ from our study, along with the results of several previous studies. In order to make a meaningful comparison, each studies results have been converted into a differential rate (using units of $R_{\oplus}$ and days, and natural log rather than base 10 or base 2). Other studies are summarized with a median and separate upper and lower uncertainties.

The most discrepant estimates are readily understood as due to using a much earlier and very different catalog of planets (i.e., Youdin, 2011; Petigura et al., 2013; Foreman-Mackey et al., 2014; Dong & Zhu, 2013). Only Mulders et al. (2018) and Bryson et al. (2019) used the Kepler DR25 planet catalog. Results for $\Gamma_{\oplus}$ from this study are intermediate, i.e., between the estimates from these studies. Equally important, we suggest that the uncertainty estimate from our study is more realistic than those of either Mulders et al. (2018); Bryson et al. (2019). While the Burke et al. (2015) analysis used a Q1-16 planet catalog, they still found a similar rate to Mulders et al. (2018).

Both the Burke et al. (2015) and Mulders et al. (2018) studies arrived at higher estimates for $\Gamma_{\oplus}$ than in this analysis, perhaps due to their use of a parametric model for the planet occurrence rate (i.e., power law or broken power-law in planet radius or planet-star radius ratio). Since the rate of planets in the $\eta_{\oplus}$ regime is only weakly constrained by data, assuming such a parametric distribution, essentially results in an extrapolation in the $\eta_{\oplus}$ regime and a large proportion of the inferred rate is based on planets that are only rarely detectable in the Kepler sample.

In contrast, a preliminary draft of Bryson et al. (2019) reports lower estimates for $\Gamma_{\oplus}$ than this study, regardless of whether they attempt to account for potential false alarms (i.e., “No Reliability” versus “With Reliability”). There is significant overlap between our posterior $\Gamma_{\oplus}$ and that of Bryson et al. (2019). We attribute the offset in the mean of our posteriors for $\Gamma_{\oplus}$ and their “No Reliability” estimate of $\Gamma_{\oplus}$ as likely due to some combination of their use of a parametric model, the different selection criteria for target stars, and/or difference in the stellar properties (based on an unpublished catalog, similar to Berger et al. (2018)). There is a more significant different in our posterior for $\Gamma_{\oplus}$ and the “With Reliability” estimate of Bryson et al. (2019). As discussed in §4.3, we caution that our posteriors for occurrence rates at $P>250$d and for $\Gamma_{\oplus}$ might overestimate the true rate by as much as a factor of $\sim~{}2$ due to potential false alarms at long orbital periods, but our upper limits should be robust, even at these long orbital periods. Using the stellar sample and stellar properties from Bryson et al. (2019) results in $2-4$ planet candidates in the $\eta_{\oplus}$ regime. Their analysis finds that these have low “reliability” and deweights these planets accordingly. We note that none of those planet candidates directly contribute to our estimate of $\Gamma_{\oplus}$, as they were either removed or assigned a radius larger than $1.5R_{\oplus}$, given our set of stellar properties and target selection criteria. According to Bryson et al. (2019), larger planets tend to have greater reliability, and thus the impact of false alarms on our planet occurrence rates are likely mitigated.

We have demonstrated that there are substantial uncertainties in occurrence rates in and near the $\eta_{\oplus}$ regime, making extrapolation even more dangerous. Therefore, we caution against overinterpreting apparent differences between this study and Burke et al. (2015), Mulders et al. (2018), or Bryson et al. (2019), as the differences in the inferred rates among these three studies are smaller than the uncertainties in $\Gamma_{\oplus}$ (which are primarily set by the effective number of stars searched).

Our upper limits are significantly more accurate and robust than previous upper limits, due to a combination of effects including the use of Kepler’s final DR25 planet candidate catalog, updated and homogeneous stellar properties from Gaia, and use of DR25 data products such as the target-specific window functions. While Burke et al. (2015), Mulders et al. (2018), and Bryson et al. (2019) report more precise values for $\Gamma_{\oplus}$, these are based on assuming a parametric model for the planet occurrence rate. They infer model parameters based on the occurrence rate of planets at much shorter orbital periods, and assume that these planets provide information about the occurrence rate in the habitable zone. For example, Burke et al. (2015) measure a differential occurrence rate between $P=50-300$d and assume that it remains unchanged at larger orbital periods. In contrast, this study reports a direct measurement of $\Gamma_{\oplus}$ between $P=237-500d$. While our upper limits should be robust, the lower limits are sensitive to the choice of prior and the unknown reliability of planet candidates in this regime (see §4.3 & Thompson et al., 2018; Bryson et al., 2019).

4.2.4 Extrapolations to Longer Periods

While the aforementioned studies have offered estimates of $\eta_{\oplus}$ for FGK stars (i.e., Petigura et al., 2013; Foreman-Mackey et al., 2014; Burke et al., 2015; Mulders et al., 2018), we caution that such estimates involved extrapolating a model that is constrained based on planets with larger sizes and/or shorter orbital periods. On one hand, we do not find compelling evidence for rapidly changing occurrence rates as a function of orbital period. On the other hand, we demonstrate that the uncertainties in planet occurrence rates in and near the eta-Earth regime are substantial, in some cases of the same order of magnitude as the median value. Extrapolation based on these rates to larger orbital periods risks amplifying errors. While it’s possible that the occurrence rate of planets as a function of planet size and period varies sufficiently slowly to render such extrapolations accurate, there are both observational and theoretical reasons to be cautious about such extrapolations. Theoretically, the origin and evolution of planets with significant gas envelopes could quite different than that of rocky planets. Indeed, studies have suggested a local minimum in the occurrence rates as a function of planet size, perhaps due to photoevaporation or core-powered mass-loss (Owen & Wu, 2017; Gupta & Schlichting, 2019). Such a feature is not captured by parametric models which assume planet occurrence rates to be proportional to a power-law (or broken power-law) in planet radius. Additionally, there may be qualitative changes in the occurrence rate and properties of planets as orbital periods approach that of the ice line. From an observational perspective, radial velocity surveys find a substantial rise in the frequency of more massive planets starting at orbital periods close to one year. Similarly, we find the occurrence rate of larger planets increases substantially as orbital periods increase beyond $\sim~{}128$d. Therefore, we recommend that scientists and mission planners be cautious in adopting estimates of $\eta_{\oplus}$ that are primarily the result of extrapolating a parametric model into the $\eta_{\oplus}$ regime.

Instead, we recommend making use of the ABC-posterior distribution for $\eta_{\oplus}$ (or $\Gamma_{\oplus}$) which encapsulates the current direct observational constraints of the occurrence rate of small planets in the habitable zone of FGK stars. This represents a direct measurement based on Kepler data for the relevant planet sizes and orbital periods and not an extrapolation. Our constraints on $\eta_{\oplus}$ are only possible thanks to the Kepler project providing DR 25 data products for both actual Kepler data and simulated data with injected transits.

For purposes of planning missions to directly image potentially Earth-like planets, mission concepts are computing expected mission yields based on the ranges $P=$ $338-788$d and $R_{p}=$ $0.8-1.4R_{\oplus}$ for a solar twin. Unfortunately, Kepler’s mission duration and sensitivity do not place meaningful constraints on small planets at such long orbital periods. If one were to extrapolate our results using the Dirichlet prior, assuming the differential occurrence rate derived over the interval $P=237-500$d and $R_{p}=$ $0.75-1.5R_{\oplus}$, then one obtains an upper limit (95th percentile) of 0.33 planets per sun-like star with $P=$ $338-788$d and $R_{p}=$ $0.8-1.4R_{\oplus}$. The inclusion of long-period planets with sizes less than 1 R⊕, results in such an extrapolation being sensitive to the choice of prior. For example, adopting independent uniform priors for the rates of each radius bin results in the 95th percentile of the posterior increasing to 1.2 planets per sun-like star with $P=$ $338-788$d and $R_{p}=$ $0.8-1.4R_{\oplus}$. This demonstrates that such estimates inevitably are sensitive to the choice of prior and have large uncertainties. As emphasized above, we recommend considerable caution when adopting occurrence rates that are based on an extrapolation from shorter orbital periods.

4.3 Future Prospects for Occurrence Rate Studies

The model used in this study could be improved further to account for additional effects. We summarize the assumptions of our current model and our rationale for postponing consideration of these effects.

Detection Model: We model the probability that a planet is detected as a function of the expected effective SNR and the number of transits observed. In principle, the detection model could be further improved to incorporate additional information (e.g., transit duration, stellar properties, sky group).

Uncertainty in Contamination: This study models the effects of contamination due to known stars, but does not account for the uncertainty in the contamination. If only a fraction of the flux observed for a given Kepler target comes from the planet-hosting star, then the transit depth is diluted by the contaminating flux. While this has been considered in detail for individual systems (Hirsch et al., 2017; Teske et al., 2018), it is unclear how much contamination affects planet occurrence rate studies. In this study, we account for the dilution of transit depths due to contamination from known stars when computing the transit depth and SNR. We adopt the contamination provided for each target in the Kepler Input Catalog (KIC). Of course, there may be additional dilution from unknown sources, particularly bound multiple star systems. We mitigate the effects of unrecognized contamination from multiple star systems by focusing our analysis on a subset of the Kepler target stars which excludes targets with either astrometric data suggestive of binary companion or that appear to reside significantly above the main sequence. While binary stars with more disparate luminosities will survive our target selection process, they have limited effect on the planet occurrence rate. If a planet transits a star significantly fainter than the primary star, then the detection probability is significantly reduced by contamination from the primary. Therefore, such planets will rarely be detected for targets in our cleaned target list. If a planet transits the primary star in a system with a significantly fainter secondary star, then the effect of dilution is modest. Therefore, we expect that the size of the effect of unrecognized dilution will have a modest effect on the occurrence rates inferred using our cleaned target list.

Multiple Transiting Planets & TCEs: For targets with multiple transiting planets, the Kepler pipeline searched for planet candidates in an iterative fashion. After each “threshold crossing event” (TCE) was identified, the pipeline masked out a portion of the lightcurve during and immediately surrounding each putative transit before searching for additional planet candidates. The pipeline selected larger MES candidates first so that the presence of a larger and/or shorter-period planet causes the effective duty cycle to be reduced when searching for lower MES signals around the same star. Further complicating matters, the pipeline identified many TCEs that were later discarded as not being strong planet candidates. Nevertheless, all TCEs result in less data being available for detecting additional planets transiting that target (Zink et al., 2019). Fully accounting for these effects would be extremely complicated due to the complexity of the pipeline, and the available data products. Further, the measurement of planet occurrence rates would no longer decouple from the characterization of planetary system architectures. An initial investigation of this effect, suggested that the detectability of planets may decrease by $\simeq~{}5.5\%$ ($15.9\%$) for planets with orbital periods less (greater) than 200 days, if the system contains another transiting planet which the pipeline would detect first (Zink et al., 2019). The DR25 pipeline run resulted in 34,032 TCEs, far greater than the $\simeq~{}8,000$ KOIs (Twicken et al., 2016). Since our combined detection and vetting efficiency is based on fitting to the results of pixel-level transit injections into light curves that already include these TCEs, our detection efficiency already accounts for this, but in an average sense, rather than using knowledge of specific TCEs for each star.

The decrease in detectability due to additional planets only affects a small fraction of stars, since the ratio of planet candidates to target stars is $\simeq~{}0.03$ in our sample. If all planetary systems were strictly coplanar, then short-period planets would transit whenever the long-period planet transited. However, for orbital periods of 237-500 days, even very modest mutual inclinations ($\simeq~{}1^{\circ}$) cause the conditional transit probability to decrease substantially. Therefore, the impact of this effect on the overall planet occurrence rate will be diluted by at least half, even if this effect is significant for analyzing the relative abundance of systems with single and multiple transiting planets. In summary, we estimate that any increase in the occurrence rate of long-period planets due to interactions of the Kepler pipeline and multiple planet systems is likely to be less than $\simeq~{}8\%$

Pipeline Timeouts: The actual Kepler DR25 transit search did not search every target for transit signals as small as the typical 7.1-$\sigma$ threshold. The DR25 stellar catalog reports the “MES threshold” and “timeout” indicators for when the search for transits was halted prematurely. The DR25 stellar catalog data indicate that only 94% of targets were searched to MES of 7.1 (the standard criterion for triggering a TCE) across all durations. When applying this only to stars in our clean FGK catalog, 96% of targets were searched to full depth. However, of the 4% that were not searched to MES of 7.1, the MES threshold was still below 8 in 80% of the cases. In practice, the detection and vetting probability is so low for such cases that our assumption that the pipeline would have found detections between MES of 8 and 7.1 around these 3.2% of stars would not affect our results significantly. The remaining 0.8% of stars with somewhat higher MES thresholds are also not a significant source of systematic error.

Reliability: Our analysis assumes that all planet candidates (identified by DR25 pipeline and associated with target stars in our cleaned sample) were correctly identified as planet candidates, rather than being some sort of false positive (e.g., diluted eclipsing binary) or false alarm (e.g., astrophysical noise). Therefore, technically, the occurrence rates we infer should be regarded as the combined occurrence rate of true planets and any astrophysical false positives or instrumental false alarms that masquerade as planets within the given range of planet sizes and orbital periods. Fortunately, the overall reliability of the Kepler planet candidate sample is quite high (Thompson et al., 2018; Bryson et al., 2019). Nevertheless, validating any individual planet candidate at high reliability is difficult (Burke et al., 2019).

Bryson et al. (2019) quantified the impact of vetting efficiency, false positives, and potential false alarms on planet occurrence rate estimates. Their results indicate that it is important to consider vetting efficiency for a wide range of planet periods and the reliability of DR25 planet candidates is very high, except for orbital planets greater than 250 days, which likely suffer from some level of false alarms. The results from our study include a detailed model for vetting efficiency, so we expect the occurrence rates we report in Figures 2 & 3 and Table 2 are accurate planet occurrence rates for all but the bins with the longest orbital period. For the last bin ($P=256-500$d) and for habitable zone estimates (see Figures 5 & 6 and Table 3), we caution that the occurrence rates we report could be interpreted as upper limits, due to the potential for a significant fraction of long-period planet candidates being false alarms. Results of Thompson et al. (2018); Bryson et al. (2019) suggest that this is likely less than a factor of $\sim~{}2$.

For bins with a small number of detected planets (e.g., small, long-period planets, but also planets with very short-periods and very large planets), a few planet candidates being a false positive or false alarm could cause our inferred rate to be larger than the true planet rate. While this could be a large fractional error, it would not contribute significantly to the total number of planets integrated over a broad range of orbital periods and planet sizes. Reliability concerns are particularly acute at long orbital periods (e.g., our 256-500d bin), as shown in Thompson et al. (2018); Bryson et al. (2019). Complicated interactions of detector noise and the spacecraft orbital period likely contribute to some of the planet candidates at long orbital periods. We have attempted to mitigate reliability concerns by analyzing a cleaned set of target stars. Nevertheless, we regard the apparent increase in the occurrence rate of R=$1.5-6$ R⊕ planets beyond 256d as highly suggestive, but worthy of a more detailed investigation that includes a detailed treatment of the reliability of long-period planet candidates. We note that reliability concerns will only lead to revisions downward in the inferred rate. Therefore, our upper limits on the occurrence rates of small planets are robust to reliability concerns.

For periods of 128-256d and 256-500d, upper limits for the sum of occurrence rate over all radii exceeds one planet per star. Previous studies have shown that the fraction of stars with planets is significantly less than the average number of planets per star, based largely on the observed distribution of targets with multiple transiting planets. If the reliability of these long-period planet candidates is shown to be high, then future studies should also investigate the extent to which long-term orbital stability helps to inform the upper limit on the occurrence rate of long-period planets. A high occurrence rate of large-radii planets could dynamically restrict the occurrence rate of smaller radius planets.

4.3.1 Recommendations

In conclusion, we recommend that mission concepts aiming to characterize potentially rocky planets in or near the habitable zone of sun-like stars prepare compelling science programs that would be robust to the true rate of $0.75-1.5$ R⊕ planets with orbital periods in 237-500 days, $f_{R,P}\simeq$ $0.03-0.40$ or a differential rate of $(d^{2}f)/[d(\ln P)~{}d(\ln R_{p})]\simeq$ $0.06-0.76$. This larger range is based on starting from the 90% credible interval derived from our calculations, but has been expanded to allow for the possibility of the reliability of such planets being as low as $\simeq~{}50\%$ or for the rate to increase by up to $\simeq~{}10\%$ due to issues discussed above.

We thank the entire Kepler team for the many years of work that has proven so successful and was critical to this study. We thank an anonymous referee for useful suggestions for improving the analysis and manuscript. We thank Gijs Mulders for feedback and sharing updated values for the Mulders et al. (2018) vetting model which allowed us to make more direct the comparisons of occurrence rates based on the different vetting models. We thank Steve Bryson for discussions of the potential impact of false alarms and reliability. D.C.H, E.B.F., D.R., and K.A. acknowledge support from NASA Origins of Solar Systems grant # NNX14AI76G and Exoplanet Research Program grant # NNX15AE21. E.B.F acknowledges support from NASA Kepler Participating Scientist Program Cycle II grant # NNX14AN76G. D.C.H, E.B.F. acknowledge support from the Penn State Eberly College of Science and Department of Astronomy & Astrophysics, the Center for Exoplanets and Habitable Worlds and the Center for Astrostatistics. The citations in this paper have made use of NASA’s Astrophysics Data System Bibliographic Services. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work made use of the gaia-kepler.fun crossmatch database created by Megan Bedell. We acknowledge the Institute for CyberScience (http://ics.psu.edu/) at The Pennsylvania State University, including the CyberLAMP cluster supported by NSF grant MRI-1626251, for providing advanced computing resources and services that have contributed to the research results reported in this paper. This material was based upon work partially supported by the National Science Foundation under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute (SAMSI). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. This study benefited from the 2013 SAMSI workshop on Modern Statistical and Computational Methods for Analysis of Kepler Data, the 2016/2017 Program on Statistical, Mathematical and Computational Methods for Astronomy, and their associated working groups.

Appendix A ExoPAG Grid DR25 Occurrence Rates

Below in Table 4 we present our DR25 occurrence rates for F, G and K stellar types as defined by the Study Analysis Group (SAG) 13 of the NASA Exoplanet Exploration Program Analysis Group (ExoPAG). The combined detection and vetting efficiency is used. \startlongtable

\listofchanges