A Probabilistic Approach to Kepler Completeness and Reliability for Exoplanet Occurrence Rates

TL;DR

This paper introduces a probabilistic method to accurately estimate exoplanet occurrence rates from Kepler data by correcting for survey incompleteness and catalog reliability, significantly refining previous estimates.

Contribution

It presents a new probabilistic framework utilizing Kepler DR25 products and Gaia data to improve completeness and reliability corrections in exoplanet occurrence studies.

Findings

Correcting for reliability halves the estimated occurrence rate.

Using Gaia stellar properties significantly impacts occurrence rate estimates.

The method demonstrates the importance of accounting for catalog reliability.

Abstract

Exoplanet catalogs produced by surveys suffer from a lack of completeness (not every planet is detected) and less than perfect reliability (not every planet in the catalog is a true planet), particularly near the survey's detection limit. Exoplanet occurrence rate studies based on such a catalog must be corrected for completeness and reliability. The final Kepler data release, DR25, features a uniformly vetted planet candidate catalog and data products that facilitate corrections. We present a new probabilistic approach to the characterization of Kepler completeness and reliability, making full use of the Kepler DR25 products. We illustrate the impact of completeness and reliability corrections with a Poisson-likelihood occurrence rate method, using a recent stellar properties catalog that incorporates Gaia stellar radii and essentially uniform treatment of the stellar population. Correcting for reliability has a significant impact: the exoplanet occurrence rate for orbital period and radius within 20% of Earth's around GK dwarf stars, corrected for reliability, is 0.015+0.011-0.007, whereas not correcting results in 0.034+0.018-0.012 - correcting for reliability reduces this occurrence rate by more than a factor of two. We further show that using Gaia-based vs. DR25 stellar properties impacts the same occurrence rate by a factor of two. We critically examine the the DR25 catalog and the assumptions behind our occurrence rate method. We propose several ways in which confidence in both the Kepler catalog and occurrence rate calculations can be improved. This work provides an example of how the community can use the DR25 completeness and reliability products.

Tables

Table 1. Table 1: Baseline occurrence rate results

Parameter	No Reliability	With Reliability	FA-only Reliability
$F_0$	${0.608}_{-0.090}^{+0.110}$	${0.432}_{-0.072}^{+0.089}$	${0.514}_{-0.083}^{+0.102}$
$\alpha$	${0.304}_{-0.496}^{+0.519}$	${0.796}_{-0.598}^{+0.635}$	${0.500}_{-0.524}^{+0.558}$
$\beta$	$-{0.557}_{-0.169}^{+0.174}$	$-{0.823}_{-0.209}^{+0.202}$	$-{0.742}_{-0.196}^{+0.192}$
${\Gamma }_{\oplus }$	${0.212}_{-0.075}^{+0.111}$	${0.094}_{-0.041}^{+0.066}$	${0.139}_{-0.055}^{+0.086}$
$F_1$	${0.190}_{-0.030}^{+0.035}$	${0.144}_{-0.027}^{+0.032}$	${0.171}_{-0.029}^{+0.034}$
${\zeta }_{\oplus }$	${0.034}_{-0.012}^{+0.018}$	${0.015}_{-0.007}^{+0.011}$	${0.023}_{-0.009}^{+0.014}$

Table 2. Table 2: Impact of planet radius uncertainties (no reliability).

Parameter	No Uncertainty	Planet Radius Uncertainty
$F_0$	${0.608}_{-0.090}^{+0.110}$	${0.663}_{-0.112}^{+0.143}$
$\alpha$	${0.304}_{-0.496}^{+0.519}$	$-{0.172}_{-0.530}^{+0.553}$
$\beta$	$-{0.557}_{-0.169}^{+0.174}$	$-{0.593}_{-0.193}^{+0.189}$
${\Gamma }_{\oplus }$	${0.212}_{-0.075}^{+0.111}$	${0.276}_{-0.104}^{+0.157}$
$F_1$	${0.190}_{-0.030}^{+0.035}$	${0.216}_{-0.036}^{+0.041}$
${\zeta }_{\oplus }$	${0.034}_{-0.012}^{+0.018}$	${0.045}_{-0.017}^{+0.026}$

Table 3. Table 3: Impact of reliability uncertainty (with reliability).

Parameter	No Uncertainty	Reliability Uncertainty
$F_0$	${0.432}_{-0.072}^{+0.089}$	${0.428}_{-0.072}^{+0.090}$
$\alpha$	${0.796}_{-0.598}^{+0.635}$	${0.807}_{-0.592}^{+0.631}$
$\beta$	$-{0.823}_{-0.209}^{+0.202}$	$-{0.832}_{-0.214}^{+0.207}$
${\Gamma }_{\oplus }$	${0.094}_{-0.041}^{+0.066}$	${0.092}_{-0.040}^{+0.066}$
$F_1$	${0.144}_{-0.027}^{+0.032}$	${0.143}_{-0.027}^{+0.032}$
${\zeta }_{\oplus }$	${0.015}_{-0.007}^{+0.011}$	${0.015}_{-0.006}^{+0.011}$

Table 4. Table 4: Comparison of occurrence rate variations.

	$F_0$		$\alpha$		$\beta$
Case	No Reliability	With Reliability	No Reliability	With Reliability	No Reliability	With Reliability
Baseline (from Table 1)	${0.608}_{-0.090}^{+0.110}$	${0.432}_{-0.072}^{+0.089}$	${0.304}_{-0.496}^{+0.519}$	${0.796}_{-0.598}^{+0.635}$	$-{0.557}_{-0.169}^{+0.174}$	$-{0.823}_{-0.209}^{+0.202}$
DR25	${0.675}_{-0.097}^{+0.115}$	${0.474}_{-0.076}^{+0.090}$	$-{0.517}_{-0.389}^{+0.402}$	$-{0.339}_{-0.444}^{+0.465}$	$-{0.552}_{-0.155}^{+0.153}$	$-{0.888}_{-0.194}^{+0.188}$
Score $>0.9$	${0.418}_{-0.084}^{+0.112}$	${0.382}_{-0.079}^{+0.105}$	${0.616}_{-0.677}^{+0.708}$	${0.780}_{-0.724}^{+0.775}$	$-{0.774}_{-0.253}^{+0.244}$	$-{0.768}_{-0.263}^{+0.253}$
No Vetting Efficiency	${0.554}_{-0.084}^{+0.101}$	${0.389}_{-0.062}^{+0.079}$	${0.388}_{-0.502}^{+0.502}$	${0.889}_{-0.603}^{+0.637}$	$-{0.609}_{-0.175}^{+0.177}$	$-{0.889}_{-0.206}^{+0.203}$
No MES Smear	${0.632}_{-0.101}^{+0.125}$	${0.434}_{-0.073}^{+0.093}$	${0.104}_{-0.485}^{+0.508}$	${0.666}_{-0.593}^{+0.632}$	$-{0.527}_{-0.178}^{+0.171}$	$-{0.802}_{-0.209}^{+0.202}$
	${\Gamma }_{\oplus }$		$F_1$		${\zeta }_{\oplus }$
Case	No Reliability	With Reliability	No Reliability	With Reliability	No Reliability	With Reliability
Baseline (from Table 1)	${0.212}_{-0.075}^{+0.111}$	${0.094}_{-0.041}^{+0.066}$	${0.190}_{-0.030}^{+0.035}$	${0.144}_{-0.027}^{+0.032}$	${0.034}_{-0.012}^{+0.018}$	${0.015}_{-0.007}^{+0.011}$
DR25	${0.334}_{-0.098}^{+0.134}$	${0.164}_{-0.059}^{+0.084}$	${0.218}_{-0.029}^{+0.032}$	${0.174}_{-0.026}^{+0.028}$	${0.054}_{-0.016}^{+0.022}$	${0.027}_{-0.010}^{+0.014}$
Score $>0.9$	${0.103}_{-0.049}^{+0.089}$	${0.087}_{-0.044}^{+0.081}$	${0.139}_{-0.030}^{+0.038}$	${0.124}_{-0.029}^{+0.036}$	${0.017}_{-0.008}^{+0.014}$	${0.014}_{-0.007}^{+0.013}$
No Vetting Efficiency	${0.178}_{-0.064}^{+0.097}$	${0.076}_{-0.033}^{+0.054}$	${0.176}_{-0.028}^{+0.031}$	${0.132}_{-0.025}^{+0.030}$	${0.029}_{-0.010}^{+0.016}$	${0.012}_{-0.005}^{+0.009}$
No MES Smear	${0.246}_{-0.089}^{+0.132}$	${0.103}_{-0.044}^{+0.072}$	${0.198}_{-0.033}^{+0.036}$	${0.146}_{-0.028}^{+0.033}$	${0.040}_{-0.014}^{+0.021}$	${0.017}_{-0.007}^{+0.012}$

Table 5. Table 5: Model definitions

Model Name	$𝜽$	Model Definition	Description
constant	$\left[c\right]$	$c$	Constant function
gaussian	$[x_0,y_0,{\sigma }_x,{\sigma }_y,b]$	$AG(x,x_0,y,y_0,{\sigma }_x,{\sigma }_y)+b$	Gaussian + background
dualBrokenPowerLaw	$[b_x,b_y,{\alpha }_x,{\beta }_x,{\alpha }_y,{\beta }_y,A]$	$AB(x,b_x,{\alpha }_x,{\beta }_x)\times B(y,b_y,{\alpha }_y,{\beta }_y)$	Broken power law in $x$ and $y$
logisticY0	$[y_0,k,A]$	$AY(y,y_0,k,1)$	$y$ Logistic
logisticY	$[y_0,k,A,b]$	$AY(y,y_0,k,1)+b$	$y$ Logistic + constant
rotatedLogisticY	$[y_0,k,\varphi ,A,b]$	$AY(y_{\mathrm{rot}}+0.5,y_0,k,1)+b$	Rotated $y$ logistic + constant
rotatedLogisticX0	$[x_0,k,\varphi ,A]$	$AY(x_{\mathrm{rot}}+0.5,x_0,-k,1)$	Rotated $x$ logistic
rotatedLogisticX02	$[x_0,k,\varphi ,\nu ,A]$	$AY(x_{\mathrm{rot}}+0.5,x_0,-k,\nu )$	Rotated $x$ logistic w/ shape parameter
rotatedLogisticX0+gaussian	$[x_0,k,\varphi ,a_0,b_0,{\sigma }_x,{\sigma }_y,\gamma ,A]$	$AY(x_{\mathrm{rot}}+0.5,x_0,-k,1)+\gamma G(x,a_0,y,b_0,{\sigma }_x,{\sigma }_y)$	Rotated $x$ logistic + Gaussian
logisticX0xlogisticY0	$[x_0,y_0,k_x,k_y,A]$	$AY(x,x_0,-k,1)\times Y(y,y_0,k,1)$	$x$ logistic times $y$ logistic
logisticX0xlogisticY02	$[x_0,y_0,k_x,k_y,{\nu }_x,{\nu }_y,A]$	$AY(x,x_0,-k,{\nu }_x)\times Y(y,y_0,k,{\nu }_y)$	$x$ logistic times $y$ logistic w/ shape parameters
logisticX0xRotatedLogisticY0	$[x_0,y_0,k_x,k_y,\varphi ,A]$	$AY(x,x_0,-k,1)\times Y(y_{\mathrm{rot}}+0.5,y_0,k,1)$	$x$ logistic times rotated $y$ logistic
logisticX0xRotatedLogisticY02	$[x_0,y_0,k_x,k_y,\nu ,\varphi ,A]$	$AY(x,x_0,-k,1)\times Y(y_{\mathrm{rot}}+0.5,y_0,k,\nu )$	$x$ logistic times rotated $y$ logistic w/ shape parameter
rotatedLogisticYXLogisticY	$[y_0,y_1,k_0,k_1,\varphi ,A,b]$	$AY(y_{\mathrm{rot}}+0.5,y_0,k_0,1)\times Y(y,y_1,k_1,1)+b$	$y$ logistic times rotated $y$ logistic + constant
rotatedLogisticX0xlogisticY0	$[x_0,y_0,k_x,k_y,{\varphi }_x,{\varphi }_y,A]$	$AY(x_{\mathrm{rot}}+0.5,x_0,-k,1)\times Y(y_{\mathrm{rot}}+0.5,y_0,k,1)$	Rotated $x$ logistic times rotated $y$ logistic
rotatedLogisticX0xlogisticY02	$[x_0,y_0,k_x,k_y,{\nu }_x,{\nu }_y,{\varphi }_x,{\varphi }_y,A]$	$AY(x_{\mathrm{rot}}+0.5,x_0,-k,{\nu }_x)\times Y(y_{\mathrm{rot}}+0.5,y_0,k,{\nu }_y)$	Rotated $x$ logistic times rotated $y$ logistic w/ shape parameters
rotatedLogisticYXFixedLogisticY	$[y_0,k,\varphi ,A,b]$	$AY(y_{\mathrm{rot}}+0.5,y_0,k_0,1)\times Y(y,0.25,33.331,1)+b$	fixed $y$ logistic times rotated $y$ logistic + constant

Table 6. Table 6: Candidate vetting completeness rate functions

Model Name	Median AIC	Minimum AIC	Relative Likelihood
logisticY0	2819.03	2819.03	2.68e-137
dualBrokenPowerLaw	2246.87	2245.45	4.70e-13
logisticX0xlogisticY0	2223.80	2223.79	4.82e-08
logisticX0xlogisticY02	2227.20	2226.76	8.78e-09
logisticX0xRotatedLogisticY0	2190.10	2190.11	1.00
logisticX0xRotatedLogisticY02	2190.06	2189.83	1.02
rotatedLogisticX0xlogisticY0	2192.08	2192.11	0.372
rotatedLogisticX0xlogisticY02	2193.69	2193.19	0.166

Table 7. Table 7: Candidate false alarm effectiveness rate functions

Model Name	Median AIC	Minimum AIC	Relative Likelihood
rotatedLogisticX0	214.39	214.06	1
rotatedLogisticX02	212.49	210.93	1.70
constant	270.55	270.55	6.35e-13
dualBrokenPowerLaw	250.08	244.08	1.77e-08
gaussian	247.10	243.48	7.87e-08
rotatedLogisticX0xlogisticY0	220.44	220.08	4.84e-02
rotatedLogisticX0+gaussian	228.03	216.74	1.09e-03
rotatedLogisticY	215.80	212.36	2.35e-02
rotatedLogisticYXLogisticY	220.60	216.53	4.47e-02
logisticY	238.92	238.45	4.69e-06
rotatedLogisticYXFixedLogisticY	215.95	212.58	0.458

Table 8. Table 8: Candidate observed false alarm rate functions

Model Name	Median AIC	Minimum AIC	Relative Likelihood
rotatedLogisticX0	307.20	307.19	1.00
rotatedLogisticX0xlogisticY0	313.20	310.49	4.98e-02
dualBrokenPowerLaw	396.48	394.01	4.10e-20
rotatedLogisticX02	302.87	301.10	8.71
rotatedLogisticX0xlogisticY02	310.98	303.80	0.151
rotatedLogisticX0+gaussian	307.02	305.55	1.10

Table 9. Table 9: Comparison of occurrence rates using different catalogs and cuts.

	$F_0$		$F_1$		${\zeta }_{\oplus }$
Case	No Reliability	With Reliability	No Reliability	With Reliability	No Reliability	With Reliability
1	${0.608}_{-0.090}^{+0.110}$	${0.432}_{-0.072}^{+0.089}$	${0.190}_{-0.030}^{+0.035}$	${0.144}_{-0.027}^{+0.032}$	${0.034}_{-0.012}^{+0.018}$	${0.015}_{-0.007}^{+0.011}$
2	${0.609}_{-0.091}^{+0.112}$	${0.393}_{-0.065}^{+0.083}$	${0.186}_{-0.028}^{+0.031}$	${0.133}_{-0.024}^{+0.029}$	${0.040}_{-0.013}^{+0.019}$	${0.016}_{-0.007}^{+0.011}$
3	${0.680}_{-0.099}^{+0.121}$	${0.470}_{-0.079}^{+0.098}$	${0.216}_{-0.029}^{+0.032}$	${0.170}_{-0.026}^{+0.030}$	${0.056}_{-0.017}^{+0.024}$	${0.027}_{-0.010}^{+0.015}$
4	${0.675}_{-0.097}^{+0.115}$	${0.474}_{-0.076}^{+0.090}$	${0.218}_{-0.029}^{+0.032}$	${0.174}_{-0.026}^{+0.028}$	${0.054}_{-0.016}^{+0.022}$	${0.027}_{-0.010}^{+0.014}$
5	${0.678}_{-0.096}^{+0.117}$	${0.476}_{-0.077}^{+0.094}$	${0.219}_{-0.028}^{+0.031}$	${0.175}_{-0.026}^{+0.029}$	${0.055}_{-0.016}^{+0.022}$	${0.027}_{-0.010}^{+0.014}$

Equations

$$P(c_{i,j}|n_{i,j},\mathbf{\theta })=(\genfrac{}{}{0ex}{}{n_{i,j}}{c_{i,j}}){\rho }_{i,j}(\mathbf{\theta }{)}^{c_{i,j}}{\left(1-{\rho }_{i,j}(\mathbf{\theta })\right)}^{n_{i,j}-c_{i,j}}.$$

$$L(\mathbf{c},\mathbf{n},\mathbf{\theta })=\prod _{i,j}\left(\genfrac{}{}{0ex}{}{n_{i,j}}{c_{i,j}}\right){\rho }_{i,j}(\mathbf{\theta }{)}^{c_{i,j}}{\left(1-{\rho }_{i,j}(\mathbf{\theta })\right)}^{n_{i,j}-c_{i,j}}$$

$$\begin{split}\log(L)=&\sum_{i,j}\biggl{[}\log\binom{n_{i,j}}{c_{i,j}}+c_{i,j}\log\left(\rho_{i,j}(\boldsymbol{\theta})\right)\\ &+\left(n_{i,j}-c_{i,j}\right)\log\left(1-\rho_{i,j}(\boldsymbol{\theta})\right)\biggr{]}.\end{split}$$

$$Y\left(x,x_0,k,\nu \right)={\left[1+\exp (-k(x-x_0))\right]}^{-\frac{1}{\nu }}$$

$$\begin{array}{rl}x=& \frac{\left(p-p_{\min }\right)}{\left(p_{\max }-p_{\min }\right)}\\ y=& \frac{\left(m-m_{\min }\right)}{\left(m_{\max }-m_{\min }\right)}\\ y_{\mathrm{rot}}=& (y-0.5)\ast \cos (\varphi )-(x-0.5)\ast \sin (\varphi )\\ \rho =& AY\left(x,x_0,-k_x,1\right)\\ & \times Y\left(y_{\mathrm{rot}}+0.5,y_0,k_y,1\right).\end{array}$$

$$\begin{array}{rl}{x}_0& =1.257_{-0.044}^{+0.056},y_0=0.136_{-0.010}^{+0.009}\\ k_x& =4.311_{-0.503}^{+0.523},k_y=15.259_{-1.305}^{+1.399}\\ \varphi & =5.566_{-0.998}^{+1.106},A=0.980_{-0.005}^{+0.006}\end{array}$$

$$R\equiv \frac{T_{\mathrm{PC}}}{N_{\mathrm{PC}}}=1-\frac{N_{\mathrm{FP}}}{N_{\mathrm{PC}}}\left(\frac{1-E}{E}\right),$$

$$R_{\mathrm{FA}}=1-\frac{N_{\mathrm{FA}}}{N_{\mathrm{notFA}}}\left(\frac{1-E_{\mathrm{FA}}}{E_{\mathrm{FA}}}\right),$$

$$R_{\mathrm{FA}}=1-\frac{F_{\mathrm{FA}}}{1-F_{\mathrm{FA}}}\left(\frac{1-E_{\mathrm{FA}}}{E_{\mathrm{FA}}}\right).$$

$$\begin{array}{rl}x=& \frac{\left(p-p_{\min }\right)}{\left(p_{\max }-p_{\min }\right)}\\ y=& \frac{\left(m-m_{\min }\right)}{\left(m_{\max }-m_{\min }\right)}\\ x_{\mathrm{rot}}=& (x-0.5)\ast \cos (\varphi )-(y-0.5)\ast \sin (\varphi )\\ E_{\mathrm{FA}}=& AY\left(x_{\mathrm{rot}}+0.5,x_0,-k_x,1\right)\end{array}$$

$$\begin{array}{rl}{x}_0& =1.159_{-0.044}^{+0.062},k_x=22.587_{-6.291}^{+8.811},\\ \varphi & =98.551_{-2.778}^{+3.834},A=0.998_{-0.002}^{+0.001}.\end{array}$$

$$\begin{array}{rl}{x}_0& =0.682_{-0.029}^{+0.028},k_x=14.120_{-1.335}^{+1.469},\\ \varphi & =-157.967_{-3.539}^{+3.608},A=0.982_{-0.004}^{+0.004}.\end{array}$$

$$\begin{array}{rl}& P\left(\mathbf{\theta }|\{p_i,r_i\},i=1\dots N_p\right)\\ & \propto P\left(\{p_i,r_i\},i=1\dots N_p|\mathbf{\theta }\right)\pi \left(\mathbf{\theta }\right).\end{array}$$

$$\begin{array}{rl}& P\left(\{p_i,r_i\},i=1\dots N_p|\mathbf{\theta }\right)\\ & =e^{-\Lambda (D)}\prod _{s=1}^{N_p}\lambda (p_i,r_i,\mathbf{\theta })\end{array}$$

$$\begin{array}{r}\lambda (p,r,\mathbf{\theta })=F_0\frac{(\alpha +1)r^{\alpha }}{r_{\max }^{\alpha +1}-r_{\min }^{\alpha +1}}\frac{(\beta +1)p^{\beta }}{p_{\max }^{\beta +1}-p_{\min }^{\beta +1}}\end{array}$$

$$\begin{array}{rl}x=& \frac{\left(p-p_{\min }\right)}{\left(p_{\max }-p_{\min }\right)}\\ y=& \frac{\left(m-m_{\min }\right)}{\left(m_{\max }-m_{\min }\right)}\\ x_{\mathrm{rot}}=& (x-0.5)\ast \cos ({\varphi }_x)-(y-0.5)\ast \sin ({\varphi }_x)\\ y_{\mathrm{rot}}=& (y-0.5)\ast \cos ({\varphi }_y)-(x-0.5)\ast \sin ({\varphi }_y)\end{array}$$

$$B(x,b,\alpha ,\beta )=\{\begin{array}{ll}{\left(\frac{x+1}{b+1}\right)}^{\alpha }& x<b\\ {\left(\frac{x+1}{b+1}\right)}^{\beta }& x\ge b\end{array}G(x,x_0,y,y_0,{\sigma }_x,{\sigma }_y)=\exp \left(-\frac{(x-x_0{)}^2}{2{\sigma }_x^2}-\frac{(y-y_0{)}^2}{2{\sigma }_y^2}\right).$$

$$P\{N\left(B_i\right)=n_i\}=\frac{{\left(\Lambda (B_i)\right)}^{n_i}}{n_i!}{e}^{-\Lambda (B_i)}$$

$$\Lambda (B_i)={\int }_{B_i}\lambda (p,r)dpdr.$$

$$P\{N\left(B_i\right)=n_i\}\approx \{\begin{array}{ll}\lambda (p_i,r_i)\Delta p\Delta r{e}^{-\Lambda (B_i)}& n_i=1\\ e^{-\Lambda (B_i)}& n_i=0.\end{array}$$

$$\begin{array}{rl}& P\{N\left(B_i\right)=n_i,i=1,\dots ,K\}\\ & =\prod _{i=1}^K\frac{{\left(\Lambda (B_i)\right)}^{n_i}}{n_i!}{e}^{-\Lambda (B_i)}\\ & \approx {\left(\Delta p\Delta r\right)}^{K_1}{e}^{-{\sum }_{i=1}^K\Lambda (B_i)}\prod _{i=1}^{K_1}\lambda (p_i,r_i)\\ & ={\left(\Delta p\Delta r\right)}^{K_1}{e}^{-{\int }_D\lambda (p,r)dpdr}\prod _{i=1}^{K_1}\lambda (p_i,r_i)\end{array}$$

$$P\{N\left(B_i\right)=n_i,i=1,\dots ,K\}={\left(\Delta p\Delta r\right)}^{K_1}{e}^{-{\int }_D{\eta }_s(p,r)\lambda dp(p,r)dr}\prod _{i=1}^{K_1}{\eta }_s(p_i,r_i)\lambda (p_i,r_i).$$

$$\begin{array}{rl}P\{N_s\left(B_i\right)=n_{s,i},s=1,\dots ,N_{\ast },i=1,\dots ,K\}& =\prod _{s=1}^{N_{\ast }}{\left(\Delta p\Delta r\right)}^{K_1}{e}^{-{\int }_D{\eta }_s(p,r)\lambda (p,r)dpdr}\prod _{i=1}^{K_1}{\eta }_s(p_i,r_i)\lambda (p_i,r_i)\\ & =V{e}^{-{\int }_D\eta (p,r)\lambda (p,r)dpdr}\prod _{s=1}^{N_{\ast }}\prod _{i=1}^{K_1}{\eta }_s(p_i,r_i)\lambda (p_i,r_i)\end{array}$$

$$\begin{array}{rl}P\{N_s\left(B_i\right)=n_{s,i},s=1,\dots ,N_{\ast },i=1,\dots ,K|\mathbf{\theta }\}& =V{e}^{-{\int }_D\eta (p,r)\lambda (p,r,\mathbf{\theta })dpdr}\prod _{s=1}^{N_{\ast }}\prod _{i=1}^{K_1}{\eta }_s(p_i,r_i)\lambda (p_i,r_i,\mathbf{\theta })\\ & =V\left(\prod _{s=1}^{N_{\ast }}{\eta }_s(p_i,r_i)\right)e^{-{\int }_D\eta (p,r)\lambda (p,r,\mathbf{\theta })dpdr}\prod _{i=1}^{K_1}\lambda (p_i,r_i,\mathbf{\theta }).\end{array}$$

$$P\{N_s\left(B_i\right)=n_{s,i},s=1,\dots ,N_{\ast },i=1,\dots ,K|\mathbf{\theta }\}=e^{-{\int }_D\eta (p,r)\lambda (p,r,\mathbf{\theta })dpdr}\prod _{i=1}^{K_1}\lambda (p_i,r_i,\mathbf{\theta }).$$

Full text

Kepler — DR25 — exoplanets — exoplanet occurrence rates — catalogs — surveys

††journal: AAS Journals††facilities: Kepler††software: Jupyter, emcee (Foreman-Mackey et al., 2013), KeplerPorts (Burke & Catanzarite, 2017)

1 Introduction

The Kepler space telescope (Borucki et al., 2010; Koch et al., 2010) has delivered unique data that enables the characterization of exoplanet population statistics, from hot Jupiters in short-period orbits to terrestrial-size rocky planets in orbits with periods up to one year111https://exoplanetarchive.ipac.caltech.edu/docs/occurrence_rate_papers.html. By observing $>$150,000 stars nearly continuously for four years looking for transiting exoplanets, Kepler detected several thousand planet candidates (PCs) (Thompson et al., 2018), leading to the confirmation or statistical validation of over 2,300 exoplanets. This rich trove of exoplanet data has delivered many insights into exoplanet structure and formation, and promises deeper insights with further analysis. One of the most exciting insights to be gained from Kepler data is the occurrence rate of temperate, terrestrial-size planets orbiting Sun-like stars (often referred to as $\eta_{\oplus}$). This occurrence rate is also a critical input to the design of future space telescopes designed to discover and characterize habitable exoplanets, such as HabEx and LUVOIR.

Fully exploiting Kepler data requires a thorough understanding of how well it reflects the underlying exoplanet population. There are several ways in which the Kepler planet candidate catalog does not directly measure the real planet population:

• •

the catalog is incomplete, missing real planets

• •

it may be unreliable, with the planet candidate catalog being polluted with false positives

• •

it may be inaccurate due to observational errors leading to incorrect planet properties.

Lack of completeness and reliability are particularly acute at the Kepler detection limit, which happens to coincide with the period and radius of Earth-Sun analog exoplanets. We therefore focus our attention on a period and radius range spanning the Kepler detection limit.

In this paper we address vetting incompleteness, a significant component of incompleteness caused by incorrectly classifying detected true planets as false positives, and vetting reliability, caused by incorrectly classifying detections as planets candidate when they are in fact not true planets. We address accuracy by using new, uniformly determined stellar properties based in part on Gaia observations, described in §3.1.

We focus our analysis on the final Kepler data release DR25 (Thompson et al., 2018) and its associated planet candidate catalog222https://exoplanetarchive.ipac.caltech.edu. DR25 contains several products designed to support the characterization of the completeness and reliability of the DR25 planet candidate catalog. The primary contribution of this paper is a new probabilistic approach to using the DR25 completeness and reliability products to characterize vetting completeness and reliability. We illustrate the impact of completeness and reliability with standard occurrence rate computations, and examine the impact on occurrence rates due to changes in various assumptions. This is the first occurrence rate computation that fully uses the DR25 completeness and reliability products to characterize vetting reliability.

1.1 Previous Work

Kepler’s survey of the Cygnus field involved four years of data collection and another four years of pipeline development, data processing, and survey characterization, culminating in the final deliveries referred to as Data Release 25 (DR25). Incremental data deliveries enabled preliminary science investigations, and several occurrence rate studies were executed as the survey progressed. The simplest, first-look estimates used Gaussian cumulative distribution functions (CDF) as proxies for pipeline completeness (Borucki et al., 2011), restricted samples where completeness was assumed to be near unity (Howard et al., 2012), and linear approximations to a Gaussian CDF (Fressin et al., 2013). Lacking a full characterization of the Kepler pipeline, others employed independent detection pipelines, including injection and recovery experiments operating on flux light curves to quantify the detection completeness (Petigura et al., 2013; Foreman-Mackey et al., 2014; Dressing & Charbonneau, 2015). Hsu et al. (2018) performed an occurrence rate calculation using approximate Bayesian computation (ABC) and Zink & Hansen (2019) computed a habitable zone occurrence rate taking into account the effect of planet multiplicity on completeness.

The performance of the Kepler pipeline was characterized incrementally as more data were collected using transit injection and recovery operating on raw pixel fluxes as described in Section 2.1 (Christiansen et al., 2013, 2015, 2016; Christiansen, 2017). These early studies provided positive feedback to the Kepler pipeline whereby deficiencies were identified and improved upon (Twicken et al., 2016). Occurrence rate calculations using 16 out of 17 quarters of data and the associated pipeline completeness offered a benchmark computation for testing methodologies and comparing independent pipelines (Burke et al., 2015). Systematic errors were explored and the tallest tent poles were identified. Among these tall tent poles were two standouts: stellar property uncertainties and catalog reliability.

The first studies to include a treatment of catalog reliability focused on identifying astrophysical false positives either deterministically through follow-up observations (Santerne et al., 2012) or probabilistically via population synthesis (Morton & Johnson, 2011; Morton & Swift, 2014; Morton et al., 2016). And while most treatments culled or weighted the planet population, others sought to model both the planet and astrophysical sources as part of the planet occurrence estimation (Fressin et al., 2013), with Farr et al. (2014) applying a mixture model approach. These efforts used idealized models of the false positive and false alarm populations. An extremely useful astrophysical false positive probability statistic was developed by Morton et al. (2016).

The most significant effort to characterize the reliability of the Kepler planet candidate catalog to date is the final Kepler DR25 catalog paper Thompson et al. (2018). The DR25 catalog includes a Robovetter score which estimates the confidence with which the Robovetter vetted a TCE. Thompson et al. (2018) suggests that restricting occurrence rate studies to high-Robovetter-score planet candidates avoids the problem of low-reliability candidates (we test this approach in §6.4.2). Hsu et al. (2018) and Mulders et al. (2018) apply this high-score approach in their occurrence rate studies.

Burke et al. (2019) analyzes DR25 reliability using the inverted and scrambled data, approaching reliability characterization via kernel density estimation. The focus of Burke et al. (2019) is how reliability impacts statistical exoplanet validation.

1.2 This Paper

The primary purpose of this paper is to present a probabilistic characterization of Kepler vetting completeness (§4.2) and reliability against false alarms (§5.1) and show the impact on standard occurrence rate computations. This probabilistic analysis is robust against sparse data and resolves detailed structure of the dependence of vetting completeness and reliability on orbital period and transit signal strength. We explore the impact of using this characterization to correct occurrence rates based on Kepler data by performing a standard Poisson-likelihood-based occurrence rate, following Burke et al. (2015). Our characterization depends on the exoplanet population, which in turn depends on the parent stellar sample that is searched for planets. We restrict our analysis to GK dwarf stars.

Our method of accounting for completeness and reliability proceeds by executing the following steps:

• •

Select a subset of the target star population, which will be our parent population of stars that are searched for planets. We apply various cuts intended to select well-behaved and well-observed stars, and we restrict our analysis to GK dwarfs, as described in §3.1.

• •

Use the injected data to characterize vetting completeness, described in §4.2.

• •

Compute the summed detection completeness, incorporating vetting completeness, described in §4.1.

• •

Use observed, inverted and scrambled data to characterize false alarm reliability, described in §5.1.

• •

Assemble the collection of planet candidates, including computing the reliability of each candidate from the false alarm reliability and False Positive Probability.

• •

Compute the desired occurrence rates, presented in §6.

We perform our completeness and reliability analysis on the planet radius range $0.5\leq\mathrm{radius}\leq 15$ $R_{\oplus}$. We perform our completeness analysis on the period range $50\leq\mathrm{period}\leq 500$ days and reliability analysis on the period range $50\leq\mathrm{period}\leq 600$ days. Our occurrence rates are focused on illustrating the impact of vetting completeness and reliability where both are low, so our occurrence rates are analyzed for $50\leq\mathrm{period}\leq 400$ days and $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$.

Following an introduction to the details of completeness, reliability and the data products that support their computation in §2, this paper has four major parts: §3 assembles the stellar and planet catalogs we use in our analysis. We choose a stellar catalog that incorporates Gaia stellar radii and features an essentially uniform treatment of the parent stellar sample. §4 describes catalog completeness and describes our characterization of vetting completeness. §5 describes our characterization of catalog reliability. In §6 we perform our baseline occurrence rate computations corrected for vetting completeness, emphasizing the difference between correcting and not correcting for reliability. We then explore alternative occurrence rate calculations, including using an alternative stellar properties catalog, demonstrating the impact of not correcting for vetting completeness and restricting our analysis to planet candidates with Robovetter score $>0.9$.

Throughout this paper we present results with confidence intervals that are the 14th and 86th percentiles of posterior distributions resulting from MCMC analysis using fixed inputs. These confidence intervals do not account for uncertainties in the inputs. We address the issue of uncertainties in the inputs in §6.3.

All results reported in this paper were produced with Python code, mostly in the form of Python Jupyter notebooks, found at the paper GitHub site333https://github.com/stevepur/DR25-occurrence-public.

2 Completeness, Reliability and Occurrence Rates,

As described above, completeness has two components. Detection completeness is the fraction of true planets that are detected by the Kepler pipeline. Vetting completeness is the fraction of detected true planets that are correctly vetted as planet candidates. Vetting reliability is the fraction of vetted planet candidates that are true planets.

During catalog creation, the reliability of the planet candidate catalog is increased by detecting and removing false positives using a variety of tests. When the transit signal is weak it can be difficult for these tests to distinguish false positives from true planets, so maximizing reliability via stringent tests can cause true planets to be classified as false alarms, reducing vetting completeness. The DR25 catalog addressed this problem with uniform automated vetting via the Robovetter, using tests that were tuned to strike a balance between vetting completeness and reliability. This uniform automated vetting made it possible to vet synthetic and modified data sets designed to statistically mimic true planets and false positives, described in §2.1, in exactly the same way that the observed data was vetted. In this way the completeness and reliability of the planet candidate catalog can be measured and corrected in occurrence rates.

We distinguish two broad classes of phenomena that pollute the planet candidate catalog:

• •

Astrophysical False Positives, such as grazing or eclipsing binaries, which produce a planet-transit-like signal with a regular ephemeris in observed light curves that are not due to planetary transits. There has been extensive effort to identify and remove such false positives from Kepler catalogs (e.g., Morton & Johnson, 2011; Bryson et al., 2013; Fressin et al., 2013; Coughlin et al., 2014; Morton et al., 2016; Thompson et al., 2018), and for high SNR transits the resulting removal of astrophysical false positives is very effective. For low SNR, however, it is more difficult to distinguish astrophysical false positives from true planetary transits. In this paper we address astrophysical false positives via the probabilistic evaluation of Morton et al. (2016).

• •

False Alarms, which trigger a transit detection with a regular ephemeris, but are not due to regularly repeating astrophysical phenomena. The dominant source of false alarms in Kepler data is instrumental artifacts. There are two important classes of instrumental artifacts that have been identified as responsible for the overwhelming majority of false alarms at long periods: rolling bands and statistical and pixel fluctuations.

• –

Rolling bands (Van Cleve & Caldwell, 2009; Caldwell et al., 2010) are thermally dependent quasi-sinusoidal electrical signals in the output of the Kepler CCDs. As the Kepler telescope slowly changes its attitude relative to the Sun, different parts of the photometer are illuminated. While the thermal insulation of the telescope and electronics is very good, it is not perfect and the resulting thermal variations cause the rolling bands to slowly move across the CCDs, often introducing signals that look very much like transits. Because Kepler’s attitude is determined by its 372-day orbit, rolling bands often induce transit-like signals that repeat with a nearly regular ephemeris with an approximately 372-day period. This is the cause of the sharp peak of TCEs in the left panel of Figure 1. Rolling bands are highly focal plane position dependent: some Kepler CCD channels have much more severe rolling bands than others.

• –

Statistical and pixel fluctuations are independent, unrelated dips in light curves due to cosmic ray hits, single transits, and statistical fluctuations that trigger transit detections when they accidentally fall into a regular ephemeris. This class of false alarms becomes much more common for long-period ephemerides because they only require three or four events to fall on a regular ephemeris, and explains the broader “shoulders” of the tall peak in the left panel of Figure 1.

In addition, stellar variability triggers false alarm transit detections on a regular ephemeris, typically at short periods. False positives and false alarms are treated differently in our analysis, as described in §5.

In an ideal world, measuring planet occurrence rates from Kepler data would be simple — divide the number of detected planets by the number of observed stars, correcting for the geometrical probability of a planet transit. To get an accurate occurrence rate, however, this approach must be corrected for completeness and reliability. This correction can be large for Earth-Sun analog systems, which are at the Kepler detection limit where both completeness and reliability are very low. Specifically, instrumental false alarms are a significant source of false transit detections in the long-period, low signal-to-noise (SNR) region of most interest for habitable zone occurrence rate studies for G and K dwarf stars. In this regime, as shown in Figure 1, the number of instrumental false alarms is very large compared to the expected population of true exoplanet detections.

2.1 DR25 Vetting and Reliability Products

The DR25 planet candidate catalog (Thompson et al., 2018) contains 4034 identified planet candidates (PCs) out of 8054 Kepler Objects of Interest (KOIs). The KOIs were extracted from a catalog of 34,032 transit detections called threshold crossing events (TCEs), which are periodic transit-like events (as identified by a matched filter; Jenkins, 2002) that have a combined signal strength above a threshold (typically $7.1\sigma$). Identification of the PCs from the KOIs was performed by a fully automated Robovetter. The Robovetter applies a variety of tests to each TCE, many of which are based on the synthetic test datasets described below, and planet candidates are TCEs that pass all tests while following a logic tree. Such automated vetting (and transit detection) is critical for the production of a statistically uniform catalog that is amenable to statistical correction for completeness and reliability.

The DR25 completeness products are based on injected data — a ground-truth of transiting planets is obtained by injecting transit signals with specific characteristics on all observed stars at the pixel level (Christiansen, 2017). This data is then analyzed by the Kepler detection pipeline to produce a catalog of detections at the injected ephemerides called injected and recovered TCEs, which are then sent through the same Robovetter used to identify planet candidates. The fraction of injected transits that are recovered as TCEs measures detection completeness, while the fraction of recovered TCEs that are vetted as planet candidates measures vetting completeness. A large number of transits were also injected on a small number of target stars to measure the dependence of completeness on transit parameters and stellar properties. This data is used to create high-resolution, per-target-star detection contours described in §4.1, providing completeness for each target star as a function of planet orbital period and radius (Burke & Catanzarite, 2017).

The rate of false alarms, measured for the first time in DR25, are characterized by manipulating observed data so that it contains no true astrophysical transiting exoplanet signals, creating a ground-truth in which any TCE or vetted planet candidate is an instrumental false alarm. There are two basic manipulations that create the data used to characterize the rate of false alarms:

• •

Data inversion flips the light curves “upside down” so that true transiting signals increase in brightness and are therefore not identified as transits. This is believed to preserve the quasi-sinusoidal rolling bands described above. The distribution of TCEs detected in the inverted data reproduces well the sharp peak at a period of 372 days that is seen in the distribution of observed TCEs in the left panel of Figure 1.

• •

Data scrambling shuffles the Kepler observational quarters in a way that destroys the regular ephemeris of astrophysical transit signals, preventing their detection by the Kepler pipeline. While this also prevents the detection of the same false alarms that are detected in the original observed data, it is believed to preserve the statistics of detections due to statistical and pixel fluctuations. The distribution of TCEs detected in the scrambled data is very similar to the broad shoulder near periods of one year in the distribution of observed TCEs seen in the left panel of Figure 1. Three different shuffles of the Kepler data are available.

The DR25 Robovetter uses a number of metrics to identify instrumental false alarms, and the inverted and scrambled data sets were used to tune their pass/fail thresholds. For an extensive discussion, see Thompson et al. (2018).

For many Robovetter metrics, the distribution of values from the inverted/scrambled data overlaps the distribution of values in the observed data, making it difficult to distinguish false alarms from true planets, particularly at long period and low SNR. This results in low reliability in some parameter spaces of the Kepler planet candidate catalog, especially near the detection limit. An effort to characterize this catalog reliability is described in Section 4 of Thompson et al. (2018). The work presented in this paper is an attempt to improve on this characterization.

3 Input Catalogs

3.1 The Stellar Catalog

Our occurrence rate starts with a parent stellar population of GK dwarf stars that is searched for planetary transits. The properties of each star determine the likelihood that a transiting planet of a given size will be observed. The radius of the planet is derived from the fitted radius ratio and stellar radius. While the most accurate stellar properties for each star is desirable for understanding the properties of the transiting planet, a statistical occurrence rate requires the most uniform stellar properties possible. This is an issue for Kepler data because target stars with actual transit detections are much better characterized than most targets stars without transit detections. This can potentially lead to unknown biases in the estimated detection completeness of §4.2.

The stellar catalog associated with the DR25 exoplanet catalog, Q1-Q17 DR25 (with supplement)444 https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_stellar_docs.html, is based on heterogeneous observations, with some stars having properties derived from asteroseismic data, others from spectral data, and most from photometric data, all fitted to Dartmouth isochrones (Mathur & Huber, 2016; Mathur et al., 2017). Berger et al. (2018) combined the DR25 stellar catalog with Gaia parallaxes (Lindegren et al., 2018) to improve stellar radii, yielding an average radius precision of less than 10% for most Kepler stars. However, the adopted effective temperatures were still heterogeneous, and no revisions of other stellar properties such as mass and surface gravities were performed.

For our parent stellar population we use the stellar catalog of Berger et al. (2020), which extends Berger et al. (2018) by deriving a full set of stellar properties from isochrone fitting using broadband photometry, Gaia parallaxes and spectroscopic metallicities where available. This catalog is based on the homogeneous derivation of temperatures and luminosities, which previously have been the dominant sources of systematic errors in stellar (and thus planet) radii. These consistently fitted stellar properties provide more uniformly derived stellar radii over the entire parent population than the DR25 stellar properties. We recompute the 4-parameter stellar limb darkening model coefficients using these stellar properties via the table tableeq5.dat in Claret & Bloemen (2011), assuming a microturbulent velocity of 2 km s*-1*. In §6.4.1 we compare the resulting baseline occurrence rates with those using the DR25 stellar properties. We address the issue of possible bias against small planets in the Berger et al. (2020) catalog in §6.4.1 and Appendix C.

Because we require information such as observational completeness from the DR25 catalog and the binary flag from Berger et al. (2018) for each target star, we merge Berger et al. (2020), the DR25 stellar catalog (with supplement), and Berger et al. (2018), keeping only the 177,798 stars that are in all three catalogs. We remove possibly poorly characterized, binary and evolved stars using the following cuts:

• •

Remove stars with Berger et al. (2020) goodness of fit (iso_gof) $<0.99$ and Gaia Renormalized Unit Weight Error (RUWE, provided by Berger et al. (2020)) (Lindegren, 2018) $>1.2$, leaving 162,219 stars. iso_gof measures the quality of the Berger et al. (2020) isochrone fitting, and RUWE combines several Gaia goodness-of-fit metrics. RUWE is expected to have a Gaussian distribution (Lindegren, 2018) for single stars. Figure  2 shows the distribution of RUWE for the Berger et al. (2020) catalog, with a Gaussian fit to those stars with RUWE $<1.15$. Above RUWE $>1.15$ there is an apparent excess in RUWE, with that excess becoming dominant ($>75\%$ of stars) at RUWE $\approx 1.2$. An excess of RUWE over a Gaussian distribution is believed to be a strong indicator of stellar multiplicity. For example, Kraus et al. (in prep) finds that few Kepler target stars with RUWE $>1.2$ are single stars. We find that the RUWE Gaussian distribution has a slight magnitude dependence: for $g\leq 13$ the fitted Gaussian has the mode at 0.98, while for $g>13$ we find the mode at 1.01. We balance the loss of “good” stars against removing “bad” stars by choosing an RUWE cutoff of 1.2, in contrast to the cutoff of 1.4 discussed in the Gaia literature, e.g. Lindegren (2018).

• •

Remove stars that, according to Berger et al. (2018), are likely binaries (Bin flag = 1 or 3; we allow Bin = 2 because that indicates a nearby companion star found via high-resolution imaging, which was only performed on a subset of the target stars). This leaves 160,633 stars.

• •

Remove stars that have evolved off the main sequence, recomputing the Evol flag described in Berger et al. (2018) using the Berger et al. (2020) stellar properties. We use the evolstate package555http://ascl.net/1905.003 to determine the evolution state of each star using the isochrone-fitted $T_{\rm eff}$, radius and logg as inputs. We remove those stars with Evol $>0$, indicating that they are likely not main sequence dwarfs. After removing these stars 105,118 stars remain.

We then remove stars whose observations were not well suited for long-period transit searches (Burke et al., 2015; Burke & Catanzarite, 2017):

• •

Remove noisy targets identified in the KeplerPorts package666https://github.com/nasa/KeplerPORTs/blob/master/DR25_DEModel_NoisyTargetList.txt, leaving 103,626 stars.

• •

Remove stars with NaN limb darkening coefficients, leaving 103,371 stars.

• •

Remove stars with NaN observation duty cycle, leaving 102,909 stars.

• •

Remove stars with a decrease in observation duty cycle $>30\%$ due to data removal from other transits detected on this star, leaving 98,672 stars.

• •

Remove stars with observation duty cycle $<60\%$, leaving 95,335 stars.

• •

Remove stars with data span $<1000$ days, leaving 87,765 stars.

• •

Remove stars with the DR25 stellar properties table timeoutsumry flag $\neq 1$, leaving 82,371 stars. This flag = 1 indicates that the Kepler pipeline completed its transit search on this star before timing out.

Finally we select our GK population using the isochrone-fitted effective temperature as $3900\mathrm{K}\leq T_{\rm eff}<6000\mathrm{K}$, using the temperature limits of Pecaut & Mamajek (2013), leaving 57,015 stars in our parent GK dwarf population. The distribution of luminosities of these stars, computed as $R_{*}^{2}T_{\rm eff}^{4}$ in Solar units, is shown in Figure 3.

The largest remaining GK star has a stellar radius of $1.536R_{\odot}$. Burke & Catanzarite (2017) states that the per-star detection completeness described in §4.1 is invalid for stars of radius $>1.25R_{\odot}$ because completeness is characterized only for transit duration below the maximum 15 hours searched by the Kepler pipeline. We extend this to $R_{*}>1.35R_{\odot}$ because we are restricting our orbital period to 400 days, which keeps transit duration under 15-hours for our GK stellar population. We do not impose the stellar radius criterion in our baseline, instead opting for the physically motivated selection based on the Evol flag (though we use the $R_{*}>1.35R_{\odot}$ radius cut when using the DR25 stellar properties catalog in §6.4.1). Our baseline stellar population has 1,043 stars, or 1.83%, with $R_{*}>1.35R_{\odot}$. The maximum transit duration across our baseline population for a 400 day period and eccentricity = 0 is 14.85 hours. The distribution of transit durations for our baseline stellar population assuming a 400-day period and eccentricity = 0 is shown in Figure 4, which shows that our population gets close to, but does not exceed, the 15-hour duration limit.

3.2 The Planetary Catalog

Our planetary catalog is the Q1–Q17 DR25 Kepler Object of Interest (KOI) table at the exoplanet archive2 (Thompson et al., 2018), restricted to planet candidates (KOIs with koi_pdisposition = CANDIDATE) on stars in the catalog from Berger et al. (2020). We accept the CANDIDATE and FALSE POSITIVE dispositions resulting from the uniform robovetter run on the TCEs.

For our baseline case, we recompute the planet radii $R_{p}$ (in Earth radii) from the stellar radii $R_{*}$ (in Solar radii) in Berger et al. (2020) and the ratio of the planet radius to the stellar radius, $A=R_{p}/R_{*}$ from the koi_ror column of the KOI table, as $R_{p}=AR_{*}R_{\odot}/R_{\oplus}$ where $R_{\odot}$ ($R_{\oplus}$) is the Solar (Earth) radius. We compute the planet radius uncertainties $\sigma_{R_{p}}$ from the stellar radius uncertainties $\sigma_{R_{*}}$ and planet radius to the stellar radius ratio uncertainties $\sigma_{A}$ via standard propagation of uncertainties: $\sigma_{R_{p}}=\sqrt{\sigma_{A}^{2}R_{*}^{2}+A^{2}\sigma_{R_{*}}^{2}}R_{\odot}/R_{\oplus}$, where the upper and lower uncertainties are computed independently.

4 The Completeness Model

As described in §1, the set of planet candidates in the DR25 KOI catalog is not expected to be complete: particularly near the Kepler detection limit we expect that some transiting planets will be detected while others will be missed, and some of those detected will be mis-classified as false positives. Detection completeness is a measure of the fraction of true transiting planets that are detected. Vetting completeness is a measure of the fraction of detected true transiting planets that are correctly classified as planet candidates. We expect detection and vetting completeness to be functions of the orbital period and the signal to noise ratio (SNR), which in the Kepler data processing pipeline is measured as the Multiple Event Statistic (MES) (Jenkins, 2002). MES measures the combined significance of all observed transits in the de-trended, whitened light curve.

Detection and vetting completeness are both measured using the DR25 transit injection data products777 https://exoplanetarchive.ipac.caltech.edu/docs/KeplerSimulated.html (Christiansen, 2017). Christiansen et al. (2013, 2015, 2016) used these injection products to produce average detection curves as a function of MES for various stellar populations, marginalized over period. While these marginalized occurrence rate curves are convenient, the Poisson likelihood method we use for our occurrence rate, described in §6, works best with completeness provided as a function of both MES and period. We will use the star-by-star detection completeness model of Burke & Catanzarite (2017), provided for each star as a two-dimensional function of MES and period that accounts for each star’s detailed observational coverage. These completeness models are derived from a comprehensive database of $1.2\times 10^{8}$ transit injection and recovery trials, which we summarize in §4.1. In §4.2 we introduce a new probabilistic approach to modeling vetting completeness.

4.1 Combined Detection and Vetting Completeness

We use the characterization of Kepler detection completeness computed by a modified version of the KeplerPorts code base3. This software computes a completeness function $\eta_{s}\left(p,r\right)$ (not to be confused with $\eta_{\oplus}$) as a function of period $p$ and planet radius $r$ for each star $s$. The completeness function is described in detail in Burke & Catanzarite (2017). We briefly summarize the main steps for calculating the completeness function and describe the augmentations that incorporate vetting completeness.

The detection completeness calculation begins with estimating the MES expected for a given planet period and radius based on the stellar properties of the host. For each period and radius, a central crossing transit depth is estimated based on the stellar properties and limb darkening provided by the stellar catalog. The central crossing transit depth is converted into an expected MES by interpolating in the tabulated values of the one-sigma depth functions for each target (Burke & Catanzarite, 2017). The one-sigma depth function corresponds to the signal depth that results in a 1-sigma value for MES, and is a function of the planet period and the expected transit duration. The resulting expected MES is mapped to detection completeness based on analysis of the injected data. We treat this pipeline detection completeness estimate and the vetting completeness as independent. Thus the detection completeness is multiplied by the vetting completeness function $\rho(p,\mathrm{expectedMES},\boldsymbol{\theta})$ described in §4.2.

This produces the combined detection and vetting completeness for a central transit. The impact of non-central transits are accounted for through MES smearing, which convolves completeness with a distribution derived from analysis of the injected data. Completeness is then multiplied by the tabulated window function, which accounts for observational gaps for this star and the requirement of the Kepler pipeline of having at least three transit events for a detection, and the geometric transit probability assuming a uniform distribution of the cosine of orbital inclination angles.

The output is a collection of completeness functions $\eta_{s}\left(p,r\right)$, one for each star $s$ which includes detection completeness, vetting completeness and geometric transit probability. We sum these functions to create $\eta\left(p,r\right)=\sum_{s=1}^{N_{*}}\eta_{s}\left(p,r\right)$ where $N_{*}$ is the number of searched stars. The summed completeness $\eta\left(p,r\right)$ is used in the occurrence rate calculations in §6.

4.2 Vetting Completeness

Vetting completeness is the fraction of detected TCEs that were correctly vetted as PCs. This vetting is uniformly performed on both the observed TCEs and on the injected data TCEs, both resulting from the Kepler data analysis pipeline, with the DR25 Robovetter (Coughlin, 2017; Thompson et al., 2018) using the same thresholds in both cases. Because in the injected data every TCE is by definition a PC, vetting completeness is simply the fraction of injected on-target TCEs that were identified as PC by the Robovetter. We study the dependence of the injected vetting completeness on TCE period and expected MES by binning the detected injected TCEs on a regular grid. Our approach treats vetting completeness as a statistical property of a stellar population, analyzed separately for each choice of stellar population or stellar properties or other choices that may change the stellar or planet population. We present our vetting completeness analysis of the baseline GK star population in detail. Other cases described in §6.4 require independent analysis, which can be found in the htmlArchive folders on the paper GitHub site3.

Previous treatments of vetting completeness, e.g. Thompson et al. (2018), partitions the expected MES-orbital period plane into cells and divides the number of injected TCEs vetted as a PC by the total number of injected TCEs in each cell, which is an estimate of the vetting completeness in that cell. Mulders et al. (2018) does the same on a radius-orbital period plane, and Coughlin (2017) does the same with multiple parameter combinations (MES, period, planet radius, stellar radius, stellar temperature, and insolation flux). Using this data, one can estimate the dependence of the vetting completeness as a function of expected MES or planet radius and orbital period based on the measured fraction in each cell via, for example, $\chi^{2}$ fitting to a parametric model assuming a Gaussian likelihood, as done in Mulders et al. (2018). When there are many TCEs and many detections, this method can be expected to work well. Near the Kepler detection limit, however, there will be few TCEs and fewer correct PC dispositions, leading to strong gridding dependence and the possibility of adjacent cells having very different values. For example, if adjacent cells have only one TCE each and one is vetted as PC while the other is vetted as false positive, then these adjacent cells will have completeness 0 and 1. Cells with no TCEs require special treatment. Addressing these problems by requiring cells large enough to contain many TCEs can result in large cells smoothing out details of the vetting completeness’ functional dependence.

Rather than fitting a parametric model of a particular functional form using $\chi^{2}$ methodology, we take a probabilistic approach using a binomial likelihood that readily handles sparsely populated regions of parameter space. Specifically, we treat the injected TCEs as a collection, with a rate $\rho$ of (correctly) vetted PCs and a rate $(1-\rho)$ of (incorrectly) vetted FPs, and vetting by the Robovetter as draws from this collection. This is a classic binomial problem, in which the probability distribution of correctly vetting a PC depends on $\rho$ and the number of TCEs in the underlying collection. For example, if there is only one TCE in a cell with a $\rho=50\%$ probability of being vetted as a PC, the probability distribution of (correctly) vetting that TCE as a PC is extremely broad, with equal probability of PC or FP. Thus it is expected that such adjacent cells with single TCEs will have vetting completeness 0 and 1, with equal likelihood. Cells with no TCEs are gracefully handled because they have zero probability of a vetted PC. This allows the use of fine grids that can detect details of the functional dependence of vetting completeness.

By partitioning the expected MES – period space with a grid and computing $\rho$ in each grid cell, we can measure the dependence of $\rho$ on expected MES and period, inferring the function $\rho(p,m)$. This is what we do in the next section.

4.2.1 Vetting Completeness for the GK Baseline

Figure 5 shows the number of TCEs in each grid cell detected by the Kepler pipeline in the injected data at the correct ephemeris in the GK baseline stellar population. The injections were with expected MES between about 8 and 15, and period less than 500 days (see Christiansen, 2017, for details). Figure 6 shows the percentage of TCEs in each cell that were vetted as PC by the robovetter. Perfect completeness is 100%. We see that for high expected MES and period $<$ 200 days the completeness in each cell is typically near 100%, while for period $>$ 300 days and expected MES $<$ 15 the completeness drops off. We will characterize this behavior using a function $\rho(p,m,\boldsymbol{\theta})$ of planet period $p$ and expected MES $m$, where the exact form of $\rho$ is specified below and $\boldsymbol{\theta}$ is the vector of function parameters. Given a specific form for $\rho$, we infer $\boldsymbol{\theta}$ from the number of PCs that are correctly vetted as PC in each cell.

We conceptualize the determination of $\rho(p,m,\boldsymbol{\theta})$ as a binomial problem, thinking of each TCE as a draw from a population of PCs (= injected TCEs), which will be correctly vetted as PC by the robovetter with probability $\rho(p,m,\boldsymbol{\theta})$. If the number of TCEs in cell $(i,j)$ is $n_{i,j}$, then the probability of vetting $c_{i,j}$ TCEs correctly as PCs, given $\rho(p_{i,j},m_{i,j},\boldsymbol{\theta})\equiv\rho_{i,j}(\boldsymbol{\theta})$, is given by the binomial distribution

[TABLE]

For a particular choice of the functional form of $\rho$, equation (1) will be used to find the $\boldsymbol{\theta}$ that is most consistent with the number of PCs in each cell.

We will infer our rate function $\rho(p,m,\boldsymbol{\theta})$ via an MCMC Bayesian inference. We treat each grid cell as independent identically distributed binomial realizations, which leads to the likelihood

[TABLE]

where $\boldsymbol{n}=\{n_{i,j}\}$ is the set of the number of injected TCEs in each cell, and $\boldsymbol{c}=\{c_{i,j}\}$ is the set of the number of injected TCEs vetted as PC in cell $(i,j)$.

We perform the MCMC inference using the emcee package888https://emcee.readthedocs.io, which requires the log likelihood

[TABLE]

We considered several functional forms for $\rho(p,m,\boldsymbol{\theta})$, described in Appendix A. Figure 6 suggests a product of functions that are approximately, but not exactly, coordinate aligned. Qualitatively, Figure 6 also suggests that a generalized logistic function

[TABLE]

may be a good fit. We construct many, though not all, of the functional forms considered in this paper from this generalized Logistic function.

Appendix A describes our use of the Akiake Information Criterion (AIC) and other considerations to select the form of $\rho$ that best fits the data. In all cases, before applying the function we transform from (period, expected MES) coordinates to homogeneous coordinates on the unit square $[0,1]\times[0,1]$, which allows rotation. Of the functions we considered, we find that a product of a non-rotated simplified logistic function in period $p$ times a rotated logistic in $p$ and expected MES $m$ to best fit the data: for $\boldsymbol{\theta}=[x_{0},y_{0},k_{x},k_{y},\phi,A]$,

[TABLE]

We used the uniform priors $-1\leq x_{0},y_{0}\leq 2$, $10^{-4}<k_{x},k_{y}<10^{4}$, $-180<\phi<180$, $0<A<1$, and initialized $\boldsymbol{\theta}$ by minimizing $-\log(L)$ using the Python optimize package. Our MCMC computation used 100 walkers, and ran for 5000 steps after 5000 steps of burn-in. Figure 9 shows the resulting posteriors, giving $\boldsymbol{\theta}=[x_{0},y_{0},k_{x},k_{y},\phi,A]$ as

[TABLE]

where the central values are the posterior median and the $+$ and $-$ errors are for the 84th and 16th percentiles. We denote the median parameter vector by $\boldsymbol{\bar{\theta}}$. The correlations between parameters apparent in Figure 9 reflect the fact that the parameters of the logistic function are not fully independent.

Figure 10 shows the resulting rate function $\rho(p_{i},m_{j},\boldsymbol{\bar{\theta}})$ evaluated at the median of the posteriors $\boldsymbol{\bar{\theta}}$, along with the underlying rates for each grid cell. Figure 11 shows two example positions on the expected MES-orbital period plane, illustrating both the dependence of vetting completeness on these parameters as well as the spread of vetting completeness due to the posterior $\boldsymbol{\theta}$ distribution. We find that the approach described in this section is robust against changes in the grid. Changing the grid resolution does not significantly change the results, so long as the resolution is sufficient to resolve features in the underlying data.

Figure 7 shows the mean and standard deviation of 1000 realizations, drawn from the posterior $\boldsymbol{\theta}$ distribution, of the fraction of correctly vetted PCs for each cell. This figure should be compared with Figures 5 and 6. As expected, where the number of TCEs per cell is low in Figure 5, the standard deviation is high. Figure 8 shows the residual of observed data in Figure 6 from the mean rate in Figure 7 in units of the standard deviation shown in Figure 7. We see that while there are isolated large outliers, as well as larger residual values where the standard deviation is high, there is no indication of a bias in $\rho(p_{i},m_{j},\boldsymbol{\theta})$. Additional details, including further characterization of the quality of our fit of $\rho(p_{i},m_{j},\boldsymbol{\theta})$, are found in the htmlArchive folders on the paper GitHub site3.

5 The Reliability Model

In this section we characterize the reliability of planet candidates in the Q1–Q17 DR25 KOI catalog. We apply the probabilistic approach of §4.2 to the problem of characterizing the probability that a DR25 planet candidate is in fact a false alarm due to instrumental systematics, or some types of stellar variability. We rely on the Q1–Q17 DR25 False Positive Probabilities table at the NASA Exoplanet Archive2 (Morton et al., 2016) to provide the probability that the planet candidate is a false positive due to astrophysical signals that imitate transits. The final reliability for each planet candidate is the product of the false alarm reliability and the false positive probability.

5.1 Vetting Reliability

Thompson et al. (2018) defined reliability as the ratio of the number of PCs which are true exoplanets, $T_{\mathrm{PC}}$, to the number of observed planet candidates $N_{\mathrm{PC}}$:

[TABLE]

where $N_{\mathrm{FP}}$ is the number of observed false positives and $E\equiv\frac{N_{\mathrm{FP}}}{T_{\mathrm{FP}}}$ is the false positive effectiveness, defined as the number of identified FPs, $N_{\mathrm{FP}}$, divided by the number of true FPs, $T_{\mathrm{FP}}$. The second equality in equation (6) is exact when all quantities are from the same population, such as the observed data analyzed by the DR25 catalog. Unfortunately, the true PCs and false positives, $T_{\mathrm{PC}}$ and $T_{\mathrm{FP}}$, are unknown for the observed data. As explained in Thompson et al. (2018), however, we can use the inverted and scrambled data sets7 described in §1, which are designed so that every detection is, by definition, a false alarm.

Astrophysical transit-like events such as KOIs or eclipsing binaries can trigger detections in the inverted and scrambled data, compromising the use of these data to measure false alarms. This happens in two ways: 1) the transits and eclipses add signals unlike the false alarms we are trying to measure, and 2) the Robovetter is not tuned to detect and remove these kinds of signals. Thompson et al. (2018) describes how the lists of inverted and scrambled detections were cleaned of signals from known transiting systems in §2.3.3. Essentially, targets that are known binaries (Kirk et al., 2016) and known KOIs (Thompson et al., 2018) are removed from the list of detections, so they do not count as either a FP or a PC. For the inverted set, because self-lensing and heartbeat star binaries type events can produce signals that look like inverted transits, 54 targets with significant periodic signals were also removed from the list. The detections dropped from the inverted and scrambled data used in this study, as well as in Thompson et al. (2018), are collected in the files kplr_droplist_inv.txt and kplr_droplist_scr.txt* (one for each scrambled data set) on the paper GitHub site3. These stars are removed from the inverted/scrambled data before the analysis described in this section.

The inverted and scrambled data are designed to measure false alarms, not all false positives. Thus, we must take care to restrict the formula for reliability in equation (6) to the population of false alarms. This implies that equation (6) becomes

[TABLE]

where FA indicates “false alarm”. $N_{\mathrm{FA}}$ is the number of identified false alarms in the observed data (determined via the NTL flag in the KOI table) and $N_{\mathrm{notFA}}$ is the number of transit detections that are not vetted as false alarms. $E_{\mathrm{FA}}=\frac{N_{\mathrm{FA}}}{T_{\mathrm{FA}}}$ is the false alarm effectiveness, the fraction of true false alarms $T_{\mathrm{FA}}$ (assumed to be all detections in the inverted/scrambled data) that are vetted as false positives in the scrambled and inverted data. Equation (7) describes only false alarms measured by the inverted and scrambled data. As described in 5.3, we will multiply this reliability against instrumental false alarms with the reliability against astrophysical false alarms (constructed using the Q1-Q17 False Positive Probabilities).

In order to apply the probabilistic approach developed in §4.2, we rephrase the reliability formula in terms of rates rather than numbers by defining the fractions $F_{\mathrm{FA}}=N_{\mathrm{FA}}/N_{\mathrm{TCE}}$ and $F_{\mathrm{notFA}}=N_{\mathrm{notFA}}/N_{\mathrm{TCE}}$, where $N_{\mathrm{TCE}}$ is the number of TCEs in the observed data. We call $F_{\mathrm{FA}}$ the observed false positive rate. We can substitute $N_{\mathrm{FA}}/N_{\mathrm{notFA}}=F_{\mathrm{FA}}/F_{\mathrm{notFA}}$ in equation (7). We further assume that notFA and FA are a complete partition of all TCEs in the observed data, so $N_{\mathrm{FA}}+N_{\mathrm{notFA}}=N_{\mathrm{TCE}}\Rightarrow F_{\mathrm{FA}}+F_{\mathrm{notFA}}=1$ and we can eliminate $N_{\mathrm{notFA}}$ from equation (7):

[TABLE]

We can now treat $F_{\mathrm{FA}}$ and $E_{\mathrm{FA}}$ as rates that determine the probability of drawing a TCE that will be vetted as a false alarm. As in §4.2 we will treat $F_{\mathrm{FA}}$ and $E_{\mathrm{FA}}$ as rates in two separate binomial problems, with functional forms that depend on period and observed MES and whose coefficients are determined via an MCMC inference. Additional details, including further characterization of the quality of our rate fits, are found in the htmlArchive folders on the paper GitHub site3.

5.1.1 Characterization of False Alarm Effectiveness $E_{\mathrm{FA}}$

To determine the false alarm effectiveness, $E_{\mathrm{FA}}$, we combine the inverted data with each of the three scrambled data sets (see Coughlin, 2017, for details on the properties of each set) to create three data sets, called “inverted/scrambled”, where every TCE should be considered a false alarm. We proceed as in §4.2, covering the period-(observed) MES plane with a regular grid, and measure the ratio of the number of false alarms to the number of TCEs in each cell. This problem is more challenging than the analysis of vetting completeness in §4.2 because the Robovetter has been tuned to do a very good job correctly identifying false alarms, resulting in relatively few cells with TCEs incorrectly vetted as PCs. We therefore combine the three inverted/scrambled data sets through concatenation in order to produce a somewhat stronger signal. This amounts to averaging the three data sets at the input level, avoiding small-number statistics issues that would arise if we fit the three data sets separately and averaging the resulting posteriors. We refer to this concatenated data set as the combined inverted/scrambled data.

Figure 12 shows the number of TCEs detected in the combined inverted/scrambled data. We see that most detected TCEs in this data is for MES $<15$ and period $\geq$ 250 days. Figure 13 shows the fraction of correctly vetted false alarms, a measure of $E_{\mathrm{FA}}$, in each cell. The signal we’re measuring is small, and is dominated by smaller fractions at low MES and period $\leq 200$ days.

We use the same likelihood as in §4.2, equation (6), where in this case $E_{\mathrm{FA}}$ plays the role of $\rho$, $n_{i,j}$ is the number of TCEs detected in the combined inverted/scrambled data in cell $(i,j)$, and $c_{i,j}$ is the number of false positives identified in cell $(i,j)$. We perform the MCMC inference as described in §4.2. We considered several functions, described in Appendix A.2, and determined that a simple rotated logistic function best describes this data set. For $\boldsymbol{\theta}=[x_{0},k_{x},\phi,A]$, $E_{\mathrm{FA}}\left(p,m,\boldsymbol{\theta}\right)$ is given by

[TABLE]

where $p$ is the orbital period, $m$ is the observed MES, and $Y$ is the logistic function from equation (4). We used the uniform priors $-1\leq x_{0}\leq 2$, $10^{-4}<k_{x}<100$, $-180<\phi<180$, $0<A<1$. The MCMC run used a hand-tuned initial condition because Python’s optimize maximum likelihood solution was physically unreasonable ($A>>1$, for example) and violated the prior. Our MCMC computation used 100 walkers, and ran for 5000 steps after 5000 steps of burn-in. Figure 14 shows the resulting posteriors, giving $\boldsymbol{\theta}=[x_{0},k_{x},\phi,A]$ as

[TABLE]

The rate function $E_{\mathrm{FA}}(p_{i},m_{j},\boldsymbol{\bar{\theta}})$ for the posterior median $\boldsymbol{\bar{\theta}}$ is shown in Figure 15. As in §4.2, 1000 realizations of the FP rate function were created, drawing from the posterior $\boldsymbol{\theta}$ distribution. The residuals of the observed false alarm fraction in Figure 13 from the mean of these realizations in units of standard deviation is shown in Figure 16, demonstrating an overall good fit to the data.

5.1.2 Characterization of Observed False Positive Rate $F_{\mathrm{FA}}$

For $F_{\mathrm{FA}}$ we count the number of false positives found in the observed data, but we must be careful to consider only non-transit-like false alarms in order to be consistent with our characterization of effectiveness. We identify such false alarms by selecting on the not-transit-like (NTL) flag = 0, indicating that the Robovetter identified this false positive as transit like, which identifies 21 candidate astrophysical false positives in our GK population inside $50\leq\mathrm{period}\leq 600$ days and $0.5\leq\mathrm{radius}\leq 15$ $R_{\oplus}$. We manually examined these FPs and identified those that show a consistent astrophysical signal in all transits as astrophysical false positives. Two TCEs with NTL=0 did not show such a consistent astrophysical signal and are deemed likely false alarms: 004371172-01 and 009394762-01. The other 19 FPs with NTL=0 were identified as astrophysical and removed from the set of FPs used in the analysis of $F_{\mathrm{FA}}$.

Figure 17 shows the number of TCEs detected in the combined inverted/scrambled data. We see that most detected TCEs in this data is for MES $<15$ and period $\geq 250$ days. The close correspondence with Figure 12 shows that most TCEs are false alarms. Figure 18 shows the fraction of identified false alarms (identified via the NTL flag as described above), a measure of $F_{\mathrm{FA}}$, in each cell.

We proceed in a very similar manner to inferring $E_{\mathrm{FA}}$ in §5.1.1. We use equation (6) as the likelihood, with $F_{\mathrm{FA}}$ playing the role of $\rho$, $n_{i,j}$ is the number of TCEs detected in the observed data in cell $(i,j)$, and $c_{i,j}$ is the number of false alarms identified in cell $(i,j)$. We perform the MCMC inference as described in §4.2. We considered several functions, described in Appendix A.3, and determined that the same simple rotated logistic function as that used in §5.1.1 best describes this data set. In equation (9), $F_{\mathrm{FA}}$ replaces $E_{\mathrm{FA}}$, providing $F_{\mathrm{FA}}\left(p,m,\boldsymbol{\theta}\right)$ for $\boldsymbol{\theta}=[x_{0},k_{x},\phi,A]$.

We used the uniform priors $-1\leq x_{0}\leq 2$, $10^{-4}<k_{x}<100$, $-180<\phi<180$, $0<A<1$, and initialized $\boldsymbol{\theta}$ by minimizing $-\log(L)$ using the Python optimize package. Our MCMC computation used 100 walkers, and ran for 5000 steps after 5000 steps of burn-in. Figure 19 shows the resulting posteriors, giving $\boldsymbol{\theta}=[x_{0},k_{x},\phi,A]$ as

[TABLE]

The rate function $F_{\mathrm{FA}}(p_{i},m_{j},\boldsymbol{\bar{\theta}})$ for the posterior median is shown in Figure 20. As in §4.2, 1000 realizations of the FP rate function were created, drawing from the posterior $\boldsymbol{\theta}$ distribution. The residuals of the observed false alarm fraction in Figure 18 from the mean of these realizations in units of standard deviation is shown in Figure 21, demonstrating an overall reasonable fit to the data.

5.1.3 Computing the False Alarm Reliability $R_{\mathrm{FA}}$

Once we have the rate functions $F_{\mathrm{FA}}$ and $E_{\mathrm{FA}}$, we can compute the false alarm reliability $R_{\mathrm{FA}}\left(p,m\right)$ from equation (8). In practice we evaluate $F_{\mathrm{FA}}$ and $E_{\mathrm{FA}}$ at a desired period and observed MES, either on a regular grid or for specific planet candidates.

Figure 22 shows the resulting reliability function in the period-MES plane. We see that for low MES there is decreased reliability around period 250 to 450 days, corresponding to the high number of TCEs in that range found in the inverted/scrambled data (see Figure 12), consistent with the excess of detections in Figure 1. Figure 23 shows the reliability function evaluated over the full posteriors of $F_{\mathrm{FA}}$ and $E_{\mathrm{FA}}$ for three example periods and observed MES. We see that for low MES near 1-year orbital periods the reliability drops to about 0.6 and has a large spread.

5.2 Astrophysical Reliability

The reliability function determined in §5.1 only provides the probability that a planet candidate is not a false alarm. To determine the probability that a candidate is not an astrophysical false alarm such as a grazing or eclipsing eclipsing binary, we use the Q1–Q17 DR25 False Positive Probabilities2 created using the technique developed in (Morton et al., 2016). These probabilities were computed for all KOIs based largely on photometric data including transit light curves and measured magnitudes. We therefore assume that they are still valid even though we are using different stellar properties. We define the astrophysical reliability of a planet candidate as 1 - the false positive probability of that candidate.

5.3 Computing the Reliability for Each Planet Candidate

We compute the reliability for each planet candidate by first evaluating $R_{\mathrm{FA}}\left(p,m\right)$ as described in §5.1.3, where $p$ is the observed orbital period and $m$ is the observed MES of the planet candidate from the KOI catalog. Then we define the reliability $R=R_{\mathrm{FA}}\left(p,m\right)\cdot\left(1-\mathrm{FPP}\right)$ where FPP is the false positive probability for that planet candidate from the Q1–Q17 False Positive Probabilities table.

6 Illustrative Occurrence Rates

We present several illustrative occurrence rates, focusing on long-period, small planets where vetting completeness and reliability have the greatest impact. We compute our occurrence rates with the method of Burke et al. (2015), modeling occurrence rates as a Poisson point process with a rate given by a product of power laws in orbital period and planet radius. We perform our occurrence rate analysis over the period and radius range of $50\leq\mathrm{period}\leq 400$ days and $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$, and integrate the resulting rate over two ranges considered by Burke et al. (2015):

• •

$F_{1}$ : $50\leq\mathrm{period}\leq 200$ days and $1\leq\mathrm{radius}\leq 2$ $R_{\oplus}$, and

• •

$\zeta_{\oplus}$ : within 20% of Earth’s orbital period and radius.

Figure 24 shows our baseline planet candidate population, with the planet markers sized and colored by that planet’s reliability and the background and contours showing the completeness function $\eta(p,r)$, including geometric transit probability. The $F_{1}$ and $\zeta_{\oplus}$ regions are indicated by boxes. We see that while the $F_{1}$ region is reasonably well populated, it has a large completeness correction of $\sim 500$. $\zeta_{\oplus}$, however, has only one low-reliability planet and a completeness correction $>10^{4}$, leading to large uncertainties in the estimate of $\zeta_{\oplus}$.

6.1 Methodology

Following Youdin (2011) and Burke et al. (2015), we study the number of planets per star as a function of orbital period $p$ and planet radius $r$, $f(p,r)$, by inferring the population rate function $\lambda(p,r)\equiv\mathrm{d}^{2}f/\mathrm{d}p\,\mathrm{d}r$ from a collection of planet detections at $(p_{i},r_{i})$ with a known completeness function $\eta(p,r)$ and reliability $R_{\mathrm{FA}}$. If $\lambda(p,r,\boldsymbol{\theta})$ is a specific function parameterized by the parameter vector $\boldsymbol{\theta}$, then our problem is to determine $\boldsymbol{\theta}$. We proceed by Bayesian inference: given a set of planet candidates with orbital period and radius $\{p_{i},r_{i}\}$, $i=1\ldots N_{p}$ where $N_{p}$ is the number of planet candidates, by Bayes’ theorem the probability of $\boldsymbol{\theta}$ is

[TABLE]

where $\pi\left(\boldsymbol{\theta}\right)$ is a prior on $\boldsymbol{\theta}$. Fixing $\pi\left(\boldsymbol{\theta}\right)$, finding the highest probability $\boldsymbol{\theta}$ amounts to maximizing the likelihood $P\left(\{p_{i},r_{i}\},i=1\ldots N_{p}|\boldsymbol{\theta}\right)$.

In Appendix B we show that maximizing the likelihood

[TABLE]

is equivalent to treating planet occurrence as a Poisson point process with rate $\lambda(p,r,\boldsymbol{\theta})$ that depends on period, planet radius, and parameters $\boldsymbol{\theta}$. Here $\Lambda(D)=\int_{D}\eta(p,r)\lambda(p,r,\boldsymbol{\theta})dp\,dr$ is the integral over the whole period-radius space $D$, where $\eta(p,r)$ is the summed completeness function from §4. However, we point out that equation (11) is not itself a Poisson probability, as is sometimes implied in the literature.

The likelihood in equation 11 accounts for completeness but not reliability. Because this likelihood is derived from a Poisson distribution, which is defined only for discrete integer counts, we cannot account for reliability by weighting a planet’s contribution by its reliability. We address reliability by performing multiple Bayesian inferences of $\boldsymbol{\theta}$ using equation (10), drawing from the planet candidates according to their reliability. For example, a planet candidate with reliability 0.9 would be included in 90% of these inferences, while another planet candidate with reliability 0.2 would be included in 20% of these inferences. Then the $\boldsymbol{\theta}$ posteriors of these inferences is concatenated to produce the posterior distribution of $\boldsymbol{\theta}$ accounting for reliability.

Following Youdin (2011) and Burke et al. (2015), we model the planet candidate population rate $\lambda(p,r,\boldsymbol{\theta})$ as a product of power laws in period and radius. Inspired by Foreman-Mackey’s implementation of Burke et al. (2015)999https://dfm.io/posts/exopop/, we adapt the form resulting from solving explicitly for the normalization $C_{n}$ from Burke’s equation (8) and using it in his unbroken power law equation (7) (Burke et al., 2015): for $\boldsymbol{\theta}=\left(F_{0},\alpha,\beta\right)$,

[TABLE]

This form ensures that $\int_{D}\lambda(p,r,\boldsymbol{\theta})dp\,dr=F_{0}$ so $F_{0}$ can be interpreted as the integrated planetary occurrence rate over the period-radius range used in the analysis.

6.2 Baseline Results

To perform our Bayesian MCMC inference, we use the emcee package. To measure the impact of correcting for reliability, we run inferences both without and with reliability correction. For our inference without reliability correction, we use 16 walkers and run for 5000 steps after 1000 steps of burn-in. For our inference with reliability correction, we run 100 inferences as described in §6.1, probabilistically sampling from the planet candidates according to their reliability, with each inference using 16 walkers and running for 2000 steps after 400 steps of burn-in. In both cases the walkers of each MCMC run are initialized in a small Gaussian distribution centered on the maximum-likelihood solution for that inference’s planet population. The posteriors from each of the 100 inferences with reliability correction were concatenated to produce the $\boldsymbol{\theta}$ posteriors. Table 1 shows the median and 16th and 84th percentiles of these posterior distributions both with and without reliability correction. We see that reliability has an overall impact of about 30% in $F_{0}$, the integrated rate over our period and radius range of $50\leq\mathrm{period}\leq 400$ days and $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$.

Figure 27 shows the marginalized population rate function $\lambda(p,r,\boldsymbol{\theta})$ for the posterior $\boldsymbol{\theta}$ distribution, accounting for uncertainty. This figure also compares the predicted number of planet detections with the binned planet candidates.

Figure 26 and Table 1 show $F_{1}$ and $\zeta_{\oplus}$, as well as $\Gamma_{\oplus}\equiv\mathrm{d}^{2}f/\mathrm{d}\log p\,\mathrm{d}\log r=p_{\oplus}r_{\oplus}\lambda\left(p_{\oplus},r_{\oplus},\boldsymbol{\theta}\right)$, with and without accounting for reliability, evaluated over all posterior values of $\boldsymbol{\theta}$. We see that even though there is significant overlap in the distributions with and without reliability, accounting for reliability has a strong impact: $\Gamma_{\oplus}$ and $\zeta_{\oplus}$ are are reduced by more than 50%, which can be understood from the very small number of low-reliability planets in the $\zeta_{\oplus}$ region in Figure 24. $F_{1}$ is the integrated rate over a region of higher reliability, but reliability still has a strong effect. $F_{0}$ is the integrated rate over our entire period-radius analysis range, but it is dominated by the fact that there are more high-reliability planet candidates, so reliability has an impact similar to $F_{1}$. Table 1 also shows the impact of accounting only for false alarm reliability, ignoring astrophysical false positive reliability, indicating that false alarm reliability accounts for about half the impact of the reliability correction.

We also computed occurrence for the SAG13 definition of $\eta_{\oplus}$101010https://exoplanets.nasa.gov/exep/exopag/sag/#sag13, $237\leq\mathrm{period}\leq 860$ days and $0.5\leq\mathrm{radius}\leq 1.5$ $R_{\oplus}$. Without accounting for reliability, we find $\eta_{\oplus}=0.302^{+0.181}_{-0.113}$, consistent with the results of Zink & Hansen (2019), while accounting for reliability yields $\eta_{\oplus}=0.126^{+0.095}_{-0.055}$. This result should be treated with caution because it involves extrapolation beyond the domain of both reliability and detection completeness characterization.

We find that the impact of accounting for reliability is significant for small planets in long-period orbits. While one can note that the median values of occurrence rates in this regime are not much more than “one $\sigma$” apart, the observed shifts in the distributions on the order of 40% are systematic, and clearly not due to statistical fluctuations.

6.3 Simple Estimates of the Impact of Input Uncertainty

A full treatment of uncertainties in occurrence rates is beyond the scope of this paper. Uncertainties in stellar properties would need to be accounted for in the selection of the parent stellar population, the modeling behind the detection completeness and impact on the Robovetter. In this work we do, however, produce uncertainties in the false alarm reliability in §5.1 through the MCMC posteriors of the fit functions, as well as planet radius uncertainties that follow from stellar radius transit fit uncertainties as described in §3.2. In this section we present simple experiments that examine the impact on our occurrence rates of the reliability and planet radius uncertainties. We study the impact of planet radius uncertainties separately from the impact of reliability uncertainty.

6.3.1 Impact of Planet Radius Uncertainty

We study the impact of planet radius uncertainties without accounting for reliability. We proceed in the same way that we study the impact of reliability, by performing several inference runs with a planet population in each run selected after applying the planet radius uncertainties. Specifically, for each run, prior to the restriction of the planet candidate population to the radius range $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$, we add to each planet’s radius an error given by a draw from a Gaussian distribution with width equal to that planet’s radius uncertainty. Each planet is randomly assigned an upper or lower errorbar with 50% probability. The planet candidate population is then restricted to the range $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$, and the inference is run.

The impact of planet radius uncertainties resulting from 1000 inference runs is shown in Table 2 and the top row of Figure 28. We see a small, consistent broadening of the width of the distributions and resulting error bars, and possibly a small systematic shift towards higher occurrence rates, but the overall impact is minor. We believe this is due to the smaller uncertainties resulting from using Gaia stellar properties, and from the fact that near the boundary of our planet size range there are many planets, so planets are equally likely to exit and enter our range due to uncertainty.

6.3.2 Impact of Reliability Uncertainty

We study the impact of uncertainties in reliability by modifying the method of computing occurrence rates with uncertainty described in §6.2. Prior to each inference run, we draw from the posteriors of the parameter vectors for false alarm efficiency (§5.1.1) and the observed false alarm rate (§5.1.2). We use these draws evaluate the FA efficiency and observed rate functions at each planet candidate’s period and observed MES, from which the false alarm reliability is computed. The reliability is then computed as described in §5.3, and each planet is included in that run with probability given by this computed reliability. In other words, the reliability rate function is realized for each run and applied to the planet candidate population.

Table 3 and the bottom row of Figures 28 compare occurrence rates with reliability correction but no reliability uncertainty, computed with the median $\boldsymbol{\bar{\theta}}$ of §5.1.1 and §5.1.2, with the reliability distribution that results from using the respective full $\boldsymbol{\theta}$ posterior distributions. We see that there is no significant impact due to reliability uncertainty apparent in 1000 inference runs. This is in spite of the broad distributions of the low-reliability planet candidates shown in Figure 29, which shows the false alarm reliability values $R_{\mathrm{R}A}$ resulting from the $\boldsymbol{\theta}$ posterior distributions. We believe this lack of impact on occurrence rates is due to the very small uncertainties for high-reliability targets (see Figure 23) combined with the facts that the low-reliability targets occur with less frequency in each inference and that most distributions in Figure 29 are close to symmetric, so are well-represented by their medians.

6.4 Variations

In this section we explore the impact of changing some of the inputs and assumptions in the baseline occurrence rate computed in §6.2. Our motivation is to understand the dependencies of the occurrence rate on these inputs and assumptions. In all cases except §6.4.2 the same models were found to be the best fit to vetting completeness, false alarm effectiveness and observed false alarm rate as in the baseline case, though the parameters of these models had different values for the different variations.

We present the resulting variation in occurrence rates in Table 4, which includes results from Table 1 for comparison. This comparison is shown graphically for $F_{1}$ and $\zeta_{\oplus}$ at the end of this section in Figure 32.

6.4.1 Using the Q1-Q17 DR25 Stellar Properties

Our baseline occurrence rates are substantially lower than several occurrence rates based on pre-Gaia stellar properties. In this section we repeat our analysis, replacing the Gaia-based catalog of Berger et al. (2020) with the pre-Gaia Q1-Q17 DR25 stellar properties from the NASA Exoplanet Archive4. We perform the same cuts as described in §3.1, with the exception that there is no cut on binary or evolved flags (these do not exist in the Q1-Q17 DR25 stellar properties) and we remove all stars with radius $>1.35R_{\odot}$. The final catalog contains 75,541 GK stars.

We perform the same analysis as in the baseline case, starting from computing the vetting completeness for this stellar catalog, computing the summed completeness function $\eta$, the reliability and occurrence rates specified in §6.2. Figure 30 shows the resulting planet population, summed completeness and reliability. Comparison with Figure 24 shows that this catalog has more planet candidates in our period-radius range than when using Berger et al. (2020). This results in the higher occurrence rates shown in the “DR25” case in Table 4.

The choice of catalog has a stronger impact on $\zeta_{\oplus}$ than on $F_{1}$: When not correcting for reliability, $F_{1}$ based on the DR25 stellar properties are about 15% higher than our baseline using Berger et al. (2020), while the DR25-based $\zeta_{\oplus}$ is about 60% higher. When correcting for reliability, $F_{1}$ based on the DR25 stellar properties are about 20% higher, while the DR25-based $\zeta_{\oplus}$ is 80% higher. Computing the SAG13 definition of $\eta_{\oplus}$10using the DR25 stellar properties without correcting for reliability yields $\eta_{\oplus}=0.499^{+0.245}_{-0.164}$, while correcting for reliability gives $\eta_{\oplus}=0.223^{+0.136}_{-0.087}$.

The Berger et al. (2020) stellar catalog used in our baseline occurrence rates differs from the DR25 stellar catalog used in this section in both the values of the stellar properties themselves and in the cuts used to define the stellar parent population. In Appendix C we study the relative impact of the difference in stellar properties vs. the impact of the different population cuts on the difference in occurrence rates. We find that the difference in occurrence rates is primarily due to the different stellar properties, primarily stellar radius and effective temperature (leading to different GK selections), and that the differing population cuts have a minor impact.

6.4.2 Baseline with a Score Cut of 0.9

The Robovetter outputs a score for each TCE, indicating the confidence with which the Robovetter vetted that TCE (Thompson et al., 2018). This score is not equivalent to reliability: for example the Robovetter confidently vetted several TCEs in the inverted/scrambled data incorrectly as PC with scores as high as 0.923. But score is roughly correlated with reliability, and Thompson et al. (2018) suggests computing high-reliability occurrence rates by considering only planet candidates with Robovetter score above some threshold. This will result in a smaller planet candidate population with lower completeness, but the resulting larger completeness correction will, in principle, correct the occurrence rate.

In this variation we impose an aggressive score cut, rejecting any planet candidate with score $<0.9$. We use the Berger et al. (2020) catalog, and compute the completeness and reliability as in the baseline case, treating any TCE with score $<0.9$ as a false positive/alarm. Mulders et al. (2018) uses this score cut in their analysis, but their analysis is on a very different period-radius range so is not directly comparable to our results.

The result is a smaller, higher reliability planet candidate population, as shown in Figure 31, with noticeably lower completeness (compare the contours in Figure 24). In this case the false alarm vetting efficiency was best fit with the constant $=0.999$, resulting in a false alarm reliability very close to 1 for the entire period-radius range. The few planet candidates with lower reliability in Figure 31 are due to their astrophysical false positive probability, which results in $F_{1}$ and $\zeta_{\oplus}$ being slightly suppressed as shown in the “Score $>0.9$” case in Table 4. This is an illustration of the fact that score cuts cannot be relied on to provide a population that is high reliability with respect to astrophysical false positives.

The agreement in $\zeta_{\oplus}$ when using this score cut and the baseline given in Table 1 is remarkable given the lack of planet candidates smaller than $2R_{\oplus}$ and orbital period $>220$ days shown in Figure 31 (compare Figure 24). We interpret this agreement as an indication that the baseline $\zeta_{\oplus}$ is dominated by extrapolation because, in the baseline population, long-period, small planets have low reliability, as discussed in §6.2. Because the baseline results and those using requiring score $>0.9$ are essentially unconstrained extrapolations from radius $>2R_{\oplus}$ and orbital period $<220$ days to smaller planets and longer periods, we believe it is premature to conclude that using this score cut provides accurate occurrence rates for radius $<2R_{\oplus}$ and orbital period $>220$ days. In §7 we propose a strategy to explore this question.

6.4.3 Baseline Without Vetting Completeness

This variation measures the impact of not including vetting completeness. This will result in a smaller completeness correction where vetting completeness is low, so we expect somewhat lower long-period, small planet occurrence rates. The “No Vetting Efficiency” case in Table 4 shows a small suppression in $\Gamma_{\oplus}$, $F_{1}$ and $\zeta_{\oplus}$ when not accounting for vetting completeness.

6.4.4 Baseline Without MES Smearing

This variation measures the impact of not smearing the MES in the calculation of detection completeness, described in §4.1. The “No MES Smear” case in Table 4 indicates an increase in small-planet, long-period occurrence rates measured by increases in $\Gamma_{\oplus}$ and $\zeta_{\oplus}$, but not smearing the MES has essentially no impact on $F_{1}$.

7 Discussion

In this paper we show that a proper characterization of vetting completeness and reliability is important, particularly near the detection limit. In particular, in §6.2 we find that characterizing Kepler reliability and completeness can impact occurrence rates by more than a factor of two near Kepler’s detection limit (see Table 1). We introduce a new approach to characterizing vetting completeness and reliability for the Kepler DR25 planet candidate population. This approach casts the problem as one of binomial probabilities via parameterized rate functions fitted to the DR25 injection, inverted and scrambled data. We develop parameterized models of completeness (described in §4), false alarm effectiveness (§5.1.1) and the observed false alarm rate (§5.1.2). The particular parametric models we choose are selected via the Akieke information criterion, which chooses the parametric model that maximizes the likelihood corrected for the number of model parameters (see Appendix A). We do not claim that our parametric models are the best or in any sense “true”, just that they are the best of the parametric models we considered, described in Table 5. But our best models do a good job of accounting for the data, and are robust against choices such as grid resolution.

We caution, however, that vetting completeness and false alarm reliability as defined in this paper are properties of the specific Robovetter metrics and vetting thresholds behind the DR25 planet candidate catalog, as well as our analysis method, rather than properties of the detections themselves. For example, a different choice of Robovetter metrics may increase completeness while decreasing reliability or vice versa. While a low reliability for a transit detection from the analysis in this paper is reason to be cautious about asserting that detection is due to a true planet, further analysis of Kepler data can potentially result in higher confidence that a transit signal is due to a true planet. For example, as described in §2, at least one major source of false alarms, rolling bands, is highly dependent on focal plane position (Van Cleve & Caldwell, 2009; Cleve et al., 2009). Though some of the DR25 vetting metrics, such as skye (Thompson et al., 2018) are focal-plane dependent, the reliability analysis in this paper largely ignores focal plane dependence by averaging over the focal plane, and potentially underestimates the reliability of a detection in a focal plane position known to have a low occurrence of, for example, rolling bands. Pixel-level analysis of transit events beyond that used in DR25 may be useful in distinguishing false alarms due to statistical fluctuations and cosmic ray events. These observations can potentially be implemented as new Robovetter metrics, which could result in a higher-reliability, more complete planet candidate catalog.

The four years of Kepler’s observation of its primary field provides a data set unlikely to be excelled in the near future. Full exploitation of this data for understanding exoplanet populations is only partially complete. This paper is an attempt to fill in a significant step in that exploitation. We deliberately chose to limit the innovations in this paper to the characterization of and correction for completeness and reliability, and the use of the uniform Gaia-based stellar properties catalog of Berger et al. (2020). We show the impact of these innovations by computing occurrence rates using standard methods from Burke et al. (2015) in order to facilitate comparison with previous occurrence rates based on similar methods. The following discussion critically examines the assumptions underlying these occurrence rates, revealing weaknesses in both the DR25 catalog and the occurrence rate calculation method, and outlines some of the directions that we believe will prove fruitful in addressing these weaknesses.

7.1 Assumptions Underlying the Baseline Occurrence Rate

We illustrate the impact of reliability by computing a variety of occurrence rates near the Kepler detection limit (see §6). We chose our specific occurrence rate method, Bayesian inference using a dual power law population model in period and radius, because it is standard and well-understood. We believe that our occurrence rates provide high-confidence insight into what the DR25 planet candidate catalog tells us about the exoplanet population for the period and radius range of $50\leq\mathrm{period}\leq 400$ days and $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$ under the following assumptions:

• •

The parent stellar population is statistically well-described by Berger et al. (2020).

• •

Detection and vetting completeness in the injected data, along with the analysis described in §4, capture the statistical behavior of detection and vetting completeness in the observed data.

• •

The false alarms in the inverted and scrambled data capture the statistical behavior of the false alarms in the observed data.

• •

The astrophysical false positives probability statistically captures the probability of astrophysical false positives in the planet candidates.

• •

Exoplanets are distributed according to a Poisson point process.

• •

Dependence of planet occurrence on period and radius is modeled by a product of power laws in period and radius.

We discuss each of these assumptions in turn.

The parent stellar population: As stated in §3.1, we choose the Berger et al. (2020) stellar properties because they are informed by Gaia radii and uniformly treat the stellar properties across the parent population. While detailed observations of individual stars may provide more accurate stellar properties for individual stars, our method requires statistically uniform analysis. This is provided by the isochrone fitting approach of Berger et al. (2020). Therefore we believe that this stellar catalog is very well suited to our analysis. We showed in §6.4.1, however, that the occurrence rate depends critically on the stellar properties in the parent population. Inaccuracies, in particular biases in stellar radius estimates, can have a strong impact on occurrence rates.

Detection and vetting completeness: The injected data and analysis from §4 makes many assumptions. The detection completeness analysis makes several empirical approximations (described in Burke & Catanzarite (2017)) that may not apply well to individual planet candidate host stars. As described in §3, some of our stellar and planet candidate population approaches the restrictions stated in Burke & Catanzarite (2017) with respect to transit duration. Because our occurrence rates include regions with very few planet candidates, it is possible that the detection completeness for long-period planets is not as well modeled as that for short-period planets. Regarding vetting completeness, we are assuming that the Robovetter vets the injected detections with the same statistical accuracy as real transiting planets in the observed data. While we have confidence in this assumption, it is possible that some true planet transit signals have properties not captured by injection which confound the Robovetter, such as asymmetric transit shapes due to non-zero eccentricity, transit timing variations, and out-of-transit flux variations.

False alarm characterization: As we discussed in §5, inverted and scrambled data is believed to statistically model three identified classes of false alarms: rolling bands, statistical fluctuations combined with cosmic ray-induced pixel events that conspire to imitate long-period small transiting planets and stellar variability. The evidence for this belief is that the distribution of detected TCEs in the inverted and scrambled data closely matches the clearly anomalous distribution of detections in the observed data centered on the Kepler orbital period (see Thompson et al. (2018)), and tuning the Robovetter to eliminate this distribution from the PC population in the inverted and scrambled data also eliminates the anomalous distribution in the observed data. While using the inverted and scrambled data to model the false alarm population is clearly effective, there are likely other types of false alarms not represented by the inverted/scrambled data, though these seem to be a minor component compared to those represented by inverted/scrambled data. Such unmeasured false alarms would cause an overestimate of false alarm reliability.

Astrophysical false positive characterization: The false positive probabilities (Morton et al., 2016) are computed making strong assumptions about the lack of evidence for stellar multiplicity associated with the transit host star. While these false positive probabilities model stellar multiplicity as candidate hypotheses, the prior used in this model strongly assumes a lack of evidence for stellar multiplicity. As described in Hirsch et al. (2017) and Ciardi et al. (2015), there is evidence that a non-trivial fraction, possibly 20%, of Kepler target stars have unknown stellar companions. Such companions could cause an overestimate of the reliability of a subset of the PC population.

Poisson Likelihood: The use of the Poisson likelihood (Equation (11)) for the distribution of exoplanets is a standard choice, but may not be correct. For example, the assumption that the probability of different planets on the same star are independent of one another (an assumption behind Equation (11)) is almost certainly not correct, as indicated by existence of many packed exoplanet systems. There is also evidence that the detection of one planet on a star can prevent the detection of other planets on the same star (Zink et al., 2019). Likelihood-free methods, such as approximate Bayesian computation as applied to occurrence rates in Hsu et al. (2018) or the population sampling method used in Zink & Hansen (2019) may yield more accurate occurrence rates.

The power law population model: Evidence is mounting against the use of a simple product of power laws in period and radius when modeling exoplanet population statistics. This is already apparent in the top-left panel of Figure 27, where the power law is a poor fit to the observed planet population as a function of radius. This is likely due to the Fulton gap (Fulton et al., 2017), though the orbital periods in our analysis are somewhat longer than in Fulton’s analysis. Further, Petigura et al. (2018) presents evidence that host star metallicity is an important parameter in exoplanet population statistics. As pointed out by Hsu et al. (2018), model mis-specification is unlikely to lead to accurate results. Several authors have avoided the use of parameterized models in occurrence rate computations (for examples, see Foreman-Mackey et al. (2014); Hsu et al. (2018); Howard et al. (2012)), which is likely to lead to more accurate occurrence rates.

We believe that whatever method of statistical analysis is applied to the planet candidate catalog, characterizing and correcting for vetting completeness and reliability is critical. In the long-period, small-planet regime we have shown that reliability can reduce occurrence rates by a factor of two. We expect that this will be the case regardless of the statistical method and model, because reliability is a property of the planet candidate catalog. The effect of vetting completeness is less dramatic in the DR25 planet candidate catalog (as opposed to detection completeness, which is very important), but should not be neglected. Other planet catalogs may increase vetting reliability at the expense of vetting completeness, in which case vetting completeness can be more significant.

For habitability studies, the common practice of grouping together a wide class of stars and computing occurrence rates as functions of period and radius is potentially misleading. For example, the large range of stellar luminosities in our GK population shown in Figure 3 means that not all stars share the same habitable zones expressed as orbital periods. But grouping such a wide class of stars is necessary to provide the required statistics due to the small number of long-period, small planet detections and the sparseness of false alarms described in §5.1.1. Occurrence rates computed as functions of insolation flux and planet radius for the same class of stars would provide the needed statistics, and are likely more informative for habitable exoplanet population studies. Recent improvements in stellar characterization of the parent stellar sample, represented by Berger et al. (2020), make insolation-flux based occurrence rates a viable alternative.

7.2 Improving the Planet Candidate Catalog

The discussion in §7.1 outlines the assumptions behind extracting our occurrence rate from the DR25 planet candidate catalog, and how those assumptions may fail. The DR25 catalog itself can likely be improved upon, particularly in the long-period small planet regime. There is evidence that several long-period, small radius detections were incorrectly classified as false positives: the Kepler False Positive Working Group (FPWG) (Bryson et al., 2015) has identified several TCEs vetted as false positive in DR25 that are viable planet candidates, identified with fpwg_disp_status = POSSIBLE PLANET in the Kepler certified false positive table at the NASA Exoplanet archive2. This is expected because the DR25 vetting process deliberately balanced statistical uniformity and accuracy for individual objects, which was required for the study in this paper but compromised accurate vetting for some objects. In principle, the lowered completeness resulting from mis-classifying true planets as false positives is corrected by characterizing vetting completeness. But in the long-period, small-planet regime there are very few, low-reliability detections and very large completeness corrections (see Figure 24), which is vulnerable to large errors due to small statistics.

Accurate characterization of completeness and reliability as developed in this paper opens an intriguing approach to addressing the problem of few detections at long period and small radius: planet candidate catalogs that have lower reliability and higher completeness. This would mitigate the small-statistics problem by providing more detections with a smaller completeness correction resulting in better statistical constraints on the extrapolations discussed in §6.4.2. We recommend an exploration of Robovetter thresholds that increase the number of detections, lowering reliability and increasing completeness. The methods to measure reliability described in this paper are a crucial step towards being able to extract more accurate occurrence rates from such a catalog.

We believe that the reliability of the planet candidate catalog can be improved by the development of metrics beyond those described in Thompson et al. (2018). We provide two examples that may prove fruitful.

• •

As described above, different regions of the Kepler focal plane have different different false alarm characteristics, which can be leveraged to more accurately evaluate the likelihood that a transit signal is due to a false alarm.

• •

Pixel-level analysis can be developed beyond the DR25 vetting metrics, based on the expectation that false alarms are likely to be significantly different from star-like transit signals at the pixel level, particularly in difference images (Bryson et al., 2013).

We expect that such improved vetting metrics will address the small statistics problem by increasing the reliability of planet candidates near the detection limit, so they have stronger statistical weight in occurrence rate calculations.

Followup observation can potentially play a role in validating the reliability characterization developed in this paper. Ground- or space-based observations of a significant number of DR25 PCs, confirming them as planets or determining them to be false positives, could provide a ground truth of the number of PCs that are true planets. This ground truth can be used to independently compute the reliability of the DR25 PC population. We caution against using such followup observations to modify the planet candidate catalog, however, as that is likely to violate the uniformity assumptions behind the completeness correction.

8 Conclusion

This paper presents a new, probabilistic approach to statistically characterizing the vetting completeness and reliability Kepler DR25 exoplanet catalog. Using a standard occurrence rate calculation, we demonstrated that correcting for reliability can have a significant impact on occurrence rates, particularly near the Kepler detection limit at orbital periods longer than 200 days and planet radius $<$1.5 R*⊕*. We also showed that the choice of stellar properties for the searched stellar sample has a significant effect on occurrence rates. The results in this paper were made possible by the uniform detection and vetting methods behind DR25 that lend themselves to statistical characterization. We believe that the results presented in this paper are directly applicable to other exoplanet surveys such as K2, TESS, and PLATO so long as they create their catalogs in a similarly uniform way and expend the effort to create test data sets that measure completeness and false positives.

We thank NASA, Kepler management, and the Exoplanet Exploration Office for continued support of and encouragement for the analysis of Kepler data. We thank Bill Borucki and the Kepler team for the excellent data that makes these studies possible. We thank Daniel Foreman-Mackey for his Python Jupyter notebooks, which kicked off and inspired the work presented in this paper. T.A.B and D.H. acknowledge support by the National Science Foundation (AST-1717000).

Appendix A Vetting Completeness and Reliability Model Selection

We investigated a variety of models for vetting completeness (§4.2), false alarm efficiency (§5.1.1) and observed false alarm rate (§5.1.2). We selected the models based largely on the lowest Akiake Information Criterion $\mathrm{AIC}=2d-2\ln(L)$, where $d$ is the number of degrees of freedom in the model and $L$ is the likelihood of the model. When comparing two models with likelihoods $L_{1}$ and $L_{2}$, the relatively likelihood of model 1 relative to model 2 is given by $\exp((\mathrm{AIC}_{2}-\mathrm{AIC}_{1})/2$. The AIC criterion is not always successful, however, particularly when the model contains parameters that do not converge. In addition, the lowest AIC sometimes results in models that are obviously not physical. In such cases we made judgement calls when making the selections, as described below.

Details of each model’s analysis is found on the GitHub website3 in the directory GKbaseline/htmlArchive.

The models we considered are defined in Table 5. All models have as input orbital period $p$, MES (either expected for vetting completeness, or observed for reliability) $m$, orbital period range $[p_{\min},p_{\max}]$, MES range $[m_{\min},m_{\max}]$ and a parameter vector $\boldsymbol{\theta}$. $\boldsymbol{\theta}$ has different elements for different models. The function evaluations start with scaling the period and MES input to $(x,y)\in[0,1]\times[0,1]$ in order to facilitate rotations. This scaling is done for consistency even if there are no rotations in the model. In some models, rotations are applied using angles $\phi_{x}$ and $\phi_{y}$ from $\boldsymbol{\theta}$. If there is only one rotation angle $\phi$ in $\boldsymbol{\theta}$, then $\phi_{x}=\phi_{y}=\phi$.

[TABLE]

Many of the functions are constructed from the logistic function $Y\left(x,x_{0},k,\nu\right)$ defined in equation (4). We also use the broken power law and non-normalized two-dimensional Gaussian

[TABLE]

A.1 Vetting Completeness Model Selection

The functions considered for fitting the observed vetting completeness rate in §4.2 are given in Table 6, along with AIC values and relative likelihoods. We chose logisticX0xRotatedLogisticY0 because it had the highest relative likelihood, best convergence behavior, and appears to be a good fit to the data.

A.2 False Alarm Effectiveness Model Selection

The functions considered for fitting the observed false alarm effectiveness rate in §5.1.1 are given in Table 7, along with AIC values and relative likelihoods. We chose rotatedLogisticX0 because it has a high relative likelihood, had the fewest parameters, and gave the most reasonable convergence results compared with other high-likelihood models.

A.3 Observed False Alarm Rate Model Selection

The functions considered for fitting the observed false alarm rate in §5.1.2 are given in Table 8, along with AIC values and relative likelihoods. We chose rotatedLogisticX0 because it gave the most reasonable results compared with other high-relative-likelihood models. We rejected rotatedLogisticX02 because one of its parameters did not converge.

Appendix B Derivation of the Likelihood from the Poisson Probability

We briefly summarize Bayesian inference using a Poisson likelihood. We will work in the period-radius parameter space.

If our planet population is described by a point process with a period and radius-dependent rate $\lambda(p,r)$, then the probability that $n_{i}$ planets occur around an individual star in some region $B_{i}$ (say a grid cell) of period-radius space is

[TABLE]

where

[TABLE]

We now cover our entire period-radius range $D$ with a sufficiently fine regular grid with spacing $\Delta p$ and $\Delta r$ so that each grid cell $i$ centered at period and radius $(p_{i},r_{i})$ contains at most one planet. Then in cell $i$

[TABLE]

We now ask: what is the probability of a specific number $n_{i}$ of planets in each cell $i$? We assume that the probability of a planet in different cells are independent, so

[TABLE]

because the $B_{i}$ cover $D$ and are disjoint. Here $K$ is the number of grid cells and ${K_{1}}$ is the number of grid cells that contain a planet = the number of planet candidates. So the grid has disappeared, and we only need to evaluate $\lambda(p,r)$ at the planet locations $(p_{i},r_{i})$ and integrate the rate function $\lambda$ over the entire domain.

We do not observe all the planets, however. We account for incompleteness, including geometric transit probability, by replacing $\lambda(p,r)$ with $\eta_{s}(p,r)\lambda(p,r)$ in equation (B1), where $\eta_{s}(p,r)$ is the completeness function for this star $s$ measured in §4. The result is the probability

[TABLE]

We now consider the probability of detecting planets around a set of $N_{*}$ stars. Assuming that the planet detections on different stars are independent of each other, then the joint probability of a specific set of detections specified by the set $\{n_{i},i=1,\ldots,N_{*}\}$ in cell $i$ on on all stars indexed by $s$ is given by

[TABLE]

where $V=\left(\Delta p\Delta r\right)^{(K_{1}N_{*})}$ and $\eta(p,r)=\sum_{s=1}^{N_{*}}\eta_{s}(p,r)$ is the sum of the completeness functions over all stars.

We now let the rate function $\lambda(p,r,\boldsymbol{\theta})$ depend on a parameter vector $\boldsymbol{\theta}$, and consider the problem of finding the $\boldsymbol{\theta}$ that maximizes the likelihood

[TABLE]

Because we are maximizing with respect to $\boldsymbol{\theta}$, we can ignore all terms that do not depend on $\boldsymbol{\theta}$. Therefore maximizing equation (LABEL:equation:poisson4) is equivalent to maximizing

[TABLE]

Appendix C Comparison of Catalog Cuts

In §6.4.1 we found that using the DR25 stellar properties catalog results in larger occurrence ratse than the baseline of §6.2, which uses the stellar properties from Berger et al. (2020). This difference can result from both the difference in the stellar properties themselves, and the fact that different stellar population cuts were used. In particular, as described in §3.1, the Berger et al. (2020) catalog contains only stars with good Gaia noise characteristics, and we impose further Gaia fit quality requirements by removing stars with qualityFlag = highRUWE.

In this appendix we explore the relative impact of the difference in stellar properties between the two catalogs compared with the impact of the different cuts. We consider the Berger et al. (2020) catalog with and without the cuts specific to this catalog. We then consider the same stars as Berger et al. (2020), but using the DR25 stellar properties and cuts. Finally we consider the DR25 stellar catalog and its restriction to those stars contained in the supplemental catalog of Mathur & Huber (2016). In all cases all steps of the occurrence rate computation are recomputed, including detection/vetting completeness and reliability.

We compute the occurrence rates $F_{1}$ and $\zeta_{\oplus}$, defined in §6. For two cases using the stellar properties of Berger et al. (2020) which differ in the population cuts:

• •

Case 1: the baseline of §6.2, starting with the Berger et al. (2020) catalog, with all the cuts described in §3.1, and planet radii corrected for Gaia stellar radii as described in §3.2. Starts with 186,548 stars and ends up with 58,974 GK stars after cuts.

• •

Case 2: Same as case 1, starting with the Berger et al. (2020) catalog, except without the highRUWE, Bin or Evol cuts described in §3.1, replacing these cuts with the cut on stellar radius removing stars with $R_{*}>1.35R_{\odot}$ described in §6.4.1. Starts with 186,548 stars, and ends up with 66,956 GK stars after cuts.

We examine three cases using the DR25 stellar properties, which differ in the population cuts:

• •

Case 3: The same cuts as case 2, starting with the Berger et al. (2020) catalog, except using DR25 stellar properties and original DR25 planet radii. Starts with 186,548 stars and ends up with 71,168 GK stars after cuts.

• •

Case 4: The DR25 stellar catalog as described in section 4.1 using original DR25 planet radii. Starts with 200,038 stars and ends up with 75,541 GK stars after cuts.

• •

Case 5: The DR25 stellar catalog as in case 4, restricted to those stars in Mathur & Huber (2016) and using original DR25 planet radii. Starts with 197,096 stars and ends up with 74,989 GK stars after cuts.

Cases 1 through 3 start with the same stars, and differ in the cuts and the use of Gaia-based vs. DR25 stellar properties.

The occurrence rates $F_{1}$ and $\zeta_{\oplus}$ for the various cases are given in Table 9 and shown in Figure 33. We see that using the same stellar properties gives similar occurrence rates, with a noticeable difference in occurrence rates computed using different stellar properties. Differing cuts using the same stellar properties apparently has a much smaller impact. We therefore conclude that stellar properties (including differences in GK classification due to differences in effective temperature) is the dominant cause of the different occurrence rates, and the population cuts play a minor role. In all cases correcting for reliability has a significant impact on $\zeta_{\oplus}$.

We can get some insight into the change in occurrence rates by examining the impact of stellar properties on the PC population in the parameter space $50\leq\mathrm{period}\leq 400$ days and $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$. We examine the difference between case 2 and case 3 because these cases start with the same parent stellar population and apply the same cuts. In case 2 these cuts are applied using the Berger et al. (2020) stellar properties, while in case 3 these cuts are applied using the DR25 stellar properties. In case 2 there are 107 PCs in the period and radius range out of 67,306 stars in the parent GK population, which yields 0.666 per star after dividing by the average completeness of 0.00239. In case 3 there are 116 PCs out of 71,168 stars in the parent GK population, which yields 0.64 planets per star after dividing by the average completeness of 0.00255. These simple estimates are consistent with the values of $F_{0}$ found in Table 9. To understand the changes in $F_{1}$ and $\zeta_{\oplus}$, we need to examine in more detail how the stellar properties effect the PC population.

There are three ways that a difference in stellar properties can change the occurrence rate:

• •

A change in stellar radius causes a change in planet radius, which will impact the period-radius dependence of the population rate function and may cause the planet to move into or out of the $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$ range.

• •

A change in stellar radius causes the star to be added to or removed from the parent population depending on whether the star becomes smaller or larger than the $1.35R_{\odot}$ cut.

• •

A change in stellar effective temperature causes the star to be added to or removed from the parent population because it is reclassified as GK or not GK.

Figure 34 shows the change in planet radius when changing from case 3 (DR25 stellar properties) to case 2, for those planet candidates that are common to both case 2 and case 3. We see that for periods between 50 and 200 days planet candidates both increase and decrease in size. For periods greater than 200 days, however, there is clear bias towards larger sizes. This effect is quantified by computing the average relative change in size in three period bins. The shortest period bin shows a near-zero average change in size, while the longest-period bin shows an average increase in size of about 8%, which is about a $2\sigma$ change.

Figure 35 shows planet candidates that either exited or entered the radius range considered in our occurrence rate in the change from case 3 to case 2. As in Figure 34, we see that at low period several planet candidates entered our planet radius domain of $0.75\leq\mathrm{radius}\leq 2.5$ $R_{\oplus}$ while other planet candidates left that domain. But for longer period, particularly $>250$ days, planets left our domain by becoming too large while no planets entered our domain. Thus there is a loss of small exoplanets due to their being larger using Berger et al. (2020) stellar properties. In addition, several stars exited or entered our parent population through reclassification due to change in effective temperature.

These figures suggest that the change in planet size when using the Berger et al. (2020) stellar properties (case 2) is a significant contributing factor in the reduced occurrence rates. But it would not be correct to conclude that this change in planet size is the “cause” of the lower occurrence rate: changes in the parent stellar population also impact occurrence rates through changes in detection completeness and changes in population due to stellar reclassification and which stars pass the stellar size cut. It is only through computing the full occurrence rate that we can measure the impact of the stellar properties.