THE OCCURRENCE OF POTENTIALLY HABITABLE PLANETS ORBITING M DWARFS ESTIMATED FROM THE FULL KEPLER DATASET AND AN EMPIRICAL MEASUREMENT OF THE DETECTION SENSITIVITY

In this paper, we focus on the population of planets orbiting small stars. Such planets are important for constraining the galactic census of exoplanets because the majority of stars in the galaxy are low-mass stars (Henry et al. 2006; Winters et al. 2015). In addition to understanding the overall occurrence rate of planets orbiting low-mass stars, we would like to know how planet occurrence depends on factors such as planet radius, orbital period, stellar insolation, and host star properties. Furthermore, small stars afford the best near-future opportunities for detailed characterization studies of small planets and their atmospheres (Charbonneau & Deming 2007). In order to prepare for these observations, we would like to know the likely distance to the closest such targets.

Computing the occurrence rate of small planets around small stars is complicated by the fact that the parameters of low-mass stars are more difficult to measure than the parameters of Sun-like stars. The main emphasis of the Kepler mission was the detection of planets around Sun-like stars, so the assumptions made in the construction of the Kepler Input Catalog (KIC) were tailored to be appropriate for Sun-like stars. Accordingly, Brown et al. (2011) cautioned against relying on KIC classifications for stars cooler than 3750 K.

Several previous studies have attempted to improve the KIC parameters for the coolest Kepler target stars. In our previous paper (Dressing & Charbonneau 2013), we used the photometry provided in the KIC to reclassify all stars cooler than 4000 K using models (Dotter et al. 2008; Feiden et al. 2011) and assumptions more appropriate for low-mass stars. Gaidos (2013) conducted a similar analysis for the population of planet candidate host stars. Other authors further constrained the properties of particular low-mass stars and associated planet candidates by acquiring follow-up spectroscopic and high resolution imaging observations (Johnson et al. 2012; Muirhead et al. 2012a, 2012b; Ballard et al. 2013; Muirhead et al. 2013; Swift et al. 2013). Recognizing the importance of characterizing the full target sample as well as the planet host stars in order to constrain the planet occurrence rate, Mann et al. (2012) acquired spectra for a subset of stars without detected planets. They found that the majority of bright (Kp < 14) Kepler target stars are giant stars (several hundred of which were classified as dwarfs in the KIC) but that 93% of fainter stars are correctly classified as dwarfs.

In this paper, we combine the current best estimates of the properties of small Kepler target stars in order to estimate the frequency of small planets around small stars. Our analysis was preceded by several studies of the planet occurrence rate based on Kepler data and we adopt some techniques from the earlier studies. In particular, we draw upon the framework established by Howard et al. (2012), Fressin et al. (2013), Dressing & Charbonneau (2013), Petigura et al. (2013b), and Petigura et al. (2013a).

Working with the first three quarters of Kepler data, Howard et al. (2012) estimated the frequency of planets around main-sequence GK stars. They found that the occurrence rate of planets increased sharply with decreasing planet size and moderately with increasing orbital period. They also found evidence for a cutoff period below which the planet occurrence rate falls off more quickly with decreasing period. The position of the cutoff period appeared to move outward from near 2 days for larger planets (R_p > 8R_⊕) to roughly 7 days for 2− 4R_⊕ planets.

Youdin (2011) used the search completeness estimates from Howard et al. (2012) to model the planet occurrence rate around Sun-like stars by a joint power law in orbital period and planet radius. He found an occurrence rate of 0.19 planets per star with periods shorter than 50 days and radii larger than 2R_⊕. He also extrapolated outward to predict an occurrence rate of roughly three Earth-like planets per star with periods shorter than a year.

Due to the low false positive rate expected for the Kepler planet candidate sample (Morton & Johnson 2011), the Howard et al. (2012) and Youdin (2011) analyses assumed that all of the candidates were bona fide transiting planets. They further assumed that Kepler would have been able to detect all transiting planets with cumulative signal-to-noise ratio (SNRs) above a set threshold of 10σ. The first assumption biased their occurrence rate estimates toward higher values while the second assumption would have resulted in an underestimate if the actual search completeness were lower.

Fressin et al. (2013) conducted a follow-up study of the Kepler planet occurrence rate incorporating both contamination from false positives and a more sophisticated model of pipeline sensitivity. In particular, they used a hierarchical approach in which they first estimated the population of Jupiter-size planet candidates that might be astrophysical false positives. They then iteratively determined the occurrence rate of small planets by modeling the fraction of larger planet candidates that might be masquerading as smaller planet candidates in diluted transit events. Fressin et al. (2013) found a global false positive rate of 9.4 ± 0.9% and noted that considering false positives is particularly important when calculating the occurrence rate of giant planets (6−22R_⊕, FP rate = 17.7% ± 2.9) and Earth-size planets (0.8–1.25R_⊕, FP rate = 12.3% ± 3.0).

Fressin et al. (2013) also used their hierarchical model to estimate the completeness threshold of the Kepler pipeline. They found that a linear ramp model in which 0% of signals with SNRs <6 and 100% of signals with SNRs >16 were detected provided a better fit to the observed planet candidate population than an abrupt step function. Accounting for both a non-zero false positive rate and a ramp sensitivity model, Fressin et al. (2013) estimated that 14.9 ± 2.4% of FGK stars host an Earth-size planet (0.8–1.25R_⊕) with a period between 0.8 and 50 days. In a similar study of Kepler data, Dong & Zhu (2013) found that roughly 20% of main sequence stars with 5000 K < T_eff < 6500 K host 1–2R_⊕ planets in periods less than 50 days.

More recently, Petigura et al. (2013b) developed their own planet search pipeline in order to search for additional planet candidates around Kepler stars. Their TERRA pipeline uses a custom light curve detrending algorithm based on principal component analysis (Petigura & Marcy 2012). After searching for planets around 42,557 relatively quiet GK stars, Petigura et al. (2013a) found that 7.7 ± 1.3% of GK stars host small planets (1–2R_⊕) in periods between 25 and 50 days. They also extrapolated to predict that 22% ± 8% of GK stars host 1–2R_⊕ planets receiving between 1/4 and 4 times the insolation received by the Earth. Their calculation incorporated a 10% correction for false positives. As a benefit of writing their own pipeline, Petigura et al. (2013a) were able to explicitly measure the completeness of their planet sample by injecting and attempting to recover transiting planets.

In a follow-up study, Foreman-Mackey et al. (2014) used the reported search completeness and planet candidates from Petigura et al. (2013a) to rederive the planet occurrence rate using a hierarchical Bayesian model. The Foreman-Mackey et al. (2014) analysis differed from the Petigura et al. (2013a) analysis in two key aspects: (1) Foreman-Mackey et al. (2014) considered measurement errors in the stellar and transit parameters and (2) they did not assume that the planet occurrence rate was flat in log period, instead using a flexible Gaussian process to model the occurrence rate assuming a smooth functional form. As a result, Foreman-Mackey et al. (2014) found an occurrence rate of potentially habitable Earth-size planets three times lower than the Petigura et al. (2013a) estimate.

Silburt et al. (2015) considered a sample of 76,711 Kepler target stars with radii of 0.8–1.2R_⊙ and estimated the search completeness using the reported Combined Differential Photometric Precision (CDPP, Christiansen et al. 2012). They employed an iterative simulation to investigate the dependence of the planet occurrence rate on planet radius without subdividing the data into bins and accounted for errors in the planet radii. For planets with periods of 20–200 days, Silburt et al. (2015) reported that the occurrence rate is higher for planets with radii of 2–2.8R_⊕ than for smaller or larger planets. In total, they estimated that a typical Sun-like star hosts 0.46 ± 0.03 planets with periods of 20–200 days and radii of 1–4R_⊕. In agreement with Petigura et al. (2013a), Silburt et al. (2015) noted that the planet occurrence rate is flat in log period. Within a broad habitable zone (HZ) extending from 0.99–1.7 AU, Silburt et al. (2015) estimated an occurrence rate of ${0.064}_{-0.011}^{+0.034}$ small (1–2R_⊕) planets per Sun-like star.

Focusing specifically on Kepler's smallest target stars, we (Dressing & Charbonneau 2013) estimated an occurrence rate of ${0.90}_{-0.03}^{+0.04}$ planets per star for 0.5–4R_⊕ planets with periods shorter than 50 days. We based our previous analysis on Q1–Q6 Kepler planet candidate list and assumed that Kepler detected all planets with cumulative SNRs > 7.1σ. Using conservative HZ limits from Kasting et al. (1993), we estimated an occurrence rate of ${0.15}_{-0.06}^{+0.13}$ potentially habitable Earth-size (0.5–1.4R_⊕) planets per small star. Kopparapu (2013) then revised this estimate to ${0.48}_{-0.24}^{+0.12}$ planets per star using the broader updated HZ boundaries from Kopparapu et al. (2013b). His result agreed well with the estimate of ${0.46}_{-0.15}^{+0.18}$ potentially habitable 0.8–2R_⊕ planets per star from Gaidos (2013). Unlike Kopparapu (2013), Gaidos (2013) adopted HZ boundaries corresponding to the 50% cloud cover case from Selsis et al. (2007).

Morton & Swift (2014) adopted a slightly different technique to estimate the frequency of small planets around small stars. They assumed that the planet radius distribution is independent of orbital period and modeled each planet using a weighted kernel density estimator when computing the occurrence rate. They found that the occurrence rate estimates from Dressing & Charbonneau (2013), Kopparapu (2013), and Gaidos (2013) for planets smaller than 1.4R_⊕ should be increased by an additional incompleteness factor of 1.6 if the assumption made by Morton & Swift (2014) about the period-independence of the planet radius distribution is correct.

Like Silburt et al. (2015), Gaidos et al. (2014) used an iterative simulation to estimate the planet occurrence rate, but they elected to focus on stars cooler than 4200 K. For orbital periods of 1–180 days and radii of 0.5–6R_⊕, Gaidos et al. (2014) calculated a cumulative occurrence rate of 2.01 ± 0.36 planets per M dwarf. Gaidos et al. (2014) also remarked that the planet occurrence rate is highest for planets with radii of approximately 1R_⊕ and lower for larger and smaller planets.

The frequency of potentially habitable planets around small stars has also been estimated from radial velocity surveys. Based on six years of observations with the HARPS spectrograph, Bonfils et al. (2013) estimated an occurrence rate of ${0.41}_{-0.13}^{+0.54}$ potentially habitable planets per M dwarf. Their definition of "potentially habitable" encompassed planets with $1\leqslant M\mathrm{sin}i\leqslant 10$ within the "early Mars" and "recent Venus" boundaries of the HZ presented in Selsis et al. (2007). A subsequent study by Robertson et al. (2014) revealed that GJ 581d, one of the two planets upon which Bonfils et al. (2013) based their occurrence rate estimate, is likely a manifestation of stellar activity. Robertson et al. (2014) reported a revised occurrence rate of 0.33 potentially habitable planets per M dwarf. The updated RV-based estimate is more similar to the estimates based on Kepler data, but accurately determining an occurrence rate with only a single planet (GJ 667Cc, Anglada-Escudé et al. 2012; Bonfils et al. 2013; Delfosse et al. 2013) is challenging. Additionally, direct comparison of planet occurrence estimates from RV and transit surveys is complicated by the need to employ a compositional model to translate planet masses into radii.

In this paper, we implement the following improvements to refine our 2013 estimate of the frequency of small planets around small stars.

• 1.  
  We use the full Q0–Q17 Kepler data set.
• 2.  
  We utilize archival spectroscopic and photometric observations to refine the stellar sample.
• 3.  
  We explicitly measure the pipeline completeness.
• 4.  
  We inspect follow-up observations of planet host stars to properly account for transit depth dilution due to light from nearby stars.
• 5.  
  We apply a correction for false positives in the planet candidate sample.
• 6.  
  We incorporate a more sophisticated treatment of the HZ.

In Section 2 we describe the selection of our stellar sample, which includes some stars whose parameters have been characterized spectroscopically via follow-up observations. We explain our planet detection pipeline in Section 3 and our procedure for vetting candidates in Section 4. We present light curve fits for the accepted planet candidates in Section 5. In Section 6, we assess the completeness of our pipeline. We then estimate and discuss the planet occurrence rate in Section 7 before concluding in Section 8.

We selected our stellar sample by first downloading a table of all 4915 stars with T_eff < 4000 and log g > 3 from the Q1–16 Kepler Stellar Catalog on the NASA Exoplanet Archive¹ (Akeson et al. 2013). This catalog is described in Huber et al. (2014) and combines the best estimates available for each star from a variety of photometric, spectroscopic, and asteroseismic analyses. The properties for the stars in the downloaded sample were primarily determined from photometry (Brown et al. 2011; Dressing & Charbonneau 2013; Gaidos 2013; Huber et al. 2014), but 2% of the sample had spectroscopically derived parameters (Mann et al. 2012, 2013b; Muirhead et al. 2012a; Martín et al. 2013). For the majority of the stars in the sample (79%), the stellar parameters were drawn from our 2013 analysis (Dressing & Charbonneau 2013).

Some of the stars in the downloaded subset had light curves indicative of binary stars, variable stars, or enhanced spot activity. In addition, some of the stars were observed only for a small number of days. In order to accurately estimate the planet occurrence rate for small stars, we wanted to select the subset of stars with the highest search completeness. We therefore performed the following series of cuts on the sample.

First, we counted the number of timestamps for which each star had "good" data (i.e., not flagged). We rejected all 2101 stars with fewer than 48940 unflagged long cadence data points. Since Kepler obtained long cadence data using 29.4-minute integration times, this cut requires 1000 days of data. One of the main goals of this paper is to measure the occurrence rate of potentially habitable planets and we wanted to ensure that Kepler would have been able to observe multiple transits of planets within the HZs (see Section 7.3) of the stars in our final sample. For reference, the median orbital period at the outer edge of the Kopparapu et al. (2013b) HZ for the stars in our final sample is 131 days and the longest period at the outer HZ is 207 days.

Second, we removed 63 stars that McQuillan et al. (2013) categorized as likely giants based on their stochastic photometric variability. The affected stars have red colors (median J − H = 0.83) consistent with their revised classification as giants. Although two of the stars had revised classifications from Dressing & Charbonneau (2013), the remaining 61 had parameters from the KIC.

For reference, we checked whether any of our target stars were known eclipsing binaries by consulting the Kepler Eclipsing Binary Catalog. We examined both the published Version 2 (Slawson et al. 2011) and the online beta version of the Third Revision.² Eleven of the target stars were listed in both versions of the catalog, six stars were listed in Version 2 only and three stars were listed in Version 3 only.

Six of the twenty targets with matches in the Eclipsing Binary Catalogs were listed in the NASA Exoplanet Archive as false positive systems (KID 5820218 = KOI 1048, KID 6620003 = KOI 1225, KID 8823426 = KOI 1259, KID 9761199 = KOI 1459, KID 9772531 = KOI 950, and KID 10002261 = KOI 959). Two were listed as planet candidate host stars (KID 5384713 = KOI 3444 and KID 11853130 = KOI 3263). Confirmed giant planet KOI 254.01 (KID 5794240, Johnson et al. 2012) was also included as an EB match because of the very large transit depth. As evidenced by the presence of KOI 254.01 in the EB catalogs, the catalogs contain both actual EBs and likely planets. Accordingly, we did not remove the twenty targets with matches in the EB catalogs from our target sample.

We then detrended all of the light curves as described in detail in Section 3.1 using smoothing lengths of 500, 1000, and 2000 minutes. We constructed a histogram of the flux distributions for each of the detrended light curves and measured the χ² of a fit to a Gaussian flux distribution. The flux distribution of a well-behaved single star should be Gaussian after detrending, but the flux distribution of an eclipsing binary can appear bimodal. We therefore flagged for visual inspection all 353 stars for which reduced χ² > 3 for any of the detrended light curves. We also measured the standard deviations σ₅₀₀, σ₁₀₀₀, and σ₂₀₀₀ of the three detrended light curves for each star and took the ratios of the standard deviations of light curves detrended using different median filters. We flagged 42 stars for which any of the ratios σ₅₀₀/σ₁₀₀₀, σ₅₀₀/σ₂₀₀₀, or σ₁₀₀₀/σ₂₀₀₀ were below 0.8. This cut was designed to pick out light curves for which the detrending algorithm failed to remove longer timescale variability.

Finally, we visually inspected the detrended and PDC-MAP photometry (Smith et al. 2012; Stumpe et al. 2012) for the 395 flagged stars with >1000 days of data and large χ² or standard deviation ratios. We rejected 207 stars with highly variable detrended light curves or strong indications of classification as an eclipsing binary. We also rejected KID 12207013 because of the unusual light curve morphology displayed in certain quarters. We accepted the remaining 187 flagged stars, increasing our final selected sample to 2543 stars. We note that the 208 stars that were rejected during the visual inspection stage are unlikely to harbor detectable planet candidates exactly because their light curves are highly variable. Similarly, our injection tests would likely recover only a small fraction of any planets injected into their light curves. Accordingly, the exclusion of these stars from our stellar sample has negligible effect on our estimated rates of planet occurrence.

Our final selected sample of 2543 stars is compared to the initial downloaded sample in Figure 1. The temperature range of the sample extends from 2661 to 3999 K, with a median stellar effective temperature of 3746 K. The median stellar radius is 0.47R_⊙ and the stars have radii spanning from 0.10R_⊙ to 0.64R_⊙ . The metallicity range is [Fe/H] = −2.5 to [Fe/H] = 0.56, with a slightly sub-solar median metallicity of [Fe/H] = −0.1. However, most of the metallicity estimates were derived from photometry (Dressing & Charbonneau 2013) and are not well-constrained. The brightest star in the sample has a Kepler magnitude Kp = 10.07, but the median brightness is Kp = 15.5. The faintest star has Kp = 16.3. The sample contains 100 known planet (candidate) host stars with 83 planet candidates and 80 confirmed planets.

Zoom In Zoom Out Reset image size

The first step in our planet detection pipeline was to clean the light curves to prepare them for the transit search. Next, we searched each light curve sequentially for planets, allowing the code to detect multiple planets per star when warranted by the data. We then accepted or rejected putative detections using the vetting procedure described in Section 4 and checked that the automatically accepted transits were not ephemeris matches with other KOIs. Finally, we visually inspected all surviving candidates and reviewed all available follow-up analyses produced by the Kepler team and the community. We discuss each step of the planet detection pipeline in more detail in the following sections.

3.1. Preparing the Light Curves

We obtained all available long cadence data for each target via anonymous ftp from the MAST.³ We then excluded all data points flagged as low quality and detrended each quarter of data independently. We produced three detrended versions of each light curve using a running sigma-clipped mean filter with widths of 500, 1000, or 2000 minutes. When calculating the mean, we excluded all points more than 3σ away from the median value of the light curve within the filtered region. We then divided the flux data by the smoothed light curve to obtain a detrended, normalized light curve for that quarter. Figure 2 provides an illustration of our light curve detrending process.

Zoom In Zoom Out Reset image size

Next, we searched for data gaps and anomalies within the detrended light curves. We defined a data gap as ≥0.75 days of missing photometry. Data gaps are frequently accompanied by sharp increases or decreases in flux that can confound searches for planets. We removed these events by excising all data points within 1 day of the start of a data gap or less than 3 days after the end of a data gap.

3.2. Searching for Transiting Planets

Our transit detection pipeline identified 3111 putative transit events associated with 534 stars. Some of those signals might have been systematics or astrophysical false positives instead of bona fide transiting planet candidates. Accordingly, we performed a series of cuts to select the events consistent with transiting planets. First, we visually inspected the candidate transit events to identify signals that were not clearly associated with spacecraft systematics or stellar activity. Of the 3111 candidate events, 511 events survived initial visual inspection. The 2600 events rejected at the visual inspection stage displayed morphologies consistent with classification as spacecraft systematics or sinusoidal brightness variations indicative of starspots rather than transiting planets. The 511 signals surviving visual inspection were associated with 246 unique host stars.

Some of the candidate events were harmonics of signals detected at integer multiples of the true period. We then ranked the accepted signals for each host star in order of decreasing Δχ² as calculated during the detection phase and iteratively fit and excised the transits of each signal in order to ensure that the lower Δχ² signals were not simply harmonics of the strongest signals. We rejected all signals with resulting Δχ² below 5σ. This "sequential vetting" step reduced the number of candidate events to 323 possible transits for 246 unique stars.

Next, we conducted a second, more intensive round of visual examination for the remaining candidate events. We compared the shapes and depths of odd and even transits, checked for the appearance of secondary eclipses, considered the depth of the putative transit relative to other possible features at the same orbital period, and investigated whether the putative transit events were dominated by a small number of deep events. After visually vetting the candidates, we checked whether putative signals had been previously identified as false positives and excluded all such events. We rejected 180 signals (including 9 known false positives) and accepted 143 signals associated with 97 unique stars.

During the vetting stage, we noticed that the phase-folded light curve for KOI 2283.01 exhibited two transit-like events with markedly different depths when folded to the 17.402 days orbital period listed in the Q1–16 KOI catalog. We therefore rejected KOI 2283.01 as a blend containing an eclipsing binary. This interpretation is consistent with the large observed centroid source offset shift of 5.3σ.

We then performed a final search for additional planets in the systems in which a previously detected planet had survived the vetting process. We executed this search by phase-folding the detrended light curve on the orbital periods of all accepted signals and removing all data points within plus or minus one best-fit duration of transit center prior to re-running the search process detailed in Section 3.2. In all cases we used the light curve that had been detrended using a 2000 minutes filter. The motivation for repeating this search after the first round of vetting was that uncertainties in the initial periods, durations, and transit centers of the accepted signals might have limited the effectiveness of the clipping performed in the initial search.

Our second round search revealed 104 candidate signals for 33 stars. We vetted the signals using the same vetting pipeline as in the first round search and accepted 14 additional candidate transiting planets associated with 11 stars. Next, we excised the transits of the signals accepted in the second round and performed a third round of transit searches. No additional signals were accepted during the third round. The full sample of 156 planet candidates included 142 signals associated with 96 stars revealed in the first round of searching and 14 signals associated with 11 stars revealed in the second round.

We compared the periods, P, and epochs, t₀, of the accepted planet candidates to the catalog of known eclipsing binaries, periodic variable stars, and KOIs compiled by Coughlin et al. (2014). We excluded the host star from the match process in order to avoid matching a signal to itself. We did not find any corresponding signals within our specified match tolerances of

$$\begin{eqnarray}\left|{P}_{\mathrm{match}}-P\right| & \leqslant & \mathrm{min}(2\;\mathrm{hr},0.001\times P)\quad \mathrm{and}\\ \left|{t}_{0,\mathrm{match}}-{t}_{0}\right| & \leqslant & \mathrm{min}(4\;\mathrm{hr},0.001\times P)\end{eqnarray} \tag{ 1 }$$

4.1. New Planet Candidate

The majority of the accepted planet candidates corresponded to signals previously identified as KOIs. We found that 155 putative planet candidates had periods and epochs matching those of known planet candidates or confirmed planets (Borucki et al. 2010, 2011a, 2011b; Batalha et al. 2013; Burke et al. 2014). We accepted one new signal in a system with previously known KOIs.

The KOI 1681 (KID 5531953) system contains three known planet candidates with periods of 6.51, 1.99, and 3.53 days. Our pipeline detected a planet candidate in the system with a radius of 1R_⊕ and a period of 21.9 days, roughly 11 times the 1.99 days orbital period of KOI 1681.02. We describe additional planet properties in Table 1. As shown in Figure 3, the transit signal is still visible in the light curve after the transits of the other three planets have been removed, suggesting that this is a new planet candidate rather than an alias of KOI 1681.02. The signal was later identified in the Q1–17 DR 24 Kepler pipeline run⁵ and is now listed as planet candidate 1681.04.

Zoom In Zoom Out Reset image size

Table 1.  New Candidate Accepted by Our Pipeline

KID	KOI	P	t0	Rp	a	Δχ2
		(days)	(days)	(R⊕)	(AU)
5531953	1681	21.913843	17.036402	1.03	0.16	41.4

Download table as:  ASCIITypeset image

4.2. Accounting for Transit Depth Dilution

4.3. False Positive Correction

4.4. Known Planet Candidates Missed by Our Pipeline

Accurate planet radius estimates are a key ingredient in the planet occurrence calculation. We therefore refined the preliminary transit parameters found in Section 3.2 by conducting a Bayesian Markov Chain Monte Carlo (MCMC) analysis with a Metropolis–Hastings acceptance criterion (Metropolis et al. 1953). We varied the orbital period P, epoch of transit center t₀, planet-to-star radius ratio R_p/R_*, semimajor axis to stellar radius ratio a/R_* and impact parameter b. We assumed that all of the orbits were circular and fixed quadratic limb darkening parameters to the values predicted from the stellar temperatures and surface gravities (Claret & Bloemen 2011). We conducted some of the planet fits using short cadence Kepler data to better constrain the shape of transit during ingress and egress.

For each planet candidate, we ran N chains starting at initial positions set by perturbing the initial solution found during the detection and validation process by up to 5σ in each parameter. We manually adjusted the step sizes for each parameter such that the acceptance fractions were between 10%–30%. We ran each chain for at least 10⁴ steps before initiating periodic convergence tests by calculating the Gelman–Rubin potential scale reduction factor $\hat{R}$ for each parameter (Gelman et al. 2004). We terminated the chains when $\hat{R}\lt 1.05$ for all parameters and then accounted for "burn-in" by removing all steps taken prior to the point at which the likelihood first became higher than the median likelihood of the chain. After merging the chains, we adopted the median values of each parameter as the best-fit value and assigned errors encompassing the 68% of values nearest to the chosen best-fit value. We provide the best-fit parameters for each detected planet candidate in Table 9 and display them in Figure 4.

Zoom In Zoom Out Reset image size

Several of the planet candidates in our sample exhibit large transit timing variations and were poorly fit by the method described above. In those cases, we adopted the transit parameters found in previous studies using combined fits to the individual transit times and the planet properties. For the candidates displaying TTVs, we incorporated fits from Rowe et al. (2014) for KOIs 248.01, 248.02, 248.03, 248.04, 314.01, 314.02, 314.03, 448.01, 448.02, 886.01, 886.02, and 886.03 and from Ioannidis et al. (2014) for KOIs 676.01 and 676.02.

In addition, we adopted the light curve parameters from the 2015 January 2 version of the NASA Exoplanet Archive for 5 candidates: 961.01, 961.03, 1681.01, 2329.01, and 3263.01. We also adopted the fits from Rowe et al. (2014) for KOIs 254.01, 430.01, 430.02, 438.01, 775.02, 868.01, and 961.02, from Swift et al. (2015) for KOIs 1902.01 and 3444.02 and from Cartier et al. (2014) for KOIs 1422.01, 1422.02, 1422.03, 1422.04, 1422.05, and 2626.01. The KOIs with parameters from Cartier et al. (2014) were detected in multi-star systems in which the identities of the host stars are unknown. As described in Section 4.2, we accounted for all possible system configurations by using fractional planets distributed around each of the possible host stars.

As depicted in Figure 4, we found that 149 of the accepted planet candidates had revised radii 0.5 < R_P < 4 R_⊕ and orbital periods 0.5 < P < 200 days. Of the remaining candidates, one had a shorter period (KOI 961.02, P = 0.45 days), one was too small (KOI 5692.01), and six were too large (KOI 254.01, 868.01, 901.01, 902.01, 1176.01, and 3263.01). KOI 868.01 also has an orbital period longer than 200 days.

In order to accurately measure the planet occurrence rate based on the results of our planet search, we needed to know the completeness of our planet candidate list. We measured the completeness of our planet detection pipeline by injecting transiting planets into the PDC-MAP light curves, detrending them, and running the detrended light curves through our detection algorithm. We did not introduce the signals at the pixel level and we are therefore unable to comment on how the initial light curve extraction process affects transiting planets. We refer instead to Christiansen et al. (2013) for a discussion of pixel-level effects. They found that transits injected at the pixel level are usually recovered with high fidelity (final SNR = 96%−98% of the expected SNR).

The transit detection process as modeled in this paper consists of two distinct stages: (1) the putative event is identified as a peak in the BLS periodogram and (2) the signal is accepted because a transit model provides a 5σ improvement to a straight-line fit. We took advantage of the two-step nature of the search process when determining the search completeness for each star in our survey.

For each star, we generated a set of 2000 trial planets with orbital periods drawn from a log uniform distribution extending from 0.5 to 200 days and uniformly distributed epochs of transit, radii (0.5–4R_⊕), and impact parameters (0–1). We then constructed transit models (Mandel & Agol 2002) for the trial planets using the assigned planetary parameters and limb darkening parameters estimated from the coefficients in Claret & Bloemen (2011) based on the stellar temperatures and surface gravities. We resampled the transit models to the 29.4-minute long cadence integration time.

The first question we wished to address was whether the transits of the trial planets would be accepted as 5σ detections in the ideal scenario in which the light curve was perfectly detrended and the orbital parameters were determined exactly. Accordingly, we first multiplied the detrended light curves by the transit models (Mandel & Agol 2002) generated using the assigned trial planet parameters. We then checked whether the difference in the χ² between the best-fit transit model (which we guessed perfectly) and a no-transit model exceeded the 5σ detection threshold of Δχ² = 30.863. If the transit model was not preferred, then we recorded the trial planet as a non-detection. For a typical star, 8% of the trial planets were rejected at this stage.

For the trial planets that would be accepted in the ideal scenario, we then conducted a more realistic test by multiplying the transit model by the raw PDC-MAP photometry before detrending using the straight line fit method described in Section 3.2. We then re-checked whether the transit model was preferred at 5σ. Trial planets that were not accepted at this stage were also recorded as non-detections. Overall, 4% of the trial planets that were accepted in the ideal case of a perfectly detrended light curve were not accepted in this more realistic test.

Finally, we tested whether the remaining trial planets would have been identified as peaks in the BLS periodograms by running a full test for at least 25 trial planets for each star. We selected the trial planets for the full test by ranking the signals detected in the second round of testing in order of increasing Δχ², where Δχ² was the value at which they were preferred to a non-transiting model. We chose the first 25 trial planets in the ranked list for which a random number draw yielded a result greater than 0.5. In other words, we thinned the sample of trial planets by 50% and selected those closest to the expected sensitivity threshold. If fewer than ten trial planets were recovered, we conducted up to 25 additional runs (for a total of 50 simulations) until at least ten planets were recovered.

For the selected trial planets, we multiplied the corresponding transit model by the raw PDC-MAP photometry for the assigned host star and detrended the light curve using a 2000 minute median filter as explained in Section 3.1. We then fed the injected light curves into the detection pipeline described in Section 3. Although we considered the full range of planet periods (0.5–200 days), we halted the search process as soon as the injected signal was detected. If the signal was not detected, we terminated the search using the usual conditions discussed in Section 3.2.

6.1. Predicting Transit Detectability

In total, we ran 83699 complete BLS injection simulations for the 2543 stars in our sample. We also injected 604278 planets that had a Δχ² below our 5σ detection threshold. The remaining 4398023 injected planets had a Δχ² above the detection threshold but were not tested in the full BLS simulation. We predicted the detectability of these trial planets by finding the fraction of BLS trial planets recovered as a function of the Δχ² computed in the second round of transit model tests. We ranked the BLS trial planets by Δχ² and computed the recovery fraction for each consecutive group of 1000 planets. Next, we smoothed the resulting histogram and predicted the likelihood of detection for the 4398023 non-BLS runs using a cubic spline interpolation based on the smoothed histogram. We limited the maximum detection likelihood to 91.2%, which was the maximum value of the cubic spline in the histogram of the recovery rate for the full BLS runs. Figure 5 displays the histogram of the recovery rate for the BLS trial planets and the extrapolated likelihoods of detection for the non-BLS runs.

Zoom In Zoom Out Reset image size

6.2. Assessing Pipeline Performance

In total, we injected 5086000 transiting planets into the light curves of the 2543 stars in our sample. We provide a catalog of injected planet parameters and recovery results in Table 3. Our pipeline successfully recovered 86% of signals injected with an expected SNR between 15 and 20. For lower SNRs, the pipeline performance decreased roughly linearly with anticipated SNR until reaching 52% recovery for signals with anticipated SNR between 5 and 7. For the purpose of assessing pipeline performance as a function of SNR, we modeled the anticipated SNR of a transiting planet as

$$\begin{eqnarray}&&\mathrm{SNR}=\displaystyle \frac{\delta }{{\mathrm{CDPP}}_{\mathrm{transit}}}\sqrt{{n}_{\mathrm{transit}}},\end{eqnarray} \tag{ 2 }$$

where δ is the median decrease in brightness during the injected transit, CDPP_transit is the CDPP on the timescale of a full transit of a planet on a circular orbit, and n_transit is the number of transits expected given the orbital period of the planet and the number of days the star was observed. As in Dressing & Charbonneau (2013), we estimated the CDPP on the timescale of a planetary transit by interpolating over the provided CDPP measured on 3, 6, and 12 hr timescales.

Table 3.  Injected Planets

	Period	Radius	Insolation	Recovery
KID	(days)	(R⊕)	(F⊕ )	Statusa
003835071	0.686	0.712	409.813	1.000
005693298	0.876	0.531	250.636	1.000
011560326	1.900	2.700	48.871	0.912
008229458	2.700	0.788	61.212	0.520
004931385	20.420	1.717	3.056	0.650
009691776	55.352	2.555	0.568	0.904
010032631	85.665	1.978	0.664	0.791
002441562	112.865	0.925	0.195	0.000
011013096	155.207	1.795	0.123	0.908
002692704	184.735	1.181	0.225	0.000

Note.

^aFor the 83699 injected planets that were tested in the complete pipeline and the 604278 planets with Δχ² below our 5σ detection threshold, the recovery status indicates whether the planet was detected (1 for recovered planets, 0 for unrecovered planets). For the remaining 4398023 injected planets that had Δχ² above the detection threshold but were not tested in the full BLS simulation, the recovery status indicates the estimated likelihood of detection (see Section 6).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

As a benefit of injecting multiple trial planets per star, we generated unique transit detectability maps for each star in our sample. For example, Figure 6 displays the transit detectability maps for KID 7104554, a Kp = 15.3, T_eff = 3957 K star with lower search completeness than the larger sample. We created the star-by-star transit detectability maps by gridding the injected planets in radius/period and radius/insolation space and calculating the fraction of detectable planets within each grid cell. For the subset of trial planets that passed both Δχ² tests yet were not selected for the full BLS search, we estimated the recovery fraction as explained in Section 6.1.

Zoom In Zoom Out Reset image size

After generating transit detectability maps in radius-period and radius-insolation space for each star independently, we created transit detectability maps for the full sample by summing the individual maps. Although the combined maps displayed in Figure 7 are useful for comparing the sensitivity of any individual star to the sensitivity of the larger sample, the binning is rather coarse. We therefore generated a second set of combined sensitivity maps by sorting the full set of 5086000 injected planets into smaller grid cells in radius/period and radius/insolation space. We then calculated the recovery fractions within each of the cells to produce the smoother sensitivity maps displayed in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in Figure 8, we found that our pipeline is very sensitive to injected planets with radii larger than 2.5R_⊕. Such planets were detected with nearly 90% efficiency out to the maximum injected orbital period of 200 days. Our pipeline had a significantly harder time detecting 1.0–1.5R_⊕ planets with periods longer than 100 days (recovery fraction = 36%) and 0.5–1.0R_⊕ planets with periods longer than 5 days (recovery fraction = 22%). Planets smaller than 1.0R_⊕ were nearly undetectable (recovery rate approximately 6%) at orbital periods longer than 150 days.

Inspecting the transit recovery map as a function of insolation revealed that the transit detectability changes sharply across the HZ. At the inner edge of the HZ (median orbital period of 50 days for the stars in our sample), we recovered 84% of 2.0R_⊕ planets and 34% of 1.0R_⊕ planets. At the outer edge of the HZ (median orbital period of 130 days), the sensitivity decreased to 80% for 2.0R_⊕ planets and 25% for 1.0R_⊕ planets. This large change in sensitivity in a very interesting region of planet radius and insolation space reveals that the pipeline sensitivity within the HZ is not well described by a single number.

6.3. Calculating Search Completeness

The overall planet search completeness depends both on the detectability of a particular transiting planet and the likelihood that a particular planet will be observed to transit. We accounted for the latter factor by determining the mean geometric probability of transit for planets orbiting the stars in our sample at particular periods or insolation levels. For a given orbital period, we computed the corresponding semimajor axis for a planet orbiting each of the stars in our sample. Next, we divided the stellar radii by the calculated semimajor axes to find the transit probability for a planet in a circular orbit.

We then multiplied the transit probability by a correction factor to account for the fact that the planets in our sample are more likely to be eccentric and were accordingly more likely to transit (Barnes 2007; Kipping 2014). We adopted a correction factor of 1.08 based on a beta distribution fit by Kipping (2013) to transiting planets with periods shorter than 382.3 days. Neglecting this correction factor would lead to an underestimate of search completeness and an overestimate of the planet occurrence rate by roughly 8% (Kipping 2014). The resulting search completeness plots are displayed in Figure 9.

Zoom In Zoom Out Reset image size

In order to estimate the planet occurrence rate, we first generated smoothed maps of the detected planet population. For each planet candidate, we counted the number of links from the MCMC posteriors that fell within each grid cell in radius/period and radius/insolation space.⁹ When converting each link of the chains from light curve parameters to physical values, we accounted for uncertainties in the stellar parameters by drawing new stellar parameters from Gaussians centered at the reported values with widths set by the reported errors. (In cases where the reported errors were asymmetric, we adopted the larger value.) We weighted each link so that the total weight equaled one minus the false positive correction (see Section 4.3) for a planet with the given radius.

The errors on the planet radii and insolation flux were large enough that the posteriors from multiple candidates overlapped to produce smoothed distributions. For the orbital periods, however, the errors were small enough that each planet appeared isolated. For the purpose of calculating the planet occurrence rate, we artificially inflated the spread of the period values so that the standard deviation of log₁₀ P distribution was equal to 0.1. We then constructed smoothed distributions of the insolation flux received by each planet by converting the periods into semimajor axes using Kepler's third law and stellar masses drawn from Gaussian distributions centered on the reported value with widths set by the reported errors.

As shown in Figure 10, the detected planet population has peaks within the region P = 1–20 days and R_p = 0.7–2.5R_⊕. There is also a noticeable lack of large planets (${R}_{p}\geqslant 1.7$R_⊕) and shorter periods ($P\leqslant 2$ days). In radius-insolation space (right panel of Figure 10), the highest peaks of the smoothed candidate distribution are located at insolations of 20–50 and 2.5–10 times the insolation received by the Earth.

Zoom In Zoom Out Reset image size

We estimated the planet occurrence rate by dividing the smoothed maps of the detected planet population in Figure 10 by the smoothed maps of the search completeness in Figure 9. The resulting maps of the planet occurrence rate are displayed in Figure 11. The division reveals that the lack of detected planets in the upper left corner of the left panel of Figure 10 is quite meaningful. That region has very high search completeness, so the lack of detected planets in that region of parameter space implies that hot mini-Neptunes and Neptunes are rare around low-mass stars. At the opposite corner of the diagram in the small-planet, long-period regime, the relatively small number of detected planets does not indicate a low occurrence rate. On the contrary, those planets were detected despite relatively low search completeness, so the underlying occurrence rate of such planets is predicted to be high.

Zoom In Zoom Out Reset image size

Consulting the right panel of Figure 11, the estimated occurrence rate of small planets is highest at insolations below roughly 0.5 F_⊕ , but that is a region of low search completeness and the occurrence rate for such planets is not well constrained. In order to more readily see trends in the planet occurrence rate as a function of planetary properties, we binned the smoothed occurrence distributions shown in Figure 11 to produce the gridded diagrams shown in Figure 12. The gridded version of the radius-period diagram (left panel of Figure 12) clearly demonstrates that planet occurrence increases with decreasing planet radius and increasing log₁₀ P.

Zoom In Zoom Out Reset image size

7.1. Dependence on Planet Radius and Period

Figures 13 and 14, and Table 4 display the planet occurrence rate as a function of planet radius and orbital period. As in previous studies (Dressing & Charbonneau 2013; Morton & Swift 2014), we found that planets with radii >3R_⊕ are rare around small stars (at least out to orbital periods of 200 days). Concentrating on planets with periods shorter than 50 days, we observed a general trend of decreasing planet occurrence with increasing planet radius between 1R_⊕ and 4R_⊕. There is an indication in Figure 14 that the planet occurrence rate may become flat in log₁₀ P for periods longer than 10 days. However, the errors on the longest orbital bins are large enough that we cannot distinguish between a brief flattening between 10–100 days and a plateau extending out to much longer orbital periods.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Number of Planets per Star Vs. Orbital Period (In Percentage)

Rp (R⊕)	0.5–1.7 days	1.7–5.5 days	5.5–18.2 days	18.2–60.3 days	60.3–200 days
0.5–1.0	${1.92}_{-0.64}^{+1.01}(70\%)$	${9.88}_{-2.46}^{+3.45}(54\%)$	${23.06}_{-5.57}^{+8.29}(33\%)$	${17.75}_{-6.54}^{+13.34}(17\%)$	(6%)
1.0–1.5	${1.83}_{-0.58}^{+0.87}(87\%)$	${9.18}_{-1.98}^{+2.62}(84\%)$	<26.54 (79%)a	${23.08}_{-6.02}^{+9.38}(64\%)$	${30.70}_{-10.54}^{+26.67}(41\%)$
1.5–2.0	${0.33}_{-0.16}^{+0.46}(90\%)$	${3.70}_{-1.18}^{+1.82}(88\%)$	${20.06}_{-4.04}^{+5.45}(87\%)$	${26.73}_{-6.08}^{+8.99}(84\%)$	${18.90}_{-6.99}^{+14.67}(74\%)$
2.0–2.5	$\lt 0.13(90\%)$	${2.25}_{-0.88}^{+1.57}(90\%)$	${15.69}_{-3.55}^{+4.91}(89\%)$	${23.65}_{-5.83}^{+8.81}(88\%)$	${14.12}_{-5.51}^{+11.55}(84\%)$
2.5–3.0	$\lt 0.11(91\%)$	${1.19}_{-0.54}^{+1.21}(90\%)$	${6.54}_{-2.18}^{+3.51}(90\%)$	${10.42}_{-3.74}^{+6.68}(89\%)$	${5.30}_{-2.52}^{+7.85}(88\%)$
3.0–3.5	<0.11 (91%)	${0.58}_{-0.29}^{+0.91}(91\%)$	${2.55}_{-1.13}^{+2.41}(91\%)$	${2.44}_{-1.19}^{+4.14}(90\%)$	${1.83}_{-0.80}^{+5.32}(89\%)$
3.5–4.0	<0.11 (91%)	${0.36}_{-0.18}^{+0.72}(91\%)$	${0.87}_{-0.43}^{+1.65}(91\%)$	${0.53}_{-0.17}^{+2.31}(90\%)$	<2.71 (90%)
0.5–1.0	${1.38}_{-0.53}^{+0.93}(66\%)$	${8.42}_{-2.39}^{+3.53}(44\%)$	${20.59}_{-5.57}^{+8.70}(26\%)$	(11%)	(3%)
1.0–1.5	${1.95}_{-0.61}^{+0.93}(86\%)$	${9.94}_{-2.13}^{+2.82}(83\%)$	<26.63 (75%)	${26.85}_{-6.79}^{+10.70}(58\%)$	${28.85}_{-10.34}^{+28.66}(32\%)$
1.5–2.0	${0.41}_{-0.20}^{+0.51}(89\%)$	${4.15}_{-1.28}^{+1.94}(88\%)$	<23.58 (86%)	${24.59}_{-6.00}^{+9.08}(82\%)$	${19.98}_{-7.42}^{+16.07}(67\%)$
2.0–2.5	<0.13 (90%)	${2.72}_{-1.01}^{+1.73}(89\%)$	${18.73}_{-3.95}^{+5.39}(89\%)$	${27.58}_{-6.31}^{+9.42}(87\%)$	${18.08}_{-6.57}^{+13.21}(82\%)$
2.5–3.0	<0.11 (90%)	${1.59}_{-0.69}^{+1.39}(90\%)$	${8.29}_{-2.55}^{+3.96}(90\%)$	${14.51}_{-4.62}^{+7.71}(89\%)$	${8.61}_{-3.84}^{+9.78}(87\%)$
3.0–3.5	<0.12 (91%)	${0.65}_{-0.32}^{+1.00}(91\%)$	${3.25}_{-1.37}^{+2.72}(90\%)$	${3.37}_{-1.62}^{+4.62}(90\%)$	${1.97}_{-0.85}^{+5.87}(89\%)$
3.5–4.0	<0.11 (91%)	${0.38}_{-0.19}^{+0.77}(91\%)$	${1.05}_{-0.52}^{+1.82}(91\%)$	${0.56}_{-0.19}^{+2.32}(90\%)$	<2.34 (89%)

Note. In this table and all subsequent occurrence rate tables, the numbers in parentheses are the fraction of injected planets that were recovered within the given intervals. In addition, the first set of entries is our estimates of the planet occurrence rate when using the stellar properties in the Huber et al. (2014) catalog. The second set of entries (below the line) are alternative estimates constructed by revising the stellar radii to lie along an empirical temperature/radius relation from Mann et al. (2013b, see Section 7.4 for details).

^aWe provide 1σ upper limits instead of two-sided errors for grid cells with large Poisson errors and very low occurrence rates.

Download table as:  ASCIITypeset image

For orbital periods shorter than 50 days, we measure an occurrence rate of ${0.56}_{-0.05}^{+0.06}$ Earth-size (1–1.5R_⊕) planets and ${0.46}_{-0.05}^{+0.07}$ super-Earths (1.5–2R_⊕) per small star. Extending the period range to 100 days, we estimate ${0.65}_{-0.05}^{+0.07}$ Earths and ${0.57}_{-0.06}^{+0.08}$ super-Earths per star. Overall, we find 2.5 ± 0.2 planets per M dwarf with radii of 1–4R_⊕ and periods shorter than 200 days. We provide cumulative planet occurrence rates for several additional choices of period and radius boundaries in Table 5.

Table 5.  Cumulative Number of Planets per Star Vs. Orbital Period (In Percentage)

Rp (R⊕)	0.5–10 days	0.5–50 days	0.5–100 days	0.5–150 days	0.5–200 days
0.5–1.0	${22.08}_{-3.79}^{+4.88}(54\%)$	${51.43}_{-6.02}^{+7.99}(44\%)$	${54.73}_{-6.13}^{+8.21}(40\%)$	${59.71}_{-6.21}^{+8.46}(37\%)$	${59.71}_{-6.21}^{+8.44}(36\%)$
1.0–1.5	${21.18}_{-3.11}^{+3.84}(84\%)$	${56.25}_{-5.01}^{+6.31}(79\%)$	${65.09}_{-5.05}^{+6.52}(75\%)$	${82.17}_{-4.24}^{+6.07}(72\%)$	${88.74}_{-3.42}^{+5.38}(70\%)$
1.5–2.0	${10.76}_{-2.40}^{+3.23}(89\%)$	${45.95}_{-5.47}^{+7.05}(87\%)$	${57.35}_{-5.76}^{+7.59}(86\%)$	${66.58}_{-5.72}^{+7.78}(85\%)$	${69.72}_{-5.63}^{+7.77}(84\%)$
2.0–2.5	${8.68}_{-2.26}^{+3.21}(90\%)$	${35.62}_{-5.35}^{+7.02}(89\%)$	${49.81}_{-6.13}^{+8.18}(89\%)$	${54.59}_{-6.26}^{+8.46}(88\%)$	${55.72}_{-6.29}^{+8.52}(88\%)$
2.5–3.0	${3.89}_{-1.38}^{+2.30}(90\%)$	${15.46}_{-3.79}^{+5.41}(90\%)$	${21.97}_{-4.86}^{+6.87}(90\%)$	${23.19}_{-5.05}^{+7.14}(90\%)$	${23.45}_{-5.09}^{+7.20}(90\%)$
3.0–3.5	${1.73}_{-0.77}^{+1.64}(91\%)$	${4.65}_{-1.75}^{+3.10}(90\%)$	${7.15}_{-2.53}^{+4.32}(90\%)$	${7.38}_{-2.60}^{+4.44}(90\%)$	${7.41}_{-2.61}^{+4.45}(90\%)$
3.5–4.0	${0.72}_{-0.36}^{+1.08}(91\%)$	${1.63}_{-0.77}^{+1.86}(91\%)$	${2.11}_{-0.98}^{+2.31}(91\%)$	${2.27}_{-1.05}^{+2.45}(90\%)$	${2.27}_{-1.05}^{+2.46}(90\%)$
0.5–1.0	${18.39}_{-3.74}^{+5.02}(49\%)$	${39.93}_{-6.08}^{+8.22}(38\%)$	${40.60}_{-6.14}^{+8.31}(34\%)$	${40.60}_{-6.12}^{+8.28}(31\%)$	${40.60}_{-6.12}^{+8.27}(30\%)$
1.0–1.5	${22.85}_{-3.31}^{+4.08}(82\%)$	${60.76}_{-5.13}^{+6.57}(75\%)$	${68.38}_{-5.06}^{+6.65}(71\%)$	${85.73}_{-3.89}^{+5.85}(68\%)$	$91.47\pm 5.04(66\%)$
1.5–2.0	${11.96}_{-2.57}^{+3.42}(88\%)$	${45.87}_{-5.48}^{+7.07}(86\%)$	${57.26}_{-5.80}^{+7.65}(85\%)$	${66.74}_{-5.76}^{+7.86}(83\%)$	${70.23}_{-5.66}^{+7.85}(82\%)$
2.0–2.5	${9.98}_{-2.48}^{+3.46}(90\%)$	${42.48}_{-5.67}^{+7.42}(89\%)$	${58.49}_{-6.13}^{+8.27}(88\%)$	${65.45}_{-6.07}^{+8.39}(88\%)$	${67.12}_{-6.03}^{+8.40}(87\%)$
2.5–3.0	${4.94}_{-1.63}^{+2.58}(90\%)$	${20.60}_{-4.44}^{+6.17}(90\%)$	${29.85}_{-5.63}^{+7.79}(90\%)$	${32.40}_{-5.92}^{+8.21}(89\%)$	${33.01}_{-5.99}^{+8.31}(89\%)$
3.0–3.5	${2.09}_{-0.90}^{+1.84}(91\%)$	${6.32}_{-2.22}^{+3.73}(90\%)$	${8.67}_{-2.88}^{+4.74}(90\%)$	${9.13}_{-3.02}^{+4.94}(90\%)$	${9.25}_{-3.05}^{+4.99}(90\%)$
3.5–4.0	${0.83}_{-0.41}^{+1.19}(91\%)$	${1.92}_{-0.89}^{+2.08}(91\%)$	${2.14}_{-0.98}^{+2.27}(90\%)$	${2.20}_{-1.01}^{+2.33}(90\%)$	${2.21}_{-1.02}^{+2.34}(90\%)$

Note. As in Table 4, the entries below the horizontal line are the estimates based on the revised stellar radii (see Section 7.4).

Download table as:  ASCIITypeset image

7.2. Dependence on Planet Radius and Insolation

In Figure 15 and Tables 6 and 7, we present the planet occurrence rate as a function of stellar insolation and planet radius. In general, we find that planet occurrence increases both with decreasing planet radius (as discussed in Section 7.1) and with decreasing log₁₀ F_P. Intriguingly, we observed that the occurrence rates of Earths, Super-Earths, and mini-Neptunes are comparable for insolations below roughly 30 F_⊕ but that planets larger than 1.5 R_⊕ are less common than smaller planets at insolations above 30 F_⊕. The error bars on the coolest insolation bin are rather large, but the divergence of the Earth and Neptune occurrence relations might be due to photo-evaporation at short orbital periods.

Zoom In Zoom Out Reset image size

Table 6.  Number of Planets per Star Vs. Insolation (In Percentage)

Rp (R⊕)	0.2–1.1 F⊕	1.1–6.3 F⊕	6.3–35.6 F⊕	35.6–200 F⊕
0.5–1.0	(11%)	${22.63}_{-6.65}^{+11.24}(24\%)$	${23.48}_{-4.71}^{+6.45}(44\%)$	${3.52}_{-1.10}^{+1.67}(60\%)$
1.0–1.5	${24.02}_{-8.69}^{+19.77}(47\%)$	<35.71 (69%)	${13.45}_{-2.78}^{+3.68}(81\%)$	${5.36}_{-1.24}^{+1.65}(85\%)$
1.5–2.0	${17.70}_{-6.18}^{+11.70}(75\%)$	<31.17 (85%)	${10.59}_{-2.49}^{+3.41}(88\%)$	${1.38}_{-0.53}^{+0.93}(89\%)$
2.0–2.5	${13.97}_{-5.07}^{+9.54}(85\%)$	${20.16}_{-4.73}^{+6.80}(88\%)$	${9.84}_{-2.40}^{+3.34}(89\%)$	${0.75}_{-0.35}^{+0.81}(90\%)$
2.5–3.0	${5.47}_{-2.53}^{+6.67}(88\%)$	${8.90}_{-3.02}^{+5.04}(90\%)$	${3.94}_{-1.38}^{+2.28}(90\%)$	${0.42}_{-0.21}^{+0.62}(90\%)$
3.0–3.5	${1.71}_{-0.78}^{+4.48}(89\%)$	${2.18}_{-1.06}^{+3.15}(90\%)$	${1.71}_{-0.77}^{+1.65}(91\%)$	${0.27}_{-0.13}^{+0.53}(91\%)$
3.5–4.0	<2.33 (90%)	${0.80}_{-0.36}^{+2.28}(91\%)$	${0.57}_{-0.28}^{+1.14}(91\%)$	${0.19}_{-0.09}^{+0.44}(91\%)$
0.5–1.0	(6%)	${10.46}_{-4.32}^{+9.36}(16\%)$	${17.22}_{-4.22}^{+6.10}(32\%)$	${3.23}_{-1.09}^{+1.75}(52\%)$
1.0–1.5	${26.87}_{-9.63}^{+23.57}(37\%)$	${29.27}_{-6.26}^{+9.12}(63\%)$	${13.19}_{-2.83}^{+3.79}(78\%)$	${5.84}_{-1.32}^{+1.75}(84\%)$
1.5–2.0	${18.99}_{-6.69}^{+13.05}(71\%)$	${24.54}_{-5.21}^{+7.34}(83\%)$	${11.35}_{-2.64}^{+3.63}(87\%)$	${1.79}_{-0.64}^{+1.07}(88\%)$
2.0–2.5	${15.29}_{-5.49}^{+10.41}(83\%)$	${22.53}_{-5.05}^{+7.21}(88\%)$	${10.93}_{-2.60}^{+3.59}(89\%)$	${1.05}_{-0.46}^{+0.94}(90\%)$
2.5–3.0	${7.68}_{-3.36}^{+7.80}(87\%)$	${11.83}_{-3.63}^{+5.75}(89\%)$	${4.40}_{-1.51}^{+2.46}(90\%)$	${0.47}_{-0.23}^{+0.65}(90\%)$
3.0–3.5	${1.90}_{-0.87}^{+4.89}(89\%)$	${2.70}_{-1.28}^{+3.34}(90\%)$	${2.11}_{-0.91}^{+1.83}(90\%)$	${0.20}_{-0.10}^{+0.47}(91\%)$
3.5–4.0	<1.81 (90%)	${0.76}_{-0.33}^{+2.22}(90\%)$	${0.75}_{-0.37}^{+1.27}(91\%)$	${0.19}_{-0.09}^{+0.45}(91\%)$

Note. As in Table 4, the entries below the horizontal line are the estimates based on the revised stellar radii (see Section 7.4).

Download table as:  ASCIITypeset image

Table 7.  Cumulative Number of Planets per Star Vs. Insolation (In Percentage)

Rp (R⊕)	0.2–200 F⊕	1.0–200 F⊕	10.0–200 F⊕	50.0–200 F⊕	100.0–200 F⊕	150.0–200 F⊕
0.5–1.0	${63.36}_{-6.34}^{+8.82}(34\%)$	${49.63}_{-6.09}^{+8.11}(41\%)$	${17.79}_{-3.46}^{+4.57}(51\%)$	${2.09}_{-0.74}^{+1.20}(58\%)$	${0.61}_{-0.28}^{+0.61}(61\%)$	${0.06}_{-0.02}^{+0.25}(62\%)$
1.0–1.5	${74.67}_{-4.85}^{+6.59}(70\%)$	${50.65}_{-4.94}^{+6.18}(77\%)$	${15.00}_{-2.46}^{+3.06}(83\%)$	${3.50}_{-0.93}^{+1.30}(85\%)$	${1.00}_{-0.39}^{+0.67}(85\%)$	${0.16}_{-0.08}^{+0.31}(86\%)$
1.5–2.0	${57.61}_{-5.75}^{+7.57}(84\%)$	${39.92}_{-5.09}^{+6.49}(87\%)$	${6.60}_{-1.66}^{+2.31}(88\%)$	${0.90}_{-0.38}^{+0.74}(89\%)$	${0.23}_{-0.11}^{+0.38}(89\%)$	${0.04}_{-0.01}^{+0.20}(89\%)$
2.0–2.5	${44.72}_{-5.76}^{+7.55}(88\%)$	${30.74}_{-4.77}^{+6.17}(89\%)$	${6.41}_{-1.74}^{+2.50}(90\%)$	${0.33}_{-0.17}^{+0.56}(90\%)$	${0.02}_{-0.00}^{+0.19}(90\%)$	${0.00}_{-0.00}^{+0.09}(90\%)$
2.5–3.0	${18.72}_{-4.22}^{+5.92}(90\%)$	${13.26}_{-3.26}^{+4.62}(90\%)$	${2.85}_{-1.03}^{+1.74}(90\%)$	${0.23}_{-0.11}^{+0.46}(90\%)$	${0.04}_{-0.01}^{+0.22}(90\%)$	<0.00 (90%)
3.0–3.5	${5.87}_{-2.10}^{+3.57}(90\%)$	${4.16}_{-1.56}^{+2.72}(90\%)$	${1.28}_{-0.57}^{+1.23}(91\%)$	${0.19}_{-0.09}^{+0.45}(91\%)$	${0.02}_{-0.00}^{+0.19}(91\%)$	<0.00 (91%)
3.5–4.0	${2.03}_{-0.94}^{+2.21}(90\%)$	${1.56}_{-0.73}^{+1.78}(91\%)$	${0.49}_{-0.24}^{+0.79}(91\%)$	${0.18}_{-0.09}^{+0.43}(91\%)$	<0.00 (91%)	<0.00 (91%)
0.5–1.0	${34.08}_{-5.78}^{+7.85}(28\%)$	${30.91}_{-5.43}^{+7.35}(34\%)$	${15.13}_{-3.46}^{+4.79}(44\%)$	${2.03}_{-0.77}^{+1.32}(52\%)$	${0.54}_{-0.25}^{+0.61}(55\%)$	${0.05}_{-0.01}^{+0.23}(57\%)$
1.0–1.5	${75.17}_{-5.00}^{+6.91}(64\%)$	${48.30}_{-5.02}^{+6.31}(73\%)$	${15.61}_{-2.55}^{+3.18}(81\%)$	${3.95}_{-1.01}^{+1.39}(83\%)$	${1.30}_{-0.47}^{+0.77}(84\%)$	${0.22}_{-0.11}^{+0.36}(84\%)$
1.5–2.0	${56.66}_{-5.84}^{+7.72}(81\%)$	${37.67}_{-5.00}^{+6.39}(85\%)$	${7.26}_{-1.77}^{+2.43}(88\%)$	${1.14}_{-0.46}^{+0.84}(88\%)$	${0.28}_{-0.14}^{+0.43}(88\%)$	${0.04}_{-0.01}^{+0.21}(88\%)$
2.0–2.5	${49.80}_{-5.89}^{+7.75}(87\%)$	${34.51}_{-5.01}^{+6.47}(89\%)$	${7.31}_{-1.90}^{+2.68}(89\%)$	${0.54}_{-0.26}^{+0.69}(90\%)$	${0.02}_{-0.00}^{+0.21}(90\%)$	<0.00 (90%)
2.5–3.0	${24.38}_{-4.94}^{+6.81}(89\%)$	${16.71}_{-3.79}^{+5.27}(90\%)$	${3.17}_{-1.13}^{+1.86}(90\%)$	${0.35}_{-0.17}^{+0.56}(90\%)$	${0.05}_{-0.01}^{+0.25}(90\%)$	<0.00 (90%)
3.0–3.5	${6.91}_{-2.36}^{+3.91}(90\%)$	${5.01}_{-1.79}^{+3.03}(90\%)$	${1.51}_{-0.67}^{+1.38}(91\%)$	${0.17}_{-0.08}^{+0.43}(91\%)$	${0.03}_{-0.00}^{+0.22}(91\%)$	<0.00 (91%)
3.5–4.0	${1.93}_{-0.89}^{+2.03}(90\%)$	${1.69}_{-0.78}^{+1.82}(91\%)$	${0.62}_{-0.31}^{+0.92}(91\%)$	${0.18}_{-0.09}^{+0.44}(91\%)$	<0.01 (91%)	${0.00}_{-0.00}^{+0.09}(91\%)$

Note. As in Table 4, the entries below the horizontal line are the estimates based on the revised stellar radii (see Section 7.4).

Download table as:  ASCIITypeset image

Across the size range we considered, the planet occurrence rate versus log₁₀ F_P rises with decreasing insolation between 100–10 F_⊕ and appears roughly flat in log₁₀ F_P between 10–0.2 F_⊕ . Figure 15 hints that the planet occurrence rate might increase again at cooler insolations, but a larger sample of long period planets will be required to test that hypothesis.

7.3. The Occurrence of Potentially Habitable Planets

The shaded regions in Figure 15 display one choice of HZ boundaries, but the definition of the HZ is still rather uncertain. Traditionally, astronomers have used the term to refer to the distance from the star at which liquid water could be present on the surface of a planet (Dole 1964; Hart 1979; Kasting et al. 1993). In theory, there could also be water-based life on worlds with subsurface oceans or non-water-based life on worlds like Titan, but any associated biosignatures would be difficult to interpret remotely. Accordingly, astronomers have concentrated thus far on the search for life as we know it, meaning surface-based lifeforms that depend on liquid water and might generate biosignatures that could alter the composition of their homeworld's atmosphere.

Even within that rather narrow definition, there are many assumptions that can affect the choice of HZ boundaries. In particular, the assumed mass and composition of the planet's atmosphere affects the surface pressure and therefore the temperature range at which water would be liquid (e.g., Vladilo et al. 2013). For instance, Pierrehumbert & Gaidos (2011) showed that planets with thick hydrogen atmospheres would have sufficient surface pressure to retain surface liquid water out to distances of 2.4 AU. It is uncertain whether biosignatures would be detectable in such atmospheres (Hu et al. 2013; Seager et al. 2013), but those worlds could still be habitable.

The presence of clouds adds an additional complication by both cooling and heating the planet. Clouds are particularly important in the case of tidally locked planets, which might be a common fate for planets orbiting within the HZs of M dwarfs. Yang et al. (2013) demonstrated that a tidally locked M dwarf planet might develop a persistent cloud patch above the sub-stellar point. That cloud patch would have a higher albedo than the planetary surface and would allow the planet to be much cooler at a given separation than a cloud-free model would predict. As a result, the HZ for a cloudy, tidally locked planet could extend to insolations as high as F_P = 1.76 F_⊕ for the moist greenhouse limit rather than the limit of F_P < 0.88 F_⊕ calculated by Kopparapu et al. (2013a) for a cloud-free model. Even for non-tidally locked planets, the presence of clouds can expand the distances corresponding to the boundaries of the HZ by roughly 40% depending on the degree of cloud cover (Selsis et al. 2007).

The orbital geometry of exoplanets is also important when assessing planetary habitability. For instance, planets with high obliquities or eccentric orbits might be partially habitable at certain latitudes or during certain times of year (Williams & Kasting 1997; Williams & Pollard 2002; Spiegel et al. 2008, 2009; Dressing et al. 2010; Cowan et al. 2012; Dobrovolskis 2013; Armstrong et al. 2014; Linsenmeier et al. 2015). Depending on the timescale for the temperature of the planet to change (which depends on factors such as the fraction of surface covered by ocean), such planets may undergo periodic global glaciations punctuated by short-lived epochs during which the surface is warm enough for liquid water (Pierrehumbert 2005; Spiegel et al. 2010). In addition, the primordial obliquities of close-in planets orbiting M dwarfs may be significantly eroded by tides (Heller et al. 2011).

Constructing a multi-dimensional HZ model for planets orbiting M dwarfs is beyond the scope of this paper, but we aspire to provide enough information so that other researchers can assess the abundance of planets within their chosen HZ boundaries. We therefore provide occurrence rates for a few possible choices of HZ boundaries in Table 8.

Table 8.  Habitable Zone Occurrence Rates (In Percentage)

FP (F⊕ )	0.25–0.88	0.23–1.54	0.25–4.00	0.25–1.76	0.25–2.78	0.25–5.85
Outer HZ:	Max GHa	Early Marsa	Fixedb	Max GHa	Max GHa	Max GHa
Inner HZ:	Moist GHa	Recent Venusa	Fixedb	Cloudy Moist GHc	Desert (a = 0.2)d	Desert (a = 0.8)d
0.5–1.0 R⊕	(13%)	(14%)	${25.28}_{-7.96}^{+14.90}(18\%)$	${15.75}_{-6.41}^{+15.34}(15\%)$	${20.61}_{-7.37}^{+15.15}(17\%)$	${27.56}_{-7.86}^{+13.70}(20\%)$
0.8–1.0 R⊕	${4.63}_{-1.78}^{+15.19}(21\%)$	${13.09}_{-5.73}^{+16.88}(22\%)$	${20.95}_{-7.41}^{+15.07}(28\%)$	${13.30}_{-5.72}^{+15.31}(24\%)$	${17.23}_{-6.66}^{+14.71}(26\%)$	${20.05}_{-6.99}^{+13.68}(30\%)$
1.0–1.5 R⊕	${15.82}_{-6.54}^{+16.60}(48\%)$	${24.28}_{-8.39}^{+17.58}(50\%)$	${46.77}_{-8.12}^{+12.33}(56\%)$	${21.82}_{-7.54}^{+15.03}(52\%)$	${31.65}_{-8.00}^{+13.20}(55\%)$	${46.90}_{-7.54}^{+11.00}(59\%)$
1.5–2.0 R⊕	${11.54}_{-4.67}^{+9.97}(75\%)$	${20.69}_{-6.32}^{+10.80}(76\%)$	${36.07}_{-6.90}^{+10.00}(79\%)$	${21.48}_{-6.22}^{+10.25}(77\%)$	${28.04}_{-6.62}^{+10.10}(78\%)$	${39.35}_{-6.64}^{+9.32}(80\%)$
2.0–2.5 R⊕	${10.25}_{-4.14}^{+8.58}(85\%)$	${17.09}_{-5.50}^{+9.48}(85\%)$	${29.17}_{-6.33}^{+9.30}(86\%)$	${18.23}_{-5.56}^{+9.27}(86\%)$	${23.99}_{-6.06}^{+9.31}(86\%)$	${30.97}_{-6.23}^{+8.93}(87\%)$
2.5–4.0 R⊕	${5.03}_{-2.37}^{+6.66}(89\%)$	${10.58}_{-4.04}^{+7.73}(89\%)$	${16.61}_{-4.87}^{+7.78}(90\%)$	${11.09}_{-4.10}^{+7.62}(89\%)$	${13.79}_{-4.53}^{+7.71}(89\%)$	${17.96}_{-4.89}^{+7.52}(90\%)$
1.0–2.0 R⊕	${27.36}_{-8.40}^{+15.77}(61\%)$	${44.97}_{-9.29}^{+15.52}(63\%)$	${82.84}_{-5.33}^{+8.99}(68\%)$	${43.31}_{-8.74}^{+14.02}(64\%)$	${59.68}_{-7.96}^{+12.18}(67\%)$	${86.25}_{-4.45}^{+7.51}(69\%)$
2.0–3.0 R⊕	${13.88}_{-5.03}^{+9.47}(86\%)$	${24.44}_{-6.53}^{+10.44}(87\%)$	${41.17}_{-6.96}^{+9.91}(87\%)$	${25.97}_{-6.55}^{+10.21}(87\%)$	${33.87}_{-6.88}^{+10.11}(87\%)$	${44.16}_{-6.76}^{+9.42}(88\%)$
3.0–4.0 R⊕	${1.40}_{-0.58}^{+4.46}(90\%)$	${3.23}_{-1.57}^{+5.08}(90\%)$	${4.61}_{-2.10}^{+5.04}(90\%)$	${3.34}_{-1.62}^{+5.00}(90\%)$	${3.91}_{-1.85}^{+5.02}(90\%)$	${4.77}_{-2.13}^{+4.87}(90\%)$
2.0–4.0 R⊕	${15.28}_{-5.33}^{+9.78}(88\%)$	${27.67}_{-6.86}^{+10.74}(88\%)$	${45.78}_{-7.07}^{+10.02}(89\%)$	${29.32}_{-6.86}^{+10.51}(88\%)$	${37.78}_{-7.09}^{+10.32}(89\%)$	${48.93}_{-6.83}^{+9.50}(89\%)$
0.5–1.4 R⊕	${19.85}_{-7.70}^{+18.42}(28\%)$	${37.34}_{-10.30}^{+20.29}(30\%)$	${67.77}_{-7.97}^{+13.06}(35\%)$	${35.60}_{-9.67}^{+18.05}(31\%)$	${49.48}_{-9.34}^{+15.53}(33\%)$	${69.67}_{-7.32}^{+11.57}(37\%)$
0.5–1.0 R⊕	(8%)	(9%)	(12%)	(10%)	(11%)	(14%)
0.8–1.0 R⊕	(14%)	${4.49}_{-2.03}^{+10.81}(15\%)$	${8.55}_{-3.93}^{+11.74}(20\%)$	${5.54}_{-2.54}^{+11.98}(16\%)$	${7.76}_{-3.61}^{+12.30}(18\%)$	${9.12}_{-4.14}^{+11.70}(22\%)$
1.0–1.5 R⊕	${19.97}_{-7.94}^{+20.84}(38\%)$	${28.20}_{-9.56}^{+21.29}(41\%)$	${49.46}_{-8.85}^{+14.16}(48\%)$	${26.25}_{-8.84}^{+18.49}(43\%)$	${35.71}_{-9.06}^{+15.79}(46\%)$	${48.68}_{-8.30}^{+12.74}(51\%)$
1.5–2.0 R⊕	${13.17}_{-5.25}^{+11.23}(69\%)$	${20.28}_{-6.59}^{+11.89}(71\%)$	${34.10}_{-7.11}^{+10.60}(75\%)$	${20.32}_{-6.36}^{+11.08}(72\%)$	${26.12}_{-6.76}^{+10.73}(74\%)$	${36.79}_{-6.80}^{+9.76}(76\%)$
2.0–2.5 R⊕	${11.62}_{-4.61}^{+9.53}(83\%)$	${18.38}_{-5.90}^{+10.29}(83\%)$	${32.50}_{-6.71}^{+9.85}(85\%)$	${19.41}_{-5.92}^{+9.97}(84\%)$	${26.23}_{-6.44}^{+9.91}(84\%)$	${34.48}_{-6.58}^{+9.41}(85\%)$
2.5–4.0 R⊕	${7.07}_{-3.17}^{+7.75}(88\%)$	${12.65}_{-4.62}^{+8.64}(89\%)$	${21.53}_{-5.61}^{+8.61}(89\%)$	${13.42}_{-4.72}^{+8.52}(89\%)$	${17.56}_{-5.26}^{+8.59}(89\%)$	${23.29}_{-5.61}^{+8.34}(89\%)$
1.0–2.0 R⊕	${33.14}_{-9.60}^{+18.53}(54\%)$	${48.48}_{-10.06}^{+17.89}(56\%)$	${83.55}_{-5.51}^{+9.92}(61\%)$	${46.57}_{-9.54}^{+16.20}(57\%)$	${61.83}_{-8.47}^{+13.69}(60\%)$	${85.47}_{-4.86}^{+8.53}(63\%)$
2.0–3.0 R⊕	${17.14}_{-5.84}^{+10.68}(85\%)$	${28.28}_{-7.16}^{+11.42}(85\%)$	${49.13}_{-7.24}^{+10.35}(86\%)$	${29.93}_{-7.13}^{+11.10}(86\%)$	${39.87}_{-7.35}^{+10.78}(86\%)$	${52.58}_{-6.95}^{+9.75}(87\%)$
3.0–4.0 R⊕	${1.54}_{-0.66}^{+4.72}(89\%)$	${2.75}_{-1.33}^{+5.11}(89\%)$	${4.89}_{-2.20}^{+5.08}(90\%)$	${2.89}_{-1.41}^{+5.02}(89\%)$	${3.92}_{-1.86}^{+5.05}(90\%)$	${5.18}_{-2.26}^{+4.95}(90\%)$
2.0–4.0 R⊕	${18.69}_{-6.13}^{+10.95}(87\%)$	${31.03}_{-7.40}^{+11.64}(87\%)$	${54.03}_{-7.20}^{+10.29}(88\%)$	${32.83}_{-7.35}^{+11.32}(88\%)$	${43.79}_{-7.46}^{+10.88}(88\%)$	${57.77}_{-6.85}^{+9.65}(88\%)$
0.5–1.4 R⊕	${18.48}_{-7.54}^{+20.46}(21\%)$	${29.91}_{-9.83}^{+21.56}(22\%)$	${53.10}_{-9.06}^{+14.77}(27\%)$	${28.96}_{-9.33}^{+19.24}(24\%)$	${39.97}_{-9.52}^{+16.65}(26\%)$	${53.79}_{-8.57}^{+13.47}(30\%)$

Notes. As in Table 4, the entries below the horizontal line are the estimates based on the revised stellar radii (see Section 7.4).

^aThese habitable zone limits are from Kopparapu et al. (2013b). ^bThese are the "simple" habitable zone boundaries adopted by Petigura et al. (2013a). ^cThis limit approximates the effect of clouds (Yang et al. 2014) by increasing the flux at the inner edge of the habitable zone by a factor of two compared to the baseline calculation by Kopparapu et al. (2013b). ^dFor these calculations, the inner edge of the habitable zone was set to the values predicted by Zsom et al. (2013) for hot "desert" worlds with low relative humidity.

Download table as:  ASCIITypeset image

In the most conservative case, we adopt the maximum greenhouse (Max GH) and moist greenhouse (Moist GH) insolation limits from Kopparapu et al. (2013b). The Max GH limit is the insolation at which adding additional CO₂ can no longer heat the surface of the planet because Rayleigh scattering begins to dominate over the greenhouse effect. At the inner edge, the Moist GH limit corresponds to the insolation at which the planet's stratosphere becomes dominated by water vapor. At that point, the planet's reservoir of hydrogen quickly escapes to space.

Table 8 also provides estimates based on the assumption that Venus and Mars were habitable at earlier times in their histories. For a Sun-like star, those constraints correspond to insolation limits of 1.776 F_⊕ and 0.321 F_⊕, respectively (Kopparapu et al. 2013a, 2013b). The insolation boundaries are lower for planets orbiting M dwarfs because the incoming radiation is redder. For a typical star in our sample (T_eff = 3748 K), the boundaries are 1.543 F_⊕ and 0.228 F_⊕, respectively.

Even more optimistically, Table 8 includes HZ occurrence rates using the cloudy inner HZ from Yang et al. (2014), two choices of desert world albedos from Zsom et al. (2013), and the convenient limits of 0.25–4 F_⊕ used by Petigura et al. (2013a). We do not provide an estimate based on the hydrogen atmosphere HZ of Pierrehumbert & Gaidos (2011) because our search completeness is very low at the maximum allowed separation of 2.4 AU.

The appropriate radius range to consider for a potentially habitable planet is more clearly defined than the appropriate insolation range. Based on radial velocity follow-up observations of Kepler planet candidates, Rogers (2015) argued that the majority of planets larger than 1.6 R_⊕ contain too many volatiles to be rocky. This result agrees with previous fits to measured exoplanet masses and radii by Weiss & Marcy (2014) and simulations by Lopez & Fortney (2014). Furthermore, Dressing et al. (2015) found that all five exoplanets smaller than 1.6 R_⊕ with masses and radii measured to a precision better than 20% have densities consistent with an Earth-like mixture of iron and silicates. Like Rogers (2015), they noted that planets larger than 1.6 R_⊕ have densities inconsistent with rocky compositions. Due to observational constraints, the population of planets with well-constrained densities is strongly biased toward highly irradiated planets. We therefore include a broader range of radius choices in Table 8 to account for the possibility that the transition between rocky and gaseous planets might occur at a slightly different radius for less irradiated planets. For instance, the 2.35 R_⊕ exoplanet Kepler-10c has a measured mass of 17.2 ± 1.9M_⊕ and a bulk density of 7.1 ± 1 g cm⁻³, higher than the densities of most 2–3 R_⊕ planets (Dumusque et al. 2014).

Adopting the most conservative assumptions (1.0 R_⊕ < R_P < 1.5 R_⊕, outer HZ = Max GH, inner HZ = Moist GH), we estimate an occurrence rate of ${0.16}_{-0.07}^{+0.17}$ potentially habitable 1–1.5 R_⊕ planets per M dwarf. The predicted occurrence rates of super-Earths (1.5–2.0 R_⊕) and larger planets (2 R_⊕ < R_P < 4 R_⊕) within the same HZ boundaries are ${0.12}_{-0.05}^{+0.10}$ super-Earths and ${0.15}_{-0.05}^{+0.10}$ larger planets per M dwarf. Expanding the radius range to 1–2 R_⊕ or increasing the HZ boundaries to the limits for recent Venus and early Mars increases the assumed occurrence rate to ${0.27}_{-0.08}^{+0.16}$ and ${0.24}_{-0.08}^{+0.18}$, respectively. In the most optimistic case, we estimate an occurrence rate of ${0.86}_{-0.04}^{+0.08}$ desert worlds with albedos of 0.8 and radii of 1–2 R_⊕ receiving insolations between the Zsom et al. (2013) inner limit and the Max GH outer limit. The assumed occurrence rate of potentially habitable M dwarf planets therefore varies by a factor of four depending on the specific choice of radius and insolation boundaries.

The range of HZ possibilities in Table 8 is particularly useful for comparing our results to those of previous studies. For instance, Petigura et al. (2013a) estimated that 22% of FGK stars host 1–2 R_⊕ planets receiving 0.25–4 F_⊕ . Within the same boundaries, we find an occurrence rate of ${0.83}_{-0.05}^{+0.09}$ small planets per M dwarf HZ.

The difference in our estimates might suggest that HZ planets are more common around lower-mass stars, but the Petigura et al. (2013a) prediction is based on an extrapolation of the occurrence rate for shorter period planets to longer period orbits assuming that planet occurrence is flat in log P. If the planet occurrence rate actually increases with log P at longer periods, then perhaps the occurrence rates of potentially habitable planets orbiting FGK and M dwarfs are more similar. Such a change in the slope of the FGK star planet occurrence rate at longer periods could be explained by the radial- and temperature-dependence of the physics governing planet formation.

7.4. Implications of Systematic Biases in Modeled Stellar Radii

The stellar parameters for the majority of the stars in our sample were estimated by fitting Dartmouth stellar models to photometric (Dressing & Charbonneau 2013; Gaidos 2013; Huber et al. 2014) or spectroscopic (Mann et al. 2012; Muirhead et al. 2012a) observations. The exceptions are one star with parameters from Mann et al. (2013b) and two very low-mass stars with parameters from Martín et al. (2013). Both Mann et al. (2013b) and Martín et al. (2013) estimated stellar radii using empirical relations based on interferometric observations of low-mass stars (Boyajian et al. 2012).

Several recent studies (e.g., Boyajian et al. 2012; Mann et al. 2013a; Newton et al. 2015) have demonstrated that theoretical stellar models do not accurately reproduce the observed radii of low-mass stars. As explained in Newton et al. (2015), there are two main issues:

• 1.  
  The radii of model stars with T_eff < 4000 K are smaller than the interferometrically measured radii by approximately 0.04–0.09R_⊙.
• 2.  
  Variations in metallicity produce significant changes in the modeled radii of low-mass stars whereas observations reveal that metallicity actually has little influence on the radii of low-mass stars.

Although the Dartmouth stellar models perform better than many alternative models, both of these effects may have caused the radii of the stars in our sample to be systematically underestimated. In order to gauge the magnitude of this effect, we recalculated the radii for our stellar sample using an empirical temperature/radius relation for main sequence stars with 3300K < T_eff (Equation (6) in Mann et al. 2013a with the additional significant figures reported by Newton et al. 2015). We did not consider changes in the stellar temperatures, but Newton et al. (2015) demonstrated that the temperatures we estimated in Dressing & Charbonneau (2013) were consistent with predictions based on empirical observations (our values were lower by 40 ± 110 K). Using the empirical temperature-radius relation to revise the radii of the 2437 stars in our sample with T_eff > 3300, we found that the median change in radius (ΔR_*) was an increase of 0.026R_⊙ (6%). The change was highly dependent on the assumed metallicity; stars with assigned [Fe/H] $\leqslant -0.5$ displayed a median size increase of 0.05R_⊙ (11%) while the estimated radii of stars with assigned [Fe/H] $\geqslant \ 0$ shrank by 0.016R_⊙ (5%).

For the planet host stars in our sample, the increase in the stellar radii leads to larger predicted radii and increased insolation fluxes for the associated planet candidates. The median planet radius increase was 6.6%, but the amplitude of the change varied considerably. The systems most strongly affected by the revision of the stellar radii were: KOI 3102 (+48%), KOI 2650 (+26%), KOI 2418 (+20%), KOI 2006 (+17%), and KOI 812 (+16%). Two of the KOIs in these systems (KOI 3102.01 and KOI 2650.02) were missed by our planet detection pipeline so they did not enter into our calculation of the planet occurrence rate.

In addition to altering the radius estimates for the detected planet candidates, the changes in the stellar radii affect the estimated survey completeness and in turn, the derived occurrence rate. If the stellar radii are typically 6% larger, then the search completeness we displayed in Figure 9 for planets between 0.50–4.00 R_⊕ actually corresponds to 0.53–4.24 R_⊕ planets. We accounted for this effect by generating new search completeness maps following the procedure outlined in Section 6 but correcting the radii and insolation flux environments of the injected planets to reflect the new radius estimates for each star. We then recalculated the planet occurrence rate using the updated search completeness maps and the revised planet properties. We present the corresponding planet occurrence maps in Figure 16 and include the resulting planet occurrence rates below the horizontal lines in Tables 4– 8.

Zoom In Zoom Out Reset image size

As expected, the most noticeable difference between the occurrence maps displayed in Figures 12 and 16 is that the ridge of high planet occurrence has moved upward to large radii. Similarly, the region of low search completeness now encompasses a slightly larger portion of our chosen parameter space. Using the revised stellar radii, we calculated occurrence rates of ${0.61}_{-0.05}^{+0.07}$ Earth-size planets and ${0.46}_{-0.05}^{+0.07}$ super-Earths per low-mass star with periods shorter than 50 days. These rates are nearly identical to the estimates presented in Section 7. Within the HZ, we estimate a frequency of ${0.20}_{-0.08}^{+0.21}$ Earths and ${0.13}_{-0.05}^{+0.11}$ super-Earths per star when adopting the moist GH inner limit and the maximum GH outer limit from Kopparapu et al. (2013b). These estimates are 26% and 14% higher, respectively, than the rates of ${0.16}_{-0.07}^{+0.17}$ Earths and ${0.12}_{-0.05}^{+0.10}$ super-Earths per HZ presented in Section 7. Although increasing the assumed stellar radii alters the inferred occurrence rates, the dominant source of error is the relatively small number of potentially habitable small planets.

In this paper, we presented an updated estimate of the planet occurrence rate for early M dwarfs based on the full four-year Kepler data set. We developed our own planet detection pipeline to search for transiting planets in the Kepler light curves. We then characterized the completeness of our pipeline by injecting simulated transiting planets into the Kepler light curves and attempting to recover them. Our search of the light curves of 2543 small stars with at least 1000 days of Kepler photometry revealed 3215 possible planetary transits. We thoroughly inspected all available follow-up observations for these objects and accounted for transit depth dilution for systems with close stellar companions. We accepted 156 planet candidates, one of which was not a previously known Kepler planet candidate.

We then measured the occurrence rate of small planets around small stars by dividing smoothed maps of the detected planet population by maps of our pipeline search completeness in radius-period space and radius-insolation space. We found that Earth-sized planets (1.0–1.5 R_⊕) are common and calculated an occurrence rate of ${0.56}_{-0.05}^{+0.06}$ Earth-sized planets with periods shorter than 50 days per early M dwarf. We also found an occurrence of ${0.46}_{-0.05}^{+0.07}$ super-Earths (1.5–2 R_⊕) with periods shorter than 50 days per early M dwarf. For orbital periods shorter than 200 days and planet radii of 1–4 R_⊕, we estimated a cumulative planet occurrence rate of 2.5 ± 0.2 planets per M dwarf.

Within a conservatively defined HZ based on the moist greenhouse and maximum greenhouse limits (Kopparapu et al. 2013a, 2013b) we estimated occurrence rates of ${0.16}_{-0.07}^{+0.17}$ Earth-size (1.0–1.5 R_⊕) planets and ${0.12}_{-0.05}^{+0.10}$ (1.5–2.0 R_⊕) super-Earths per small star. Adopting a wider planet size range of 1–2 R_⊕ and considering the effects of clouds (Yang et al. 2013) increased our estimate to ${0.43}_{-0.09}^{+0.14}$ potentially habitable planets per star. Considering desert worlds (Zsom et al. 2013) would increase the measured occurrence rate to nearly one potentially habitable planet per M dwarf. These estimates span the range of previous estimates of the occurrence rate of potentially habitable M dwarf planets.

An order of magnitude calculation multiplying the occurrence rate of potentially habitable 1–1.5 R_⊕ planets between the empirical early Mars outer boundary and the recent Venus inner boundary by an estimate of the number density of small stars in the galaxy from the RECONS survey (Henry et al. 2006; Winters et al. 2015) therefore suggests that the nearest potentially habitable planet is most likely 2.6 ± 0.4 pc away and is within 3.5 pc with 95% confidence. This estimate assumes that the occurrence rate of potentially habitable planets orbiting later M dwarfs is identical to that for early M dwarfs, which is consistent with the results of Berta et al. (2013).

Correcting for the geometric probability of transit (assuming that 1.4% of potentially habitable M dwarf planets transit), the nearest transiting potentially habitable planet is likely to be ${10.6}_{-1.8}^{+1.6}$ pc away and is within 14.6 pc with 95% confidence.¹⁰ Early M dwarfs at distances of 3.5 and 14.6 pc would have apparent K band magnitudes of 2.9 and 6.0, respectively, well within the magnitude range probed by current and upcoming planet surveys of nearby, bright stars such as CARMENES (Quirrenbach et al. 2010), CHEOPS (Broeg et al. 2013), ExoplanetSat (Smith et al. 2010), ExTrA (Bonfils et al. 2014), HPF (Mahadevan et al. 2010), MEarth (Nutzman & Charbonneau 2008; Berta et al. 2012), PLATO (Rauer et al. 2014), K2 (Howell et al. 2014), SPECULOOS (Gillon et al. 2013), SPIRou (Thibault et al. 2012), and TESS (Ricker et al. 2014).

C.D. is supported by a National Science Foundation Graduate Research Fellowship. Support for this work was provided through the NASA Kepler Mission Participating Scientist Program grants NNX09AB53G and NNX12AC77G, which were awarded to D.C. This publication was made possible through the support of a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. We thank the anonymous referee for providing feedback that improved the quality of this paper and the Kepler team for providing the community with a fantastic collection of data. We are grateful to Jonathan Irwin for sharing a fast implementation of the transit model (Mandel & Agol 2002) and for providing valuable advice. We thank Jessie Christiansen for answering questions about the Kepler pipeline completeness and providing helpful suggestions. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This publication made use of the Kepler Community Follow-Up Observing Program website (https://cfop.ipac.caltech.edu), the NASA Exoplanet Archive, and the Mikulski Archive for Space Telescopes (MAST). The NASA Exoplanet Archive is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.