The Occurrence of Rocky Habitable-zone Planets around Solar-like Stars from Kepler Data

Attempts to measure η_⊕ and general occurrence rates with Kepler have been made since the earliest Kepler catalog releases (Borucki et al. 2011). Youdin (2011) and Howard et al. (2012) were two influential early studies, in which planets found in only the first four months of data (Q0–Q2) were used to constrain Kepler occurrence rates. Youdin (2011) developed a maximum likelihood method to fit an underlying planetary distribution function, which later influenced the Poisson likelihood function method adopted by, e.g., Burke et al. (2015) and Bryson et al. (2020a). Howard et al. (2012) took an alternative approach of estimating occurrence rates in bins of planets, defined over a 2D grid of planet radius and orbital period. In each bin, nondetections were corrected for by weighting each planet by the inverse of its detection efficiency. This inverse detection efficiency method is one of the most popular approaches in the literature.

Catanzarite & Shao (2011) and Traub (2012) (Q0–Q5; Borucki et al. 2011) were among the first papers to focus on the η_⊕ question specifically. Later η_⊕ papers from Dressing & Charbonneau (2013) (Q1–Q6; Batalha et al. 2013), Kopparapu et al. (2013) (Q1–Q6), Burke et al. (2015) (Q1–Q16; Mullally et al. 2016), and Silburt et al. (2015) (Q1–Q16) were able to take advantage of newer planet catalogs based on increased amounts of data. Other papers have used custom pipelines to search Kepler light curves to estimate η_⊕ with independently produced planet catalogs: namely, Petigura et al. (2013) (Q1–Q15), Foreman-Mackey et al. (2014) (Q1–Q15), Dressing & Charbonneau (2015) (Q1–Q16), and Kunimoto & Matthews (2020) (Q1–Q17). Still more have been meta-analyses of results from the exoplanet community based on different Kepler catalogs (Kopparapu 2018; Garrett et al. 2018).

Comparisons between these η_⊕ studies are challenging due to the wide variety of catalogs used, some of which are based on only a fraction of the data of the others. Characterization of completeness has also varied between authors, with some assuming a simple analytic model of detection efficiency (e.g., Youdin 2011; Howard et al. 2012), others empirically estimating detection efficiency with transit injection/recovery tests (e.g., Petigura et al. 2013; Burke et al. 2015), and still others simply assuming a catalog is complete beyond some threshold (e.g., Catanzarite & Shao 2011). Borucki et al. (2011) provided a comprehensive analysis of completeness bias, reliability against astrophysical false positives, and reliability against statistical false alarms based on manual vetting and simple noise estimates. Fully automated vetting has been implemented via the Robovetter (Coughlin 2017) for the Kepler DR24 (Coughlin et al. 2016) and DR25 catalogs. The final Kepler data release (DR25), based on the full set of Kepler observations and accompanied by comprehensive data products for characterizing completeness, has been essential for alleviating issues of completeness and reliability. The DR25 catalog is now the standard used by occurrence rate studies (e.g., Mulders et al. 2018; Hsu et al. 2018; Zink et al. 2019; Bryson et al. 2020a).

DR25 was the first catalog to include data products that allowed for the characterization of catalog reliability against false alarms due to noise and systematic instrumental artifacts, which are the most prevalent contaminants in the η_⊕ regime. Thus nearly all previous works did not incorporate reliability against false alarms in their estimates. Bryson et al. (2020a) were the first to directly take into account reliability against both noise/systematics and astrophysical false positives and in doing so found that occurrence rates for small planets in long-period orbits drop significantly after reliability correction. Mulders et al. (2018) attempted to mitigate the impact of contamination by using a DR25 Disposition Score cut (see Section 7.3.4 of Thompson et al. 2018) as an alternative to reliability correction. As shown in Bryson et al. (2020b), while this approach does produce a higher-reliability planet catalog, explicit accounting for reliability is still necessary for accurate occurrence rates.

Studies have also varied in stellar property catalogs used, and exoplanet occurrence rates have been shown to be sensitive to such choices. For instance, the discovery of a gap in the radius distribution of small planets, first uncovered in observations by Fulton et al. (2017), was enabled by improvements in stellar radius measurements by the California Kepler Survey (Johnson et al. 2017; Petigura et al. 2017). The use of Gaia DR2 parallaxes, which have resulted in a reduction in the median stellar radius uncertainty of Kepler stars from ≈27% to ≈4% (Berger et al. 2020b), has been another significant improvement with important implications for η_⊕. Bryson et al. (2020a) showed that the occurrence rates of planets near Earth's orbit and size can drop by a factor of 2 if one adopts planet radii based on Gaia stellar properties rather than on pre-Gaia Kepler Input Catalog stellar properties.

Measuring η_⊕ requires a definition of what it actually means to be considered a rocky planet in the HZ. Different authors use different definitions, including regarding whether η_⊕ refers to the number of rocky HZ planets per star or to the number of stars with rocky HZ planets. In this paper, for reasons detailed in Section 2, we define η_⊕ as the average number of planets per star with a planet radius between 0.5 and 1.5 Earth radii in the star's HZ, where the average is taken over stars with effective temperatures between 4800 and 6300 K. We compute η_⊕ for both conservative and optimistic HZs, denoted respectively as ${\eta }_{\oplus }^{{\rm{C}}}$ and ${\eta }_{\oplus }^{{\rm{O}}}$.

Most of the existing literature on HZ occurrence rates characterize the HZ using orbital period, where a single period range is adopted to represent the bounds of the HZ for the entire stellar population considered. However, no single period range covers the HZ for a wide variety of stars. Figure 2 shows two example period ranges used for HZ occurrence rate studies relative to the HZ of each star in our stellar parent sample. The SAG13⁴⁹ HZ range of 237 ≤ period ≤ 860 days is shown in blue, and ζ_⊕, defined in Burke et al. (2015) as within 20% of Earth's orbital period, is shown in orange. While these period ranges cover much of the HZ for G stars, they miss significant portions of the HZs of K and F stars and include regions outside the HZ even when restricted to G stars. This will be true for any fixed choice of orbital period range for the range of stellar effective temperatures required for good statistical analysis. Such coverage will not lead to accurate occurrence rates of planets in the HZ. Given that the period ranges of many HZ definitions also extend beyond the detection limit of Kepler, computing η_⊕ requires extrapolation of a fitted population model to longer orbital periods. Lopez & Rice (2018) and Pascucci et al. (2019) presented evidence and theoretical arguments that inferring the population of small rocky planets at low instellation from a population of larger planets at high instellation can introduce significant overestimates of η_⊕.

Zoom In Zoom Out Reset image size

For these reasons, we choose to work in terms of the instellation flux, measuring the photon flux incident on the planet from its host star, rather than the orbital period. In Section 3 we describe how we adopt existing occurrence rate methods and completeness characterizations to use instellation flux instead of orbital period. We address concerns with extrapolating completeness to long orbital periods by providing bounds on the impact of the limited coverage of completeness data (Section 3.3). Following Howard et al. (2012), Youdin (2011), and Burke et al. (2015), among others, we compute the number of planets per star f. As in Youdin (2011) and Burke et al. (2015), we first compute a population model in terms of the differential rate λ ≡ d²f/dr dI, where r is the planet radius and I is the instellation flux. We consider several possible functional forms for λ and allow λ to depend on the stellar host's effective temperature. We compute λ over the radius range 0.5 R_⊕ ≤ r ≤ 2.5 R_⊕ and the instellation flux range 0.2 I_⊕ ≤ I ≤ 2.2 I_⊕, averaged over the effective temperatures of the stellar population used for the computation (Section 3.4). Occurrence rates will be computed by integrating λ over the desired planet radius and instellation flux range and averaging over the desired effective temperature range to give f, the average number of planets per star (Section 3.5).

By restricting our analysis to planets with r ≤ 2.5 R_⊕ in regions of instellation flux close to the HZ, we believe we are avoiding the biases pointed out by Lopez & Rice (2018) and Pascucci et al. (2019). As seen in Figure 1, there are more detected planets with 1.5 R_⊕ ≤ r ≤ 2.5 R_⊕ than planets with 0.5 R_⊕ ≤ r ≤ 1.5 R_⊕, so our results in Section 4 will be driven by these larger planets, but all planets we consider are at similar low levels of instellation. In Figure 2 of Lopez & Rice (2018) we note that for instellation flux between 10 and 20 there is little change in the predicted relative sizes of the 1.5 R_⊕ ≤ r ≤ 2.5 R_⊕ and 0.5 R_⊕ ≤ r ≤ 1.5 R_⊕ planet populations. Naively extrapolating this to instellation of <2 in the HZ, we infer that the sizes of these larger- and smaller-planet populations in the HZ are similar. Therefore by working at low instellation flux we are likely less vulnerable to overestimating the HZ rocky planet population.

We use both Poisson likelihood–based inference with Markov Chain Monte Carlo (MCMC) and likelihood-free inference with approximate Bayesian computation (ABC). The Poisson likelihood method is one of the most common approaches to calculating exoplanet occurrence rates (e.g., Burke et al. 2015; Zink et al. 2019; Bryson et al. 2020a), while ABC was only applied for the first time in Hsu et al. (2018, 2019). As described in Section 3.4.2, these methods differ in their treatment of reliability and input uncertainty and allow us to assess the possible dependence of our result on the assumption of a Poisson likelihood function. We present our results in Section 4. Recognizing the importance of reliability correction in the η_⊕ regime and confirming the same impact of reliability of Bryson et al. (2020a), we choose to report only our reliability-incorporated results. We also present results both with and without incorporating uncertainties in planet radius, instellation flux, and host star effective temperature in our analysis. Our final recommended population models and implications for η_⊕, as well as how our results relate to previous estimates, are discussed in Section 5.

All results reported in this paper are produced with Python code, mostly in the form of Python Jupyter notebooks, found at the paper's GitHub site.⁵⁰

When summarizing a distribution, we use the notation ${m}_{-e2}^{+e1}$ to refer to the 68% credible interval [m − e2, m + e1], where m is the median and "n% credible interval" means that the central n% of the values fall in that interval. We caution our distributions are typically not Gaussian, so this interval should not be treated as "1σ." We will also supply 95% and 99% credible intervals for our most important results.

We use the following notation throughout the paper:

• 1.  
  r: planet radius in units of Earth radii R_⊕
• 2.  
  I: planet instellation flux in units of Earth instellation I_⊕
• 3.  
  T_eff: stellar effective temperature in kelvins. When referring to a planet, this is the effective temperature of that planet's host star.
• 4.  
  f: the number of planets per star, typically a function of r, I, and T_eff
• 5.  
  λ: the differential rate population model ≡ d²f/dr dI, typically a function of r, I, and T_eff. λ is defined by several parameters—for example, exponents when λ is a power law.
• 6.  
  θ: the vector of parameters that define λ, whose contents depend on the particular form of λ

η_⊕ without modification refers to the average number of HZ planets per star with 0.5 R_⊕ ≤ r ≤ 1.5 R_⊕ and a host star effective temperature between 3900 and 6300 K, with the conservative or optimistic HZ specified in that context. ${\eta }_{\oplus }^{{\rm{C}}}$ and ${\eta }_{\oplus }^{{\rm{O}}}$, respectively, specifically refer to occurrence in the conservative and optimistic HZs. Additional subscripts on η_⊕ refer to different stellar populations. For example ${\eta }_{\oplus ,\mathrm{GK}}^{{\rm{C}}}$ is the occurrence of conservative HZ planets with 0.5 R_⊕ ≤ r ≤1.5 R_⊕ around GK host stars.

A key aspect in computing HZ planet occurrence rates is the location and width of the HZ. Classically, it is defined as the region around a star in which a rocky-mass/size planet with an Earth-like atmospheric composition (CO₂, H₂O, and N₂) can sustain liquid water on its surface. Surface liquid water is important for the development of life as we know it, and the availability of water on the surface assumes that any biological activity on the surface alters the atmospheric composition of the planet, betraying the presence of life when observed with remote detection techniques.

Various studies estimate the limits and the width of the HZ in the literature (see Kopparapu 2018 and Kopparapu et al. 2019 for a review) and explore the effect of physical processes such as tidal locking, the rotation rate of the planet, the combination of different greenhouse gases, planetary mass, obliquity, and eccentricity on HZ limits. These effects necessitate a more nuanced approach to identifying habitability limits and are particularly useful for exploring those environmental conditions where habitability could be maintained. However, for the purpose of calculating the occurrence rates of planets in the HZ, it is best to use a standard for HZ limits as a first attempt, such as Earth-like conditions. One reason is that it would become computationally expensive to estimate the occurrence rates of HZ planets considering all the various HZ limits arising from these planetary and stellar properties. Furthermore, future flagship mission concept studies like LUVOIR (The LUVOIR Team 2019), HabEX (Gaudi et al. 2020), and OST (Meixner et al. 2019) used the classical HZ limits as their standard case to estimate exoEarth mission yields and identify associated biosignature gases. Therefore, in this study we use the conservative and optimistic HZ estimates from Kopparapu et al. (2014), where the conservative inner and outer edges of the HZ are defined by the "runaway greenhouse" and "maximum greenhouse" limits and the optimistic inner and outer HZ boundaries are the "recent Venus" and "early Mars" limits. By using these HZ limits, we are able to (1) make a consistent comparison with already published occurrence rates of HZ planets in the literature that have also used the same HZ limits, (2) provide a range of values for HZ planet occurrence, and (3) obtain an "average" occurrence rate of planets in the HZ, as the conservative and optimistic HZ limits from Kopparapu et al. (2014) span the range of HZ limits from more complex models and processes.

We consider planets in the 0.5–1.5 R_⊕ size range to calculate rocky planet occurrence rates, as studies have suggested that planets that fall within these radius bounds are most likely to be rocky (Rogers 2015; Wolfgang et al. 2016; Chen & Kipping 2017; Fulton et al. 2017). While some studies have indicated that the rocky regime can extend to as large as 2.5 R_⊕ (Otegi et al. 2020), many of these large–radius regime planets seem to be highly irradiated planets, receiving stellar fluxes much larger than those of planets within the HZ. Nevertheless, we have also calculated the occurrence rates of planets with radii up to 2.5 R_⊕. We note that Kane et al. (2016) also used Kopparapu et al. (2014) HZ estimates to identify HZ planet candidates using DR24 planet candidate catalog and DR25 stellar properties.

We also limit the host stellar spectral types to stars with 4800 ≤ T_eff ≤ 6300 K, covering mid-K to late F. The reason for limiting to T_eff > 4800 K is twofold: (1) the inner working angle (the smallest angle on the sky at which a direct-imaging telescope can reach its designed ratio of planet to star flux) for the LUVOIR coronagraph instrument ECLIPS falls off below 48 mas at 1 micron (3λ/D) for a planet at 10 pc for T_eff ≤ 4800 K, and (2) planets are likely tidal-locked or synchronously rotating below 4800 K, which could potentially alter the inner HZ limit significantly (Yang et al. 2013, 2014; Godolt et al. 2015; Way et al. 2015; Wolf & Toon 2015; Kopparapu et al. 2016, 2017; Bin et al. 2018). The upper limit of 6300 K is a result of planets in HZs having longer orbital periods around early F stars, where Kepler was not capable of detecting these planets, as described in Section 1.

The width of the HZ for hotter stars is larger than the width for cooler stars, implying that the HZ occurrence rate may be dependent on the host star's effective temperature. In this section we derive an approximate form for this effective temperature dependence, which we refer to as the "geometric effect."

We compute the instellation flux I on a planet orbiting a particular star as $I={R}_{* }^{2}{T}^{4}/{a}^{2}$, where R_* is the stellar radius in solar radii, T = T_eff/T_⊙ is the effective temperature divided by the solar effective temperature, and a is the semimajor axis of the planet orbit in astronomical units. We assume the orbit is circular. Then the size of the HZ Δa is determined by the instellation flux at the inner and outer HZ boundaries, I_inner and I_outer, as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}a & = & {a}_{\mathrm{outer}}-{a}_{\mathrm{inner}}\\ & = & \ {R}_{* }{T}^{2}\left(\displaystyle \frac{1}{\sqrt{{I}_{\mathrm{outer}}}}-\displaystyle \frac{1}{\sqrt{{I}_{\mathrm{inner}}}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

The factor $\left(1/\sqrt{{I}_{\mathrm{outer}}}-1/\sqrt{{I}_{\mathrm{inner}}}\right)$ has a weak T_eff dependence, ranging from 1.25 at 3900 K to 0.97 at 6300 K, which we crudely approximate as constant in this paragraph. We also observe that, for the main-sequence dwarf stellar populations we use in our computations (described in Section 3.1), R_* has an approximately linear dependence on T, which we write as $\left(\tau T+{R}_{0}\right)$ (τ ≈ 1.8 and R₀ ≈ −0.74). Therefore

$$\begin{eqnarray}&&{\rm{\Delta }}a\propto \left(\tau T+{R}_{0}\right){T}^{2}.\end{eqnarray} \tag{ 2 }$$

So even if the differential occurrence rate λ has no dependence on a and therefore no dependence on I, the HZ occurrence rate may depend on T_eff simply because hotter stars have larger HZs.

Several studies, such as Burke et al. (2015) and Bryson et al. (2020a), have examined planet occurrence in terms of the orbital period p and have shown that df/dp is well approximated by a power law p^α. In Appendix A we show that this power law, combined with the relationship between instellation flux and period, implies that the instellation flux portion of the differential rate function λ, df/dI, has the form

$$\begin{eqnarray}&&{df}/{dI}\approx {{CI}}^{\nu }{\left(\left(\tau T+{R}_{0}\right){T}^{4}\right)}^{\delta }\end{eqnarray} \tag{ 3 }$$

where $\nu =-\tfrac{3}{4}\left(\alpha -\tfrac{7}{3}\right)$ and δ = −ν − 1. This form incorporates the combined effects of the size of the HZ increasing with T_eff and dependence on the period power law p^α. The derivation in Appendix A uses several crude approximations, so Equation (3) is qualitative rather than quantitative.

In Section 3.4 we consider forms of the population model λ that separate the geometric effect in Equation (2) from a possible more physical dependence on T_eff, and compare them with a direct measurement of the T_eff dependence. To separate the geometric effect, we incorporate a geometric factor g(T_eff) inspired by Equation (2). Because of the crude approximations used to derive Equations (2) and (3) we use an empirical fit to the HZ width Δa for all stars in our stellar sample. Because we will use this fit in models of the differential population rate df/dI in Section 3.4.1, we perform the fit computing Δa for each star using a fixed flux interval ΔI ∈ [0.25, 1.8]. Because $a={R}_{* }^{2}{T}^{4}/\sqrt{I}$, Δa is just a scaling of each star's luminance ${R}_{* }^{2}{T}^{4}$ by the factor $1/\sqrt{0.25}-1/\sqrt{1.8}$. As shown in Figure 3, Δa is well fit, with well-behaved residuals, by the broken power law

$$\begin{eqnarray}g({T}_{\mathrm{eff}})=\left\{\begin{array}{ll}{10}^{-11.84}\,{T}_{\mathrm{eff}}^{3.16} & \mathrm{if}\ {T}_{\mathrm{eff}}\leqslant 5117\,{\rm{K}},\\ {10}^{-16.77}\,{T}_{\mathrm{eff}}^{4.49} & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

If the semimajor axes of planets are uniformly distributed in our stellar sample, then we expect that HZ planet occurrence would have a T_eff dependence due to Equation (4). In individual planetary systems, however, there is evidence of constant spacing in $\mathrm{log}(a)$ (Weiss et al. 2018), implying spacing proportional to Δa/a. In this case there would be no impact of larger HZs with increasing T_eff: taking a to be the average of the inner and outer semimajor axes, the star-dependent terms cancel each other, so Δa/a is the same for all stars, independently of T_eff. This would imply that HZ occurrence has no T_eff dependence due to the increasing size of the HZ. However, Weiss et al. (2018) also pointed out that the characteristic spacing of systems is at most weakly correlated with stellar type, with stronger correlation of the characteristic spacing with planet radius. A more detailed treatment of planet spacing is necessary to explore how multiplicity influences HZ planet occurrence.

As in Bryson et al. (2020a), our stellar catalog uses the Gaia-based stellar properties from Berger et al. (2020b) combined with the DR25 stellar catalog at the NASA Exoplanet Archive (see Footnote ⁴⁹), with the cuts described in the baseline case of Bryson et al. (2020a). We summarize these cuts here for convenience.

We begin by merging the catalog from Berger et al. (2020b), the DR25 stellar catalog (with its supplement), and the catalog from Berger et al. (2018), keeping only the 177,798 stars that are in all three catalogs. We remove poorly characterized, binary, and evolved stars, as well as stars whose observations were not well suited for long-period transit searches (Burke et al. 2015; Burke & Catanzarite 2017), with the following cuts:

• 1.  
  Remove stars with a Berger et al. (2020b) goodness of fit (iso_gof) of <0.99 and a Gaia renormalized unit weight error (RUWE; Lindegren 2018), as provided by Berger et al. (2020b), of >1.2, leaving 162,219 stars.
• 2.  
  Remove stars that, according to Berger et al. (2018), are likely binaries, leaving 160,633 stars.
• 3.  
  Remove stars that have evolved off the main sequence, recomputing the Evol flag described in Berger et al. (2018) using the Berger et al. (2020b) stellar properties, leaving 105,118 stars.
• 4.  
  Remove noisy targets identified in the KeplerPorts package,⁵¹ leaving 103,626 stars.
• 5.  
  Remove stars with NaN limb darkening coefficients, leaving 103,371 stars.
• 6.  
  Remove stars with a NaN observation duty cycle, leaving 102,909 stars.
• 7.  
  Remove stars with a decrease in observation duty cycle of >30% due to data removal from other transits detected on these stars, leaving 98,672 stars.
• 8.  
  Remove stars with an observation duty cycle of <60%, leaving 95,335 stars.
• 9.  
  Remove stars with a data span of <1000 days, leaving 87,765 stars.
• 10.  
  Remove stars with a DR25 stellar properties table timeoutsumry flag ≠ 1, leaving 82,371 stars.

Selecting FGK stars with effective temperatures between 3900 and 7300 K, which are a superset of the stellar populations we consider in this paper, we have 80,929 stars.

We are primarily interested in HZ occurrence rates for stars with effective temperatures hotter than 4800 K, and Kepler observational coverage is very poor above 6300 K (see Section 2). We fit our population model using two stellar populations and examine the solutions to determine which range of stellar temperatures is best for computing the desired occurrence rate. These stellar populations are as follows:

• 1.  
  hab: stars with effective temperatures between 4800 and 6300 K (61,913 stars)
• 2.  
  hab2: stars with effective temperatures between 3900 and 6300 K (68,885 stars)

The effective temperature distribution of the stars in these populations is shown in Figure 4. This distribution has considerably fewer cooler stars than we believe is the actual distribution of stars in the Galaxy. Our analysis is weighted by the number of stars as a function of effective temperature.

Zoom In Zoom Out Reset image size

There are two stellar population cuts recommended in Burke & Catanzarite (2017) that we do not apply. The first is the requirement that stellar radii be less than 1.35 R_⊙ (1.25 R_⊙ in Burke & Catanzarite 2017, but Burke now recommends 1.35 R_⊙ (private communication)). We do not impose this stellar radius cut, instead opting for a physically motivated selection based on the Evol flag. After our cuts, 6.8% of the hab2 population are stars that have radii larger than 1.35 R_⊙. The completeness analysis for these stars is not expected to be as accurate as that for smaller stars.

The second recommended cut that we do not apply is the requirement that the longest transit duration be less than 15 hr. This cut is due to the fact planet search in the Kepler pipeline does not examine transit durations longer than 15 hr (Twicken et al. 2016). For the hab2 population, assuming circular orbits, the transit durations of planets at the inner optimistic HZ boundary exceed 15 hr for 2.7% of the stars. The transit durations of planets at the outer optimistic HZ boundary exceed 15 hr for 35% of the stars, with the duration being less than 25 hr for 98.7% of the stars. While the exact impact on completeness analysis of transit durations longer than 15 hr is unknown, there is evidence that the impact is small. KOI 5236.01, for example, has a transit duration of 14.54 hr, an orbital period of 550.86 days, and a transit signal-to-noise ratio (S/N) of 20.8. KOI 5236.01 is correctly identified in the Kepler pipeline in searches for transit durations of 3.5–15 hr. KOI 7932.01 has a transit duration of 14.84 hr, an orbital period of 502.256 days, and a transit S/N of 8.1, among the lowest transit S/N for planet candidates with periods of >450 days. KOI 7932.01 is correctly identified in searches using transit durations of 9–15 hr. So even for low-S/N transits, the transit can be identified in searches for transit durations 9–15 times the actual duration. If these examples are typical, we can expect that transit durations of up to 25 hr will be detected. While these examples do not show that the impact of long transits on completeness is actually small, the bulk of these long durations occur in orbits beyond 500 days, so they are absorbed by the upper and lower bounds on completeness we describe in Section 3.3.2. We are confident that long transit durations, as well as large stars, cause completeness to decrease, so their impact falls within these upper and lower bounds.

We use the planet properties from the Kepler DR25 threshold-crossing event (TCE) catalog (Twicken et al. 2016), with the Gaia-based planet radii and instellation fluxes from Berger et al. (2020a).

Three DR25 small-planet candidates that are near their host star's HZ (planet radius of ≤2.5 R_⊕ and instellation flux between 0.2 and 2.2 I_⊕) are not included in our planet sample. KOI 854.01 and KOI 4427.01 are orbiting host stars with effective temperatures of ≤3900 K, and KOI 7932.01's host star is cut from our stellar populations because it is marked "evolved" (see Bryson et al. 2020a).

3.3.1. Detection and Vetting Completeness

3.3.2. Completeness Extrapolation

For most stars in our parent sample, there are regions of the HZ that require orbital periods beyond the 500 day limit of the period–radius completeness contours. Figure 5 shows the distribution of orbital periods at the inner and outer optimistic HZ boundaries for FGK stars in our stellar sample relative to the 500 day limit. We see that a majority of these stars will require some completeness extrapolation to cover their HZs and a small fraction of stars have no completeness information at all. It is unknown precisely how the completeness contours will extrapolate out to longer periods, but we believe that the possible completeness values can be bounded.

Zoom In Zoom Out Reset image size

We assume that completeness is, on average, a decreasing function of orbital period. Therefore, the completeness beyond 500 days will be less than the completeness at 500 days. While this may not be a correct assumption for a small number of individual stars due to local completeness minima in the period due to the window function (Burke & Catanzarite 2017), we have high confidence that this assumption is true on average. We therefore bound the extrapolated completeness for each star by computing the two extreme extrapolation cases:

• 1.  
  Constant-completeness extrapolation, where, for each radius bin, completeness for periods greater than 500 days is set to the completeness at 500 days. This extrapolation will have higher completeness than reality, resulting in a smaller completeness correction and lower occurrence rates, which we take to be a lower bound. In the tables below we refer to this lower bound as "low" values. Here "low" refers to the resulting occurrence rates, and some population model parameters in the this case will have higher values.
• 2.  
  Zero-completeness extrapolation, where, for each radius bin, completeness for periods greater than 500 days is set to zero. Zero completeness will have lower completeness than reality, resulting in a larger completeness correction and higher occurrence rates, which we take to be an upper bound. In the tables below we refer to this upper bound as "high" values. Here "high" refers to the resulting occurrence rates, and some population model parameters in this case will have lower values.

We solve for population models and compute occurrence rates for both extrapolation cases. Figure 6 shows the relative difference in the completeness contours summed over all stars. We see that for effective temperatures below ∼4500 K the difference between constant- and zero-completeness extrapolation is very close to zero, because these cooler stars are well covered by completeness data (see Figure 1). We therefore expect the upper and lower occurrence rate bounds to converge for these stars.

Zoom In Zoom Out Reset image size

The Poisson likelihood we use requires the completeness contours to be summed over all stars, while the ABC method requires the completeness to be averaged over all stars. We observe a significant dependence of summed completeness on effective temperature, shown in Figure 7. We address this dependence of completeness on effective temperature by summing (for the Poisson likelihood) or averaging (for ABC) the completeness contours in effective temperature bins, as described in Section 3.4.2.

Zoom In Zoom Out Reset image size

3.3.3. Reliability

As described in Section 1.2, we develop a planet population model using a parameterized differential rate function λ and use Bayesian inference to find the model parameters that best explain the data. To test the robustness of our results, we use both the Poisson likelihood MCMC method of Burke et al. (2015) and the ABC method of Kunimoto & Matthews (2020) to compute our population model. Both methods are modified to account for vetting completeness and reliability, with the Poisson likelihood method described in Bryson et al. (2020a) and the ABC method described in Kunimoto & Bryson (2020) and Bryson et al. (2020b). In contrast to previous work, we also take into account uncertainties in planet radius, instellation flux, and host star effective temperature, described in Section 3.4.2.

3.4.1. Population Model Choices

We consider three population models for the differential rate function λ(I, r, T). These models are functions of instellation flux I, planet radius r, and stellar effective temperature T_eff. These models depend on possibly different sets of parameters, which we describe with the parameter vector θ. For each model, we solve for the θ that best describes the planet candidate data:

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{1}(I,r,T,{\boldsymbol{\theta }}) & = & {F}_{0}{C}_{1}{r}^{\alpha }{I}^{\beta }{T}^{\gamma }g(T),{\boldsymbol{\theta }}=\left({F}_{0},\alpha ,\beta ,\gamma \right)\\ {\lambda }_{2}(I,r,T,{\boldsymbol{\theta }}) & = & {F}_{0}{C}_{2}{r}^{\alpha }{I}^{\beta }{T}^{\gamma },{\boldsymbol{\theta }}=\left({F}_{0},\alpha ,\beta ,\gamma \right)\\ {\lambda }_{3}(I,r,T,{\boldsymbol{\theta }}) & = & {F}_{0}{C}_{3}{r}^{\alpha }{I}^{\beta }g(T),{\boldsymbol{\theta }}=\left({F}_{0},\alpha ,\beta \right)\end{array}\end{eqnarray} \tag{ 5 }$$

where g(T) is given by Equation (4). The normalization constants C_i in Equation (5) are chosen so that the integral of λ from r_min to r_max and I_min to I_max, averaged over T_min to T_max, = F₀, so F₀ is the average number of planets per star in that radius, instellation flux, and effective temperature range.

λ₁ allows for dependence on T_eff beyond the geometric dependence described in Section 2.2, breaking possible degeneracy between any intrinsic T_eff and the geometric dependence by fixing the geometric dependence as g(T). So, for example, if the planet population rate's dependence is entirely due to the larger HZ for hotter stars, captured in λ₁ by g(T), then there is no additional T_eff dependence and γ = 0. λ₂ does not separate out the geometric T_eff dependence. λ₃ assumes that there is no T_eff dependence beyond the geometric effect.

All models and inference calculations use uniform uninformative priors: 0 ≤ F₀ ≤ 50,000, −5 ≤ α ≤ 5, −5 ≤ β ≤ 5, and −500 ≤ γ ≤ 50. The computations are initialized to a neighborhood of the maximum likelihood solution obtained with a standard nonlinear solver.

3.4.2. Inference Methods

Both the Poisson likelihood and ABC inference methods use the same stellar and planet populations and the same characterization of completeness and reliability computed using the approach of Bryson et al. (2020a). These steps are as follows:

• 1.  
  Select a 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. We consider two such populations, defined by the effective temperature range as described in Section 3.1, in order to explore the dependence of our results on the choice of stellar population.
• 2.  
  Use the injected data to characterize vetting completeness.
• 3.  
  Compute the detection completeness using a version of KeplerPorts⁵² modified for vetting completeness and insolation interpolation, incorporating vetting completeness and geometric probability for each star, and sum over the stars in effective temperature bins, as described in Section 3.3.1.
• 4.  
  Use observed, inverted, and scrambled data to characterize false-alarm reliability, as described in Section 3.3.3.
• 5.  
  Assemble the collection of planet candidates, including by computing the reliability of each candidate from the false-alarm reliability and false-positive probability.
• 6.  
  For each model in Equation (5), use the Poisson likelihood or ABC method to infer the model parameters θ that are most consistent with the planet candidate data for the selected stellar population.

Because vetting completeness and reliability depend on the stellar population and the resulting planet catalog, all steps are computed for each choice of stellar population.

A full study of the impact of uncertainties in stellar and planet properties would include the impact of uncertainties on detection contours and is beyond the scope of this paper. However, we study the impact of uncertainties in planet radius, instellation flux, and host star effective temperature, shown Figure 1, on our occurrence rates. For both the Poisson likelihood and ABC methods we perform our inference computation both with and without uncertainties, which allows us to estimate the approximate contribution of input planet property uncertainties to the final occurrence rate uncertainties. In Bryson et al. (2020a) it was shown that uncertainties in reliability characterization have effectively no impact.

For the Poisson likelihood inference of the parameters in Equation (5) without input uncertainty, reliability is implemented by running the MCMC computation 100 times, with the planets removed with a probability given by their reliability. The likelihood used in the Poisson method is Equation (17) in Appendix B. For details see Bryson et al. (2020a).

We treat input uncertainties in a manner similar to how we treat reliability: we run the Poisson MCMC inference 400 times, each time selecting the planet population according to reliability. We then sample the planet instellation flux, radius, and star effective temperature from the two-sided normal distribution with the width given by the respective catalog uncertainties. We perform this sampling prior to restricting to our period and instellation flux ranges of interest so planets whose median property values are outside the range may enter the range due to their uncertainties. The posteriors from the 400 runs are concatenated together to give the posterior distribution of the parameters θ for each model. This approach to uncertainty does not recompute the underlying parent stellar population with resampled effective temperature uncertainties, because that would require recomputation of the completeness contours with each realization, which is beyond our computational resources. Shabram et al. (2020) performed a similar uncertainty study, properly resampling the underlying parent population, and observed an impact of uncertainty similar to ours (see Section 4.2). Our analysis of uncertainty should be considered an approximation. While the result is not technically a sample from a posterior distribution, in Section 4 we compare the resulting sample to the posterior sample from the model neglecting these uncertainties and find that the population parameter values and resulting occurrence rates change in a predictable way.

The ABC-based inference of the parameters in Equation (5) is computed using the approach of Kunimoto & Bryson (2020), with some modifications to accommodate temperature dependence and uncertainties on planet radius, instellation flux, and temperature.

In the ABC method, the underlying Kepler population is simulated in each completeness effective temperature bin separately. N_p = F₀N_sh(T) planets are drawn for each bin, where N_s is the number of stars in the bin and h(T) collects the model-dependent temperature terms from Equation (5), averaged over the temperature range of the bin and normalized over the entire temperature range of the sample. Following the procedure of Mulders et al. (2018), we assign each planet an instellation flux between 0.2 and 2.2 I_⊕ from the cumulative distribution function of I^β and a radius between 0.5 and 2.5 R_⊕ from the cumulative distribution function of r^α. The detectable planet sample is then simulated from this underlying population by drawing from a Bernoulli distribution with a detection probability averaged over the bin's stellar population. We compare the detected planets to the observed planet candidate population using a distance function, which quantifies agreement between the flux distributions, radius distributions, and sample sizes of the catalogs. For the distances between the flux and radius distributions, we choose the two-sample Anderson–Darling (AD) statistic, which has been shown to be more powerful than the commonly used Kolmogorov–Smirnoff test (Engmann & Cousineau 2011). The third distance is the modified Canberra distance from Hsu et al. (2019),

$$\begin{eqnarray}&&\rho =\displaystyle \sum _{i=1}^{N}\frac{| {n}_{s,i}-{n}_{o,i}| }{\sqrt{{n}_{s,i}+{n}_{o,i}}}\end{eqnarray} \tag{ 6 }$$

where n_{s, i} and n_o,i are the number of simulated and observed planets within the ith bin's temperature range and the sum is over all N bins. For more details, see Bryson et al. (2020b).

These simulations are repeated within a population Monte Carlo ABC algorithm to infer the parameters that give the closest match between the simulated and observed catalogs. With each iteration of the ABC algorithm, model parameters are accepted when each resulting population's distance from the observed population is less than the 75th quantile of the previous iteration's accepted distances. Following the guidance of Prangle (2017), we confirm that our algorithm converges by observing that the distances between the simulated and observed catalogs approach zero with each iteration, and see that the uncertainties on the model parameters flatten out to a noise floor.

This forward model is appropriate for estimating the average number of planets per star in a given flux, radius, and temperature range, similar to the Poisson likelihood method. However, rather than requiring many inferences on different catalogs to incorporate reliability or input uncertainty, we take a different approach. For reliability, we modify the distance function as described in Bryson et al. (2020b). In summary, we replace the two-sample AD statistic with a generalized AD statistic developed in Trusina et al. (2020) that can accept a weight for each data point, and set each observed planet's weight equal to its reliability. We also alter the third distance (Equation (6)) so that a planet's contribution to the total number of planets in its bin is equal to its reliability. As demonstrated in Kunimoto & Bryson (2020) and Bryson et al. (2020b), this weighted distance approach gives results consistent with the Poisson likelihood function method with reliability. Meanwhile, to account for input uncertainty, the observed population is altered for every comparison with a simulated population by randomly assigning each observed planet a new radius, instellation flux, and host star effective temperature from the two-sided normal distribution with the width given by their respective uncertainties.

Once the population rate model λ has been chosen and its parameters determined as described in Section 3.4, we can compute the number of HZ planets per star. For planets with radii between r₀ and r₁ and instellation fluxes between I₀ and I₁, for a star with effective temperature T_eff the number of planets per star is

$$\begin{eqnarray}&&f({T}_{\mathrm{eff}})={\int }_{{r}_{0}}^{{r}_{1}}{\int }_{{I}_{0}}^{{I}_{1}}\lambda (I,r,{T}_{\mathrm{eff}},{\boldsymbol{\theta }})\,{dI}\,{dr}.\end{eqnarray} \tag{ 7 }$$

For a collection of stars with effective temperatures ranging from T₀ to T₁, we compute the average number of planets per star, assuming a uniform distribution of stars in that range, as

$$\begin{eqnarray}&&f=\displaystyle \frac{1}{{T}_{1}-{T}_{0}}{\int }_{{T}_{0}}^{{T}_{1}}f(T)\,{dT}.\end{eqnarray} \tag{ 8 }$$

We typically compute Equation (8) for every θ in the posterior of our solution, giving a distribution of occurrence rates.

The HZ is not a rectangular region in the I–T_eff plane (see Figure 1), so to compute occurrence in the HZ for a given T_eff, we integrate I from the inner HZ flux I_out(T_eff) to the outer flux I_in(T_eff):

$$\begin{eqnarray}&&{f}_{\mathrm{HZ}}({T}_{\mathrm{eff}})={\int }_{{r}_{\min }}^{{r}_{\max }}{\int }_{{I}_{\mathrm{out}}({T}_{\mathrm{eff}})}^{{I}_{\mathrm{in}}({T}_{\mathrm{eff}})}\lambda (I,r,{T}_{\mathrm{eff}},{\boldsymbol{\theta }})\,{dI}\,{dr}.\end{eqnarray} \tag{ 9 }$$

The functions I_out(T) and I_in(T) are given in Kopparapu et al. (2014) and depend on the choice of conservative versus optimistic HZ. f_HZ(T_eff) will be a distribution of occurrence rates if we use a distribution of θ. For a collection of stars with effective temperatures ranging from T₀ to T₁, we compute f_HZ(T_eff) for a sampling of T ∈ [T₀, T₁] and concatenate these distributions together to make a distribution of HZ occurrence rates f_HZ for that radius, flux, and temperature range. When we compute f_HZ to determine the occurrence rate for a generic set of stars, we uniformly sample over [T₀, T₁] (in practice we use all integer kelvin values of T ∈ [T₀, T₁]). The resulting distribution is our final result.

Figure 8 shows the impact of uncertainty in stellar effective temperature on the HZ boundaries. For each star we compute the uncertainty in the HZ boundaries with 100 realizations of that star's effective temperature with uncertainty, modeled as a two-sided Gaussian. The gray regions in Figure 8 show the 86% credible intervals of the uncertainty of the HZ boundaries. These intervals are small relative to the size of the HZ and are well centered on the central value. For example, consider the inner optimistic HZ boundary, which has the widest error distribution in Figure 8. The median of the difference between the median HZ uncertainty and the HZ boundary without uncertainty is less than 0.002%, with a standard deviation less than 0.9%. Therefore, we do not believe that uncertainties in HZ boundaries resulting from stellar effective temperature uncertainties have a significant impact on occurrence rates.

Zoom In Zoom Out Reset image size

For each choice of population differential rate model from Equation (5) and stellar population from Section 3.1, we determine the parameter vector θ with zero and constant extrapolated completeness, giving high and low bounds on the occurrence rates. These solutions are computed over the radius range 0.5 R_⊕ ≤ r ≤ 2.5 R_⊕ and the instellation flux range 0.2 I_⊕ ≤ I ≤ 2.2 I_⊕ using the hab and hab2 stellar populations described in Section 3.1. We perform these calculations both without and with input uncertainties in the planet radius, instellation flux, and T_eff shown in Figure 1. Examples of the resulting θ posterior distributions are shown in Figure 9 and the corresponding rate functions λ are shown in Figure 10. The solutions for models 1–3 for the hab and hab2 stellar populations computed using the Poisson likelihood method are given in Table 1, and those computed using ABC are given in Table 2. An example of the sampled planet population using input uncertainties is shown in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Parameter Fits with 68% Confidence Limits for Models 1–3 from Equation (5) for the hab and hab2 Stellar Populations from Section 3.1, Computed with the Poisson Likelihood Method

	With Uncertainty		Without Uncertainty
	Based on hab Stars	Based on hab2 Stars	Based on hab Stars	Based on hab2 Stars
	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound
	Model 1
F0	$+{1.08}_{-0.57}^{+1.56}$ – $+{1.97}_{-1.17}^{+3.73}$	$+{1.11}_{-0.44}^{+0.88}$ – $+{1.59}_{-0.7}^{+1.56}$	$+{0.70}_{-0.32}^{+0.77}$ – $+{1.41}_{-0.77}^{+2.12}$	$+{1.02}_{-0.37}^{+0.66}$ – $+{1.46}_{-0.59}^{+1.18}$
α	$-{1.05}_{-1.2}^{+1.41}$ – $-{1.09}_{-1.18}^{+1.36}$	$-{1.08}_{-0.85}^{+0.94}$ – $-{1.18}_{-0.87}^{+0.96}$	$-{0.29}_{-1.18}^{+1.39}$ – $-{0.51}_{-1.15}^{+1.35}$	$-{0.96}_{-0.74}^{+0.83}$ – $-{1.03}_{-0.77}^{+0.83}$
β	$-{0.56}_{-0.42}^{+0.48}$ – $-{1.18}_{-0.56}^{+0.6}$	$-{0.84}_{-0.3}^{+0.32}$ – $-{1.19}_{-0.36}^{+0.37}$	$-{0.43}_{-0.4}^{+0.46}$ – $-{1.13}_{-0.5}^{+0.54}$	$-{0.78}_{-0.28}^{+0.3}$ – $-{1.15}_{-0.33}^{+0.34}$
γ	$-{1.84}_{-3.39}^{+3.33}$ – $+{0.91}_{-3.88}^{+3.87}$	$-{2.67}_{-1.57}^{+1.59}$ – $-{1.38}_{-1.78}^{+1.84}$	$-{2.13}_{-3.13}^{+3.06}$ – $+{0.25}_{-3.47}^{+3.39}$	$-{2.33}_{-1.46}^{+1.47}$ – $-{1.03}_{-1.64}^{+1.66}$
	Model 2
F0	$+{1.04}_{-0.55}^{+1.52}$ – $+{1.96}_{-1.18}^{+3.72}$	$+{1.13}_{-0.46}^{+0.92}$ – $+{1.50}_{-0.66}^{+1.43}$	$+{0.80}_{-0.38}^{+0.9}$ – $+{1.33}_{-0.71}^{+1.98}$	$+{1.00}_{-0.38}^{+0.7}$ – $+{1.43}_{-0.6}^{+1.21}$
α	$-{1.03}_{-1.23}^{+1.41}$ – $-{1.05}_{-1.18}^{+1.45}$	$-{1.13}_{-0.86}^{+0.97}$ – $-{1.12}_{-0.86}^{+0.94}$	$-{0.48}_{-1.17}^{+1.36}$ – $-{0.47}_{-1.16}^{+1.33}$	$-{0.92}_{-0.78}^{+0.88}$ – $-{1.01}_{-0.8}^{+0.88}$
β	$-{0.56}_{-0.42}^{+0.48}$ – $-{1.20}_{-0.56}^{+0.6}$	$-{0.85}_{-0.3}^{+0.32}$ – $-{1.18}_{-0.35}^{+0.37}$	$-{0.51}_{-0.4}^{+0.45}$ – $-{1.09}_{-0.51}^{+0.55}$	$-{0.80}_{-0.28}^{+0.31}$ – $-{1.18}_{-0.33}^{+0.34}$
γ	$+{2.60}_{-3.61}^{+3.56}$ – $+{5.33}_{-4.02}^{+3.94}$	$+{1.19}_{-1.63}^{+1.64}$ – $+{2.31}_{-1.85}^{+1.91}$	$+{2.03}_{-3.23}^{+3.1}$ – $+{4.66}_{-3.51}^{+3.43}$	$+{1.26}_{-1.54}^{+1.54}$ – $+{2.73}_{-1.71}^{+1.73}$
	Model 3
F0	$+{1.13}_{-0.58}^{+1.52}$ – $+{1.83}_{-1.0}^{+2.76}$	$+{1.41}_{-0.59}^{+1.14}$ – $+{1.89}_{-0.78}^{+1.51}$	$+{0.89}_{-0.4}^{+0.94}$ – $+{1.24}_{-0.6}^{+1.49}$	$+{1.25}_{-0.5}^{+0.93}$ – $+{1.75}_{-0.7}^{+1.26}$
α	$-{1.08}_{-1.18}^{+1.39}$ – $-{1.06}_{-1.18}^{+1.38}$	$-{1.37}_{-0.83}^{+0.91}$ – $-{1.28}_{-0.82}^{+0.9}$	$-{0.60}_{-1.12}^{+1.3}$ – $-{0.38}_{-1.18}^{+1.34}$	$-{1.21}_{-0.78}^{+0.88}$ – $-{1.17}_{-0.76}^{+0.85}$
β	$-{0.56}_{-0.41}^{+0.48}$ – $-{1.16}_{-0.51}^{+0.55}$	$-{0.89}_{-0.29}^{+0.32}$ – $-{1.29}_{-0.32}^{+0.35}$	$-{0.49}_{-0.39}^{+0.45}$ – $-{1.06}_{-0.48}^{+0.53}$	$-{0.83}_{-0.28}^{+0.3}$ – $-{1.26}_{-0.31}^{+0.33}$

Notes. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. "With Uncertainty" means the planet candidate radius, instellation flux, and host star effective temperature uncertainties are taken into account.

Download table as:  ASCIITypeset image

Table 2.  Parameter Fits with 68% Confidence Limits for Models 1–3 from Equation (5) for the Four Stellar Populations from Section 3.1, Computed with the ABC Method

	With Uncertainty		Without Uncertainty
	Based on hab Stars	Based on hab2 Stars	Based on hab Stars	Based on hab2 Stars
	Low Bound–High Bound	Low Bound—High Bound	Low Bound–High Bound	Low Bound–High Bound
	Model 1
F0	$+{1.18}_{-0.56}^{+0.95}$ – $+{2.04}_{-0.99}^{+1.44}$	$+{1.17}_{-0.44}^{+0.78}$ – $+{1.61}_{-0.65}^{+1.05}$	$+{0.73}_{-0.29}^{+0.54}$ – $+{1.37}_{-0.61}^{+1.08}$	$+{0.94}_{-0.32}^{+0.45}$ – $+{1.41}_{-0.57}^{+0.99}$
α	$-{1.14}_{-0.89}^{+1.02}$ – $-{0.95}_{-0.86}^{+0.99}$	$-{1.14}_{-0.77}^{+0.75}$ – $-{1.18}_{-0.67}^{+0.72}$	$-{0.17}_{-0.97}^{+1.19}$ – $-{0.11}_{-0.88}^{+1.17}$	$-{0.71}_{-0.68}^{+0.67}$ – $-{0.83}_{-0.74}^{+0.77}$
β	$-{0.69}_{-0.38}^{+0.41}$ – $-{1.32}_{-0.44}^{+0.51}$	$-{0.90}_{-0.26}^{+0.31}$ – $-{1.26}_{-0.31}^{+0.35}$	$-{0.67}_{-0.35}^{+0.38}$ – $-{1.30}_{-0.41}^{+0.43}$	$-{0.89}_{-0.25}^{+0.25}$ – $-{1.26}_{-0.32}^{+0.30}$
γ	$-{0.84}_{-4.11}^{+3.81}$ – $+{2.16}_{-3.68}^{+3.91}$	$-{2.60}_{-1.84}^{+1.74}$ – $-{1.14}_{-2.02}^{+2.15}$	$-{1.93}_{-3.37}^{+3.54}$ – $+{1.82}_{-3.77}^{+3.94}$	$-{2.27}_{-1.71}^{+1.65}$ – $-{0.78}_{-1.91}^{+2.11}$
	Model 2
F0	$+{1.06}_{-0.48}^{+0.90}$ – $+{1.88}_{-0.87}^{+1.38}$	$+{1.14}_{-0.43}^{+0.74}$ – $+{1.66}_{-0.70}^{+1.26}$	$+{0.70}_{-0.28}^{+0.44}$ – $+{1.22}_{-0.57}^{+0.87}$	$+{0.96}_{-0.34}^{+0.45}$ – $+{1.35}_{-0.51}^{+0.85}$
α	$-{1.07}_{-0.93}^{+1.11}$ – $-{0.96}_{-0.90}^{+0.90}$	$-{1.15}_{-0.76}^{+0.76}$ – $-{1.25}_{-0.72}^{+0.80}$	$-{0.06}_{-0.88}^{+1.02}$ – $-{0.09}_{-0.98}^{+1.07}$	$-{0.75}_{-0.71}^{+0.70}$ – $-{0.82}_{-0.66}^{+0.76}$
β	$-{0.67}_{-0.35}^{+0.41}$ – $-{1.30}_{-0.44}^{+0.45}$	$-{0.90}_{-0.27}^{+0.28}$ – $-{1.24}_{-0.36}^{+0.33}$	$-{0.64}_{-0.36}^{+0.39}$ – $-{1.19}_{-0.42}^{+0.43}$	$-{0.89}_{-0.23}^{+0.26}$ – $-{1.27}_{-0.31}^{+0.32}$
γ	$+{3.09}_{-3.42}^{+3.64}$ – $+{5.68}_{-3.42}^{+2.60}$	$+{1.34}_{-2.04}^{+1.95}$ – $+{2.85}_{-2.24}^{+2.21}$	$+{2.72}_{-3.78}^{+3.19}$ – $+{5.03}_{-3.58}^{+2.90}$	$+{1.44}_{-1.64}^{+1.68}$ – $+{3.04}_{-1.97}^{+2.01}$
	Model 3
F0	$+{1.20}_{-0.57}^{+0.93}$ – $+{1.75}_{-0.77}^{+1.19}$	$+{1.57}_{-0.58}^{+0.93}$ – $+{1.90}_{-0.69}^{+1.08}$	$+{0.82}_{-0.31}^{+0.56}$ – $+{1.21}_{-0.49}^{+0.82}$	$+{1.25}_{-0.47}^{+0.69}$ – $+{1.62}_{-0.59}^{+0.81}$
α	$-{1.16}_{-0.88}^{+1.14}$ – $-{0.93}_{-0.81}^{+0.98}$	$-{1.54}_{-0.67}^{+0.68}$ – $-{1.32}_{-0.65}^{+0.71}$	$-{0.27}_{-0.98}^{+1.02}$ – $-{0.09}_{-0.99}^{+1.06}$	$-{1.14}_{-0.70}^{+0.76}$ – $-{0.95}_{-0.67}^{+0.70}$
β	$-{0.71}_{-0.36}^{+0.41}$ – $-{1.24}_{-0.41}^{+0.42}$	$-{0.97}_{-0.23}^{+0.27}$ – $-{1.33}_{-0.27}^{+0.26}$	$-{0.66}_{-0.34}^{+0.32}$ – $-{1.24}_{-0.39}^{+0.42}$	$-{0.96}_{-0.23}^{+0.25}$ – $-{1.35}_{-0.26}^{+0.27}$

Notes. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. "With Uncertainty" means the planet candidate radius, instellation flux, and host star effective temperature uncertainties are taken into account.

Download table as:  ASCIITypeset image

We compare models 1–3 using the Akaike information criterion (AIC), but the AIC does not indicate that one of the models is significantly more consistent with the data than another. The resulting relative likelihoods from the AIC analysis relative to model 1 are 0.31 for model 2 and 2.07 for model 3. Such low relative likelihood ratios are not considered compelling.

The F₀ parameter, giving the average number of planets per star in the solution domain (see Section 5), indicates that solutions using zero-completeness extrapolation (see Section 3.3.2) yield higher occurrence than those using constant-completeness extrapolation. This is because zero-completeness extrapolation induces larger completeness corrections. The zero-completeness extrapolation solution provides an upper bound on the HZ occurrence rate, and the constant extrapolation solution provides a lower bound. Reality will be somewhere in between and is likely to be closer to the zero extrapolation case for hotter stars, which have lower completeness in their HZs.

The hab stellar population is a subset of the hab2 population, so one may expect that they give similar solutions. But the hab stellar population contains significant regions of extrapolated completeness and low-reliability planet candidates for the flux range considered in our solution (see Figure 1). While the hab2 population contains the same regions, hab2 also contains many cooler stars that host higher-reliability planet candidates. These stars provide a better constraint on the power laws we use to describe the population. The result is that the hab2 solution has significantly smaller uncertainties than the hab solution, as seen in Tables 1 and 2.

Table 3 gives η_⊕, computed using the Poisson likelihood method, for the optimistic and conservative HZs for the hab and hab2 stellar populations and models 1–3. The low and high values correspond to the solutions using constant- and zero-completeness extrapolation, which bound the actual occurrence rates (see Section 3.3.2). We see the expected behavior of zero-completeness extrapolation leading to higher occurrence due to a larger completeness correction. Table 4 gives the same occurrence rates computed using the ABC method. The distributions of η_⊕ using these models and the Poisson likelihood method are shown in Figure 12. We see that for each model, when incorporating the input uncertainties, the hab and hab2 stellar populations yield consistent values of η_⊕. Without using the input uncertainties the hab population yields consistently lower values for η_⊕, though the difference is still within the 68% credible interval. Model 3 with hab2 gives generally higher occurrence rates than model 1. We also see that the median occurrence rates are about ∼10% higher when incorporating input uncertainties, qualitatively consistent with those of Shabram et al. (2020), who also see higher median occurrence rates when incorporating uncertainties. It is not clear what causes this increase in occurrence rates: on the one hand the sum of the inclusion probability, defined in Appendix C, for the planet candidates in Table 8 is 53.6, compared with the 54 planet candidates in the analysis without uncertainty, indicating that, on average, more planets exit the analysis than enter it when incorporating uncertainties. On the other hand, the sum of the inclusion probability times the planet radius is 106.6, compared with 105.0 for the planet candidates in the analysis without uncertainty, indicating that, on average, larger planets are entering the analysis. This may have an impact on the power-law model, leading to higher occurrence rates.

Zoom In Zoom Out Reset image size

Table 3.  η_⊕, Computed Using Population Models Based on the hab and hab2 Stellar Populations, with and without Uncertainties for Models 1–3 and Using the Poisson Method

	With Uncertainty		Without Uncertainty
	Based on hab Stars	Based on hab2 Stars	Based on hab Stars	Based on hab2 Stars
	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound
	Conservative HZ
Model 1	${0.30}_{-0.21}^{+0.69}$ – ${0.54}_{-0.39}^{+1.46}$	${0.37}_{-0.21}^{+0.48}$ – ${0.60}_{-0.36}^{+0.90}$	${0.15}_{-0.11}^{+0.32}$ – ${0.34}_{-0.25}^{+0.83}$	${0.34}_{-0.18}^{+0.37}$ – ${0.54}_{-0.30}^{+0.69}$
Model 2	${0.28}_{-0.20}^{+0.66}$ – ${0.53}_{-0.39}^{+1.46}$	${0.39}_{-0.23}^{+0.51}$ – ${0.56}_{-0.33}^{+0.83}$	${0.19}_{-0.13}^{+0.39}$ – ${0.32}_{-0.23}^{+0.78}$	${0.33}_{-0.18}^{+0.38}$ – ${0.53}_{-0.30}^{+0.70}$
Model 3	${0.31}_{-0.22}^{+0.69}$ – ${0.55}_{-0.39}^{+1.22}$	${0.59}_{-0.34}^{+0.74}$ – ${0.79}_{-0.44}^{+0.95}$	${0.21}_{-0.15}^{+0.42}$ – ${0.30}_{-0.21}^{+0.64}$	${0.50}_{-0.28}^{+0.60}$ – ${0.72}_{-0.39}^{+0.80}$
	Optimistic HZ
Model 1	${0.50}_{-0.35}^{+1.09}$ – ${0.80}_{-0.57}^{+2.07}$	${0.58}_{-0.33}^{+0.73}$ – ${0.88}_{-0.51}^{+1.27}$	${0.26}_{-0.18}^{+0.52}$ – ${0.51}_{-0.36}^{+1.17}$	${0.54}_{-0.28}^{+0.56}$ – ${0.80}_{-0.44}^{+0.99}$
Model 2	${0.48}_{-0.33}^{+1.06}$ – ${0.78}_{-0.58}^{+2.05}$	${0.61}_{-0.35}^{+0.77}$ – ${0.83}_{-0.48}^{+1.17}$	${0.32}_{-0.22}^{+0.62}$ – ${0.48}_{-0.33}^{+1.11}$	${0.52}_{-0.28}^{+0.58}$ – ${0.78}_{-0.44}^{+1.00}$
Model 3	${0.53}_{-0.37}^{+1.10}$ – ${0.81}_{-0.57}^{+1.73}$	${0.92}_{-0.52}^{+1.12}$ – ${1.14}_{-0.63}^{+1.35}$	${0.36}_{-0.25}^{+0.68}$ – ${0.46}_{-0.31}^{+0.92}$	${0.79}_{-0.44}^{+0.92}$ – ${1.04}_{-0.56}^{+1.14}$

Notes. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. "With Uncertainty" means the planet candidate radius, instellation flux, and host star effective temperature uncertainties are taken into account.

Download table as:  ASCIITypeset image

Table 4.  η_⊕, Computed Using Population Models Based on the hab and hab2 Stellar Populations, with and without Uncertainties for Models 1–3 and Using the ABC Method

	With Uncertainty		Without Uncertainty
	Based on hab Stars	Based on hab2 Stars	Based on hab Stars	Based on hab2 Stars
	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound
	Conservative HZ
Model 1	${0.33}_{-0.20}^{+0.46}$ – ${0.50}_{-0.31}^{+0.69}$	${0.40}_{-0.21}^{+0.45}$ – ${0.61}_{-0.33}^{+0.63}$	${0.16}_{-0.10}^{+0.24}$ – ${0.26}_{-0.16}^{+0.44}$	${0.29}_{-0.15}^{+0.28}$ – ${0.49}_{-0.27}^{+0.61}$
Model 2	${0.30}_{-0.19}^{+0.41}$ – ${0.52}_{-0.31}^{+0.63}$	${0.40}_{-0.21}^{+0.43}$ – ${0.64}_{-0.35}^{+0.73}$	${0.14}_{-0.09}^{+0.19}$ – ${0.25}_{-0.15}^{+0.36}$	${0.31}_{-0.15}^{+0.28}$ – ${0.47}_{-0.24}^{+0.50}$
Model 3	${0.35}_{-0.22}^{+0.47}$ – ${0.50}_{-0.29}^{+0.56}$	${0.69}_{-0.35}^{+0.64}$ – ${0.81}_{-0.40}^{+0.69}$	${0.18}_{-0.11}^{+0.27}$ – ${0.27}_{-0.16}^{+0.36}$	${0.50}_{-0.26}^{+0.48}$ – ${0.62}_{-0.31}^{+0.55}$
	Optimistic HZ
Model 1	${0.54}_{-0.33}^{+0.72}$ – ${0.73}_{-0.45}^{+0.98}$	${0.62}_{-0.32}^{+0.66}$ – ${0.89}_{-0.47}^{+0.89}$	${0.26}_{-0.16}^{+0.37}$ – ${0.39}_{-0.23}^{+0.62}$	${0.45}_{-0.23}^{+0.42}$ – ${0.71}_{-0.38}^{+0.86}$
Model 2	${0.48}_{-0.30}^{+0.66}$ – ${0.75}_{-0.44}^{+0.90}$	${0.62}_{-0.32}^{+0.63}$ – ${0.92}_{-0.49}^{+1.02}$	${0.24}_{-0.14}^{+0.30}$ – ${0.37}_{-0.22}^{+0.52}$	${0.47}_{-0.23}^{+0.41}$ – ${0.68}_{-0.34}^{+0.70}$
Model 3	${0.56}_{-0.36}^{+0.72}$ – ${0.73}_{-0.41}^{+0.81}$	${1.05}_{-0.52}^{+0.96}$ – ${1.15}_{-0.57}^{+0.99}$	${0.30}_{-0.18}^{+0.43}$ – ${0.39}_{-0.23}^{+0.51}$	${0.75}_{-0.39}^{+0.72}$ – ${0.89}_{-0.44}^{+0.78}$

Notes. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. "With Uncertainty" means the planet candidate radius, instellation flux, and host star effective temperature uncertainties are taken into account.

Download table as:  ASCIITypeset image

Table 5 gives the occurrence rates for a variety of planet radius and host star effective temperature ranges, computed using the hab2 stellar population and models 1–3. We see that the uncertainties for the 1.5–2.5 R_⊕ planets are significantly smaller than those for the 0.5–1.5 R_⊕ planets, indicating that the large uncertainties in η_⊕ are due to the small number of observed planets in the 0.5–1.5 R_⊕ range. The distributions of occurrence for the two bounding extrapolation types are shown in Figure 13. The difference between these two bounding cases is smaller than the uncertainties. Table 6 gives the 95% and 99% intervals for the 0.5–1.5 R_⊕ planets using model 1 computed with hab2.

Zoom In Zoom Out Reset image size

Table 5.  Number of Habitable Zone Planets per Star for Various Ranges of Planet Radii and Host Star Effective Temperatures, Computed Using the Population Model Based on the hab2 Stellar Population with the Poisson Likelihood Method and Incorporating Uncertainties in Planet Radius, Instellation Flux, and Host Star Effective Temperature

Planet Radius	4800–6300 K	3900–6300 K	3900–5300 K (K Dwarfs)	5300–6000 K (G Dwarfs)
	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound
	Conservative HZ
	Model 1
0.5–1.5 R⊕	${\bf{0}}.{{\bf{37}}}_{-{\bf{0.21}}}^{+{\bf{0.48}}}$ – ${\bf{0}}.{{\bf{60}}}_{-{\bf{0.36}}}^{+{\bf{0.90}}}$	${0.35}_{-0.19}^{+0.43}$ – ${0.50}_{-0.29}^{+0.73}$	${0.32}_{-0.17}^{+0.35}$ – ${0.42}_{-0.23}^{+0.50}$	${0.38}_{-0.22}^{+0.50}$ – ${0.63}_{-0.38}^{+0.94}$
1.5–2.5R⊕	${0.16}_{-0.05}^{+0.07}$ – ${0.24}_{-0.08}^{+0.14}$	${0.15}_{-0.05}^{+0.06}$ – ${0.20}_{-0.07}^{+0.12}$	${0.14}_{-0.04}^{+0.05}$ – ${0.17}_{-0.06}^{+0.07}$	${0.17}_{-0.05}^{+0.06}$ – ${0.26}_{-0.09}^{+0.13}$
0.5–2.5R⊕	${0.54}_{-0.24}^{+0.52}$ – ${0.85}_{-0.42}^{+0.99}$	${0.51}_{-0.22}^{+0.46}$ – ${0.71}_{-0.34}^{+0.80}$	${0.46}_{-0.19}^{+0.37}$ – ${0.60}_{-0.26}^{+0.52}$	${0.56}_{-0.25}^{+0.53}$ – ${0.90}_{-0.45}^{+1.01}$
	Model 2
0.5–1.5R⊕	${\bf{0}}.{{\bf{39}}}_{-{\bf{0.23}}}^{+{\bf{0.51}}}$ – ${\bf{0}}.{{\bf{56}}}_{-{\bf{0.33}}}^{+{\bf{0.83}}}$	${0.36}_{-0.20}^{+0.44}$ – ${0.47}_{-0.27}^{+0.67}$	${0.33}_{-0.18}^{+0.36}$ – ${0.40}_{-0.22}^{+0.47}$	${0.40}_{-0.23}^{+0.53}$ – ${0.59}_{-0.35}^{+0.86}$
1.5–2.5R⊕	${0.16}_{-0.05}^{+0.06}$ – ${0.24}_{-0.08}^{+0.13}$	${0.15}_{-0.05}^{+0.06}$ – ${0.20}_{-0.07}^{+0.11}$	${0.14}_{-0.04}^{+0.05}$ – ${0.17}_{-0.06}^{+0.07}$	${0.17}_{-0.05}^{+0.06}$ – ${0.25}_{-0.08}^{+0.12}$
0.5–2.5R⊕	${0.56}_{-0.26}^{+0.54}$ – ${0.81}_{-0.39}^{+0.91}$	${0.51}_{-0.23}^{+0.47}$ – ${0.68}_{-0.32}^{+0.75}$	${0.47}_{-0.20}^{+0.38}$ – ${0.58}_{-0.25}^{+0.50}$	${0.57}_{-0.27}^{+0.56}$ – ${0.85}_{-0.42}^{+0.94}$
	Model 3
0.5–1.5R⊕	${\bf{0}}.{{\bf{59}}}_{-{\bf{0.34}}}^{+{\bf{0.74}}}$ – ${\bf{0}}.{{\bf{79}}}_{-{\bf{0.44}}}^{+{\bf{0.95}}}$	${0.43}_{-0.26}^{+0.63}$ – ${0.59}_{-0.35}^{+0.82}$	${0.31}_{-0.17}^{+0.37}$ – ${0.43}_{-0.24}^{+0.51}$	${0.64}_{-0.35}^{+0.72}$ – ${0.85}_{-0.46}^{+0.93}$
1.5–2.5R⊕	${0.20}_{-0.07}^{+0.11}$ – ${0.29}_{-0.10}^{+0.15}$	${0.15}_{-0.06}^{+0.13}$ – ${0.22}_{-0.09}^{+0.17}$	${0.11}_{-0.04}^{+0.05}$ – ${0.16}_{-0.05}^{+0.07}$	${0.22}_{-0.06}^{+0.07}$ – ${0.32}_{-0.09}^{+0.11}$
0.5–2.5R⊕	${0.81}_{-0.39}^{+0.79}$ – ${1.11}_{-0.52}^{+1.02}$	${0.60}_{-0.32}^{+0.71}$ – ${0.84}_{-0.44}^{+0.92}$	${0.42}_{-0.20}^{+0.39}$ – ${0.60}_{-0.28}^{+0.54}$	${0.87}_{-0.38}^{+0.75}$ – ${1.18}_{-0.51}^{+0.97}$
	Optimistic HZ
	Model 1
0.5–1.5R⊕	${\bf{0}}.{{\bf{58}}}_{-{\bf{0.33}}}^{+{\bf{0.73}}}$ – ${\bf{0}}.{{\bf{88}}}_{-{\bf{0.51}}}^{+{\bf{1.28}}}$	${0.54}_{-0.29}^{+0.64}$ – ${0.73}_{-0.41}^{+1.02}$	${0.49}_{-0.26}^{+0.51}$ – ${0.60}_{-0.33}^{+0.69}$	${0.60}_{-0.34}^{+0.75}$ – ${0.93}_{-0.55}^{+1.32}$
1.5–2.5R⊕	${0.25}_{-0.06}^{+0.09}$ – ${0.35}_{-0.11}^{+0.19}$	${0.23}_{-0.06}^{+0.09}$ – ${0.29}_{-0.10}^{+0.17}$	${0.21}_{-0.06}^{+0.07}$ – ${0.25}_{-0.08}^{+0.09}$	${0.26}_{-0.07}^{+0.09}$ – ${0.37}_{-0.11}^{+0.17}$
0.5–2.5R⊕	${0.84}_{-0.37}^{+0.78}$ – ${1.25}_{-0.59}^{+1.39}$	${0.78}_{-0.33}^{+0.68}$ – ${1.04}_{-0.47}^{+1.13}$	${0.71}_{-0.28}^{+0.53}$ – ${0.86}_{-0.36}^{+0.72}$	${0.87}_{-0.38}^{+0.79}$ – ${1.33}_{-0.63}^{+1.42}$
	Model 2
0.5–1.5R⊕	${\bf{0}}.{{\bf{61}}}_{-{\bf{0.35}}}^{+{\bf{0.77}}}$ – ${\bf{0}}.{{\bf{83}}}_{-{\bf{0.48}}}^{+{\bf{1.17}}}$	${0.55}_{-0.31}^{+0.66}$ – ${0.68}_{-0.39}^{+0.95}$	${0.50}_{-0.27}^{+0.54}$ – ${0.57}_{-0.31}^{+0.66}$	${0.63}_{-0.36}^{+0.79}$ – ${0.87}_{-0.51}^{+1.22}$
1.5–2.5R⊕	${0.25}_{-0.06}^{+0.09}$ – ${0.34}_{-0.10}^{+0.17}$	${0.23}_{-0.07}^{+0.09}$ – ${0.29}_{-0.10}^{+0.16}$	${0.21}_{-0.06}^{+0.07}$ – ${0.25}_{-0.08}^{+0.09}$	${0.26}_{-0.07}^{+0.08}$ – ${0.37}_{-0.11}^{+0.16}$
0.5–2.5R⊕	${0.87}_{-0.39}^{+0.81}$ – ${1.19}_{-0.56}^{+1.28}$	${0.79}_{-0.34}^{+0.70}$ – ${0.99}_{-0.45}^{+1.05}$	${0.72}_{-0.29}^{+0.55}$ – ${0.84}_{-0.35}^{+0.69}$	${0.90}_{-0.40}^{+0.83}$ – ${1.25}_{-0.59}^{+1.32}$
	Model 3
0.5–1.5R⊕	${\bf{0}}.{{\bf{92}}}_{-{\bf{0.52}}}^{+{\bf{1.12}}}$ – ${\bf{1}}.{{\bf{14}}}_{-{\bf{0.63}}}^{+{\bf{1.35}}}$	${0.67}_{-0.40}^{+0.96}$ – ${0.85}_{-0.50}^{+1.17}$	${0.46}_{-0.26}^{+0.55}$ – ${0.61}_{-0.33}^{+0.70}$	${1.00}_{-0.54}^{+1.08}$ – ${1.23}_{-0.65}^{+1.32}$
1.5–2.5R⊕	${0.31}_{-0.11}^{+0.16}$ – ${0.42}_{-0.15}^{+0.21}$	${0.22}_{-0.10}^{+0.20}$ – ${0.31}_{-0.13}^{+0.25}$	${0.16}_{-0.05}^{+0.07}$ – ${0.23}_{-0.07}^{+0.10}$	${0.34}_{-0.08}^{+0.10}$ – ${0.46}_{-0.12}^{+0.14}$
0.5–2.5R⊕	${1.26}_{-0.59}^{+1.19}$ – ${1.60}_{-0.73}^{+1.43}$	${0.92}_{-0.49}^{+1.07}$ – ${1.20}_{-0.62}^{+1.30}$	${0.63}_{-0.29}^{+0.58}$ – ${0.85}_{-0.38}^{+0.75}$	${1.35}_{-0.57}^{+1.11}$ – ${1.71}_{-0.71}^{+1.36}$

Notes. The η_⊕ values for model 1 are shown in boldface. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. As explained in Section 5.2 we recommend model 1 as the baseline model. Results from the other models are included for comparison.

Download table as:  ASCIITypeset image

Table 6.  Credible Intervals of the Upper and Lower Bounds on HZ Occurrence for Model 1 Computed Using the Population Model Based on the hab2 Stellar Population (see Section 1) and Accounting for Input Uncertainty

	Low	High	Total
95% Credible Interval
${\eta }_{\oplus }^{{\rm{C}}}$	[0.07, 1.91]	[0.10, 3.77]	[0.07, 3.77]
${\eta }_{\oplus }^{{\rm{O}}}$	[0.11, 2.88]	[0.16, 5.29]	[0.11, 5.29]
${\eta }_{\oplus ,{\rm{G}}}^{{\rm{C}}}$	[0.07, 1.92]	[0.10, 3.76]	[0.07, 3.76]
${\eta }_{\oplus ,{\rm{G}}}^{{\rm{O}}}$	[0.11, 2.90]	[0.16, 5.26]	[0.11, 5.26]
${\eta }_{\oplus ,{\rm{K}}}^{{\rm{C}}}$	[0.07, 1.34]	[0.09, 1.92]	[0.07, 1.92]
${\eta }_{\oplus ,{\rm{K}}}^{{\rm{O}}}$	[0.11, 1.96]	[0.13, 2.66]	[0.11, 2.66]
99% Credible Interval
${\eta }_{\oplus }^{{\rm{C}}}$	[0.04, 3.19]	[0.06, 6.91]	[0.04, 6.91]
${\eta }_{\oplus }^{{\rm{O}}}$	[0.06, 4.76]	[0.09, 9.58]	[0.06, 9.58]
${\eta }_{\oplus ,{\rm{G}}}^{{\rm{C}}}$	[0.04, 3.13]	[0.06, 6.57]	[0.04, 6.57]
${\eta }_{\oplus ,{\rm{G}}}^{{\rm{O}}}$	[0.06, 4.65]	[0.10, 9.09]	[0.06, 9.09]
${\eta }_{\oplus ,{\rm{K}}}^{{\rm{C}}}$	[0.04, 2.06]	[0.05, 3.06]	[0.04, 3.06]
${\eta }_{\oplus ,{\rm{K}}}^{{\rm{O}}}$	[0.07, 2.97]	[0.08, 4.20]	[0.07, 4.20]

Note. See Section 1.3 for the definitions of the different types of η_⊕.

Download table as:  ASCIITypeset image

Figure 14 shows the dependence of the HZ occurrence rate on effective temperature for models 1–3 based on the hab2 stellar population and for model 1 for the hab stellar population. For each model, the occurrence using zero and constant extrapolation is shown. Models 1 and 2 show a weak increase in occurrence with increasing effective temperature. Model 3 shows a stronger increase in occurrence with effective temperature, consistent with model 3's assumption that the only temperature dependence is the geometric effect described in Section 2.2. However, as shown in Figure 6, the difference between constant and zero extrapolated completeness is near zero for T_eff ≤ 4500 K, so we would expect the difference in occurrence rates to be close to zero in that temperature range. This is true for models 1 and 2 but not true for model 3. We take this as evidence that models 1 and 2 are correctly measuring a T_eff dependence beyond the geometric effect. We recognize that the statistical evidence for this T_eff dependence is not compelling, since the overlapping 68% credible intervals for the two completeness extrapolations would allow an occurrence rate independent of T_eff.

Zoom In Zoom Out Reset image size

An issue that arises with zero-completeness extrapolation (= high bound) is that, strictly speaking, planet candidates with orbital periods of >500 days are in a region of zero completeness around their host star and would not contribute to the Poisson likelihood (see Equation (16) in Appendix B). There is one such planet candidate in the hab2 stellar sample with reliability = 0.67. Performing our Poisson inference, we remove this planet with a period of >500 days for model 1 and incorporate input uncertainties, yielding an optimistic η_⊕ = ${0.70}_{-0.41}^{+1.01}$, compared with ${0.88}_{-0.51}^{+1.27}$ (from Table 3) when including the planet. While this result is well within the 68% credible interval of the result with the planet included, removing this planet has a noticeable impact on the upper bound for the optimistic η_⊕. However this planet is in fact detected, implying that completeness is not zero for periods of >500 days, at least for this planet's host star. If the actual completeness is very close to zero, a planet detection implies a large population. We therefore leave this planet in the analysis, thinking of "zero completeness" as a limit of the completeness going to zero when the HZ includes orbital periods of >500 days, summed or averaged over the stellar population for our computations.

We choose to compute our occurrence rates as a function of instellation flux for two major reasons: this allows a more direct characterization of each star's HZ, and it allows a more constrained extrapolation to longer orbital periods than working directly in terms of orbital period.

By considering instellation flux, we can measure HZ occurrence by including observed planets from the HZ of their host stars, which is not possible across a wide range of stellar temperatures when using orbital period (see Figure 2). Instellation flux also allows a direct measurement of the impact of uncertainties in stellar effective temperature and planetary instellation flux.

The HZ of most G and F stars includes orbital periods that are beyond those periods well covered by Kepler observations (see Figures 1 and 2), requiring significant extrapolation of orbital period–based planet population models to long orbital periods. Such extrapolation is poorly known or constrained, leading to possible significant and unbounded inaccuracies. In instellation flux, however, there is planet data throughout those regions of our domain of analysis that have reasonable completeness (see Figure 1), so no extrapolation in instellation flux is required. In this sense, replacing orbital period with instellation flux (determined by the orbital period for each star) moves the problem of extrapolating the population model to longer orbital periods to extrapolating the completeness data to lower instellation flux. In Section 3.3.2 we argue that completeness, on average, decreases monotonically with decreasing instellation flux. This allows us to bound the extrapolated completeness between no decrease at all (constant extrapolation) and zero completeness for instellation flux for orbital periods beyond the 500 day limit where completeness is measured. We then perform our analysis for the two extrapolation cases and find that their difference in HZ occurrence rates is small relative to our uncertainties. In this way we provide a bounded estimate of HZ occurrence rates using instellation flux, rather than the unbounded extrapolation resulting from using orbital period.

Our approach to measuring η_⊕ is to compute the planet population rate model λ(r, I, T_eff) ≡ d²f(r, I, T_eff)/dr dI, integrate over r and I, and average over T_eff. We compute the population model λ using the hab and hab2 stellar populations (Section 3.1) to measure the sensitivity of our results to stellar type, and we consider several possible functional forms for λ (Equation (5)).

5.2.1. Comparing Population Rate Function Models

5.2.2. Comparing the hab and hab2 Stellar Populations

We find reasonable consistency in η_⊕ across the models for both the hab and hab2 stellar populations as shown in Tables 3 and 4. Table 5 gives the occurrence rates for several planet radius and stellar effective temperature ranges, using the hab2 population model. The uncertainties reflecting our 68% credible intervals for our η_⊕, counting HZ planets with radii of 0.5–1.5 R_⊕, are large, with positive uncertainties being nearly 150% of the value when using input uncertainties. Comparing occurrence rate uncertainties with and without input uncertainties in Tables 3 and 4, we see that the bulk of the uncertainties occur without input uncertainties, while using input uncertainties increases the output uncertainty by nearly 20%. The much smaller uncertainties for larger planets (1.5 R_⊕ ≤ r ≤ 2.5 R_⊕) in Table 5 suggest the large uncertainties in η_⊕ are driven by the very small number of detections in the η_⊕ range, combined with the very low completeness (see Figure 1). Low completeness will cause large completeness corrections, which will magnify the Poisson uncertainty due to few planet detections. There are more planets with larger radii that are in a regime of higher completeness, resulting in lower uncertainties for larger planet occurrence rates.

Estimates of η_⊕ are useful in calculating exoEarth yields from direct-imaging missions, such as in the flagship concept studies for LUVOIR and HabEX. These mission studies assumed an occurrence rate of ${\eta }_{\oplus }={0.24}_{-0.16}^{+0.46}$ for Sun-like stars based on the NASA ExoPAG SAG13 meta-analysis of Kepler data (Kopparapu et al. 2018). The expected exoEarth candidate yields from the mission study reports are ${54}_{-34}^{+61}$ for LUVOIR-A (15 m), ${28}_{-17}^{+30}$ for LUVOIR-B (8 m), and 8 for the HabEx 4 m baseline configuration, which combines a more traditional coronagraphic starlight suppression system with a formation-flying starshade occulter. Table 5 provides the η_⊕ values based on the three models for G (Sun-like) and K dwarfs. If we assume the range of η_⊕,G from the conservative to the optimistic HZ from Table 5 for planets of 0.5–1.5 R_⊕ from the "low" end, say, for model 1, η_⊕,G would be between ${0.38}_{-0.22}^{+0.50}$ and ${0.60}_{-0.34}^{+0.75}$. While these η_⊕ values appear to be larger than the ${0.24}_{-0.16}^{+0.46}$ occurrence rate assumed by the mission studies, it should be noted that these studies adopted radius range of 0.82 to 1.4 R_⊕, and a lower limit of $0.8\,\times \,{a}^{-0.5}$, where a is the HZ-corrected semi-major axis. This is slightly a smaller HZ region and lower radius than the one used in our study. As a result, it is possible that we might be agreeing with their assumed η⊕ value if we use the same bounding boxes. Computing the conservative HZ as described in Section 3.5 but replacing the planet radius range with $0.82\le r\le 1.4\,{R}_{\oplus }$ gives a lower bound of ${0.18}_{-0.09}^{+0.16}$ and an upper bound of ${0.28}_{-0.14}^{+0.30}$, nicely bracketing the value assumed in mission studies.

η_⊕ can also be used to estimate, on average, the nearest HZ planet around a G and K dwarf assuming the planets are distributed randomly. Within the solar neighborhood, the stellar number density ranges from 0.0033 to 0.0038 pc⁻³ for G dwarfs and 0.0105 to 0.0153 pc⁻³ for K dwarfs (Mamajek & Hillenbrand 2008; Kirkpatrick et al. 2012).⁵³ For G dwarfs, multiplying with the conservative (${0.38}_{-0.22}^{+0.50}$) model 1 "low" end of the η_⊕,G values (i.e., the number of planets per star), we get between ${0.0013}_{-0.0007}^{+0.0016}$ and ${0.0014}_{-0.0008}^{+0.0019}$ HZ planets pc⁻³. The nearest HZ planet around a G dwarf would then be expected to be at a distance of d = (3/(4π × n_p))^1/3, where n_p is the planet number density in pc⁻³. Substituting, we get a d between ${5.9}_{-1.26}^{+1.76}$ pc and ${5.5}_{-1.34}^{+1.83}$ pc, or essentially around ∼6 pc away. A similar calculation for K dwarfs assuming the model 1 conservative HZ η_⊕,K values from Table 5 indicates that, on average, the nearest HZ planet could be between ${4.1}_{-0.90}^{+1.19}$ pc and ${3.6}_{-0.79}^{+1.05}$ pc, or around ∼4 pc away.

An additional speculative calculation one could do is to take the number of G dwarfs in the solar neighborhood within 10 pc—19 from RECONS⁵⁴ —and multiply it with the "low" conservative η_⊕,G value from model 1, ${0.38}_{-0.22}^{+0.50}$. We then get ${7.2}_{-4.2}^{+9.5}$ HZ planets around G dwarfs (of all sub–spectral types) within 10 pc. A similar calculation for K dwarfs from the same RECONS data with 44 stars and a model 1 "low" value, conservative HZ ${\eta }_{\oplus ,{\rm{K}}}={0.32}_{-0.17}^{+0.35}$ indicates that there are ${14}_{-7.5}^{+15}$ HZ planets around K dwarfs within 10 pc. It should be noted that the number for the nearest HZ planet and the number of HZ planets in the solar neighborhood use the "low" end of the rocky planet (0.5–1.5 R_⊕) occurrence rate values from Table 5. As such, these represent the lower-bound estimates. In other words, there may potentially be an HZ rocky planet around a G or K dwarf that is closer, and maybe more HZ planets, than the values quoted above.

This can be quantified from the numbers shown in Table 6, which provides the 95% and 99% credible intervals of the upper and lower bounds on HZ occurrence for model 1 computed with hab2 and accounting for input uncertainty. If we use only the "low" values and the lower end of the conservative HZ occurrence values from this table (0.07 for 95% credible interval, 0.04 for 99%), then the nearest HZ planet around a G or K dwarf star is within ∼6 pc away with 95% confidence and within ∼7.5 pc away with 99% confidence. Similarly, there could be ∼4 HZ planets within 10 pc with 95% confidence and ∼3 HZ planets with 99% confidence. The lower bound of the 95% credible interval for G dwarfs is also 0.07, so assuming 100 billion stars in our Galaxy and a G-dwarf fraction of 0.041, there are likely at least 287 million rocky conservative HZ planets orbiting G-dwarf stars in our Galaxy.

The numbers provided in this section are first-order estimates to simply show a meaningful application of η_⊕, given its uncertainties. Weiss et al. (2018), Mulders et al. (2018), and He et al. (2020) have shown that there is strong evidence for multiplicity and "clustering'' (planets in multiplanet systems are correlated in size to their neighbors). However, the similar-size trend has not been characterized for orbits corresponding to the HZs of G and K dwarfs. If the pattern of similarly sized planets extends to HZ, the nearest sun-like star hosting an Earth-size planet in the HZ would be farther from our sun than we estimated in this paper, but it could have multiple Earth-sized planets in the HZ.

Our work is the first to compute HZ occurrence rates using the incident stellar flux for HZs around subsets of FGK stars based on the DR25 catalog and Gaia-based stellar properties and using the HZ definition of Kopparapu et al. (2014). Other works in the literature have produced occurrence rates for orbital periods related to FGK HZs, but as discussed in Section 1.2 and Section 3.5, we measure occurrence between HZ boundaries as a function of T_eff. This is a different quantity from occurrence rates based on orbital period. The few occurrence rate estimates in the literature based on instellation flux, such as that of Petigura et al. (2013), use a rectangular region in terms of radius and instellation flux and only approximate the HZ. Therefore, we do not directly compare our occurrence rates with those in previous literature.

We provide a formal computation of Γ_⊕, which is commonly used to compare η_⊕ estimates (see, for example, Figure 14 of Kunimoto & Matthews 2020). For a planet of period p and radius r, we define ${\rm{\Gamma }}\equiv {d}^{2}f/d\mathrm{log}p\,d\mathrm{log}r$ = pr d²f/dp dr, and Γ_⊕ is Γ evaluated at Earth's period and radius. We need to express Γ in terms of instellation flux I, which we will do using ${d}^{2}f/{dp}\,{dr}=\left({d}^{2}f/{dI}\,{dr}\right)\left({dI}/{dp}\right)$.

For a particular star, the instellation on a planet at period p in years is given by $I={R}_{* }^{2}{T}^{4}{M}_{* }^{-2/3}{p}^{-4/3}$, where R_* is the stellar radius in solar radii, M_* is the stellar mass in solar masses, and T = T_eff/T_⊙ is the star's effective temperature divided by the solar effective temperature. Then

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }} & = & p\,r\,\displaystyle \frac{{d}^{2}f}{{dI}\,{dr}}\left(\displaystyle \frac{{dI}}{{dp}}\right)\\ & = & \ -\displaystyle \frac{4}{3}\,\displaystyle \frac{{R}_{* }^{2}{T}^{4}}{{M}_{* }^{\tfrac{2}{3}}}r\,{p}^{-\tfrac{4}{3}}\,\lambda (I,r,T,{\boldsymbol{\theta }})\end{array}\end{eqnarray} \tag{ 10 }$$

because d²f/dI dr = λ(I, r, T, θ), one of the differential population rate models from Equation (5). To compute Γ_⊕ for a particular star, we evaluate Equation (10) at r = 1  R_⊕, p = 1 year, and $I={R}_{* }^{2}{T}^{4}{M}_{* }^{-2/3}$, the instellation a planet with a 1 year orbital period would have from that star. The result is the Γ_⊕ in radius and period implied by our differential population rate function in radius and instellation for that star and may be compared directly with Γ_⊕ from period-based occurrence studies.

We compute Γ_⊕ using model 1 from Equation (5) with input uncertainty on the hab2 stellar population. For each star in hab2, we evaluate Equation (10) using the posterior θ distribution and concatenate the resulting Γ_⊕ distributions from all the stars. We do this for both completeness extrapolations in Section 3.3.2, giving low and high bounds. This results in a Γ_⊕ between ${0.45}_{-0.24}^{+0.46}$ and ${0.50}_{-0.26}^{+0.46}$. While this is a formal mathematical exercise that has not been demonstrated to be truly equivalent to Γ_⊕ defined in period space, the match between this value and our conservative HZ ${\eta }_{\oplus }^{{\rm{C}}}={0.37}_{-0.21}^{+0.48}$ – ${0.60}_{-0.36}^{+0.90}$ for model 1 for hab2 with input uncertainty in Table 3 is remarkable.

Our value of Γ_⊕ is somewhat higher than values using post-Gaia stellar and planetary data (see, for example, Figure 14 of Kunimoto & Matthews 2020), but not significantly so. For example, Bryson et al. (2020a) found ${{\rm{\Gamma }}}_{\oplus }={0.09}_{-0.04}^{+0.07}$ when correcting for reliability in the period–radius space. Using the same population model, Bryson et al. (2020a) found a SAG13 η_⊕ value of ${0.13}_{-0.06}^{+0.10}$. It is not clear by how much we should expect Γ_⊕ to correspond to η_⊕.

As described in Section 4.2 and Figure 14, our results indicate a weak, and thus not compelling, increase in HZ planet occurrence with increasing stellar effective temperature. This T_eff dependence is weaker than would be expected if planet occurrence were uniformly spaced in the semimajor axis (see Section 2.2) because hotter stars have larger HZs. This can be seen quantitatively in the median differential population rates for models 1 and 2 using the hab2 population in Tables 1 and 2. In model 2 we observe a median T_eff exponent γ < 3, compared with the prediction of γ ≈ 3 to 4.5 due to the larger HZ for hotter stars from Equation (4). This is reflected in model 1, which includes the correction for the larger HZ so if T_eff dependence were due only to the larger HZ, then γ would equal 0. The high bound of model 1 finds a median γ < −1, indicating that we see fewer HZ planets than expected in the larger HZ of hotter stars if the planets were uniformly spaced. However the upper limits of γ's 68% credible interval in Tables 1 and 2 are consistent with the prediction of uniform planet spacing in larger HZs for hotter stars. For example, the posterior of γ for model 1 (hab2 population, high bound) has γ ≥ 0 22.3% of the time.

Our detection of a weaker T_eff dependence than expected from larger HZs is qualitatively consistent with the increasing planet occurrence with lower T_eff found in Garrett et al. (2018) and Mulders et al. (2015). But our uncertainties do not allow us to make definitive conclusions about this T_eff dependence.

To study the dependence of our result on the planet sample, we perform a bootstrap analysis. We run a Poisson likelihood inference using hab2 and model 1 with zero extrapolation (high bound) 400 times, resampling the planet candidate population with replacement. Each resampled run removes planets according to their reliability as described in Section 3.4.2 but does not consider input uncertainty. The concatenated posterior of these resampled runs gives ${F}_{0}={1.404}_{-0.680}^{+1.768}$, $\alpha =-{0.920}_{-1.072}^{+1.236}$, $\beta =-{1.175}_{-0.444}^{+0.465}$, and $\gamma =-{1.090}_{-2.217}^{+2.446}$. These parameters yield ${\eta }_{\oplus }^{{\rm{C}}}={0.483}_{-0.324}^{+0.997}$ and ${\eta }_{\oplus }^{{\rm{O}}}={0.716}_{-0.472}^{+1.413}$. Comparing these with the hab2 model 1 high value without uncertainty in Tables 1 and 3, we see that the central values from the bootstrap study are well within the 68% credible interval of our results and the uncertainties are as much as 50% higher.

A similar study of dependence on the stellar sample is not feasible because each resampled stellar population would require a full recomputation of detection and vetting completeness and reliability. Performing hundreds of these computations is beyond our available resources.

All results presented in this paper are computed with corrections for planet catalog completeness and reliability (see Section 3.3). Figure 15 shows an example of what happens when there is no correction for catalog reliability. We compute ${\eta }_{\oplus }^{{\rm{C}}}$, the occurrence in the conservative HZ, with model 1, zero-completeness extrapolation (high value), accounting for input uncertainty and using the hab2 stellar population. With reliability correction, we have ${\eta }_{\oplus }^{{\rm{C}}}={0.60}_{-0.36}^{+0.90}$, and without reliability correction we have ${\eta }_{\oplus }^{{\rm{C}}}={1.25}_{-0.60}^{+1.40}$. In this typical case reliability has a factor-of-two impact, consistent with Bryson et al. (2020a), though because of the large uncertainties the difference is less than the 68% credible interval.

Zoom In Zoom Out Reset image size

Our definition of η_⊕, restricted to stars with effective temperatures between 4800 K and 6300 K, varies somewhat from the literature. To connect with other occurrence rate studies we repeat our analysis using the GK (3900 K ≤ T_eff ≤ 6000 K) and FGK (3900 K ≤ T_eff ≤ 7300 K) stellar populations to compute planet population models, with the results listed in Table 7. We provide our η_⊕ derived from these stellar populations as well as the HZ occurrence for the GK and FGK T_eff ranges. The values for our definition of η_⊕ are consistent with the values in Tables 3 and 4. We caution, however, that the FGK population extends well into stellar effective temperatures where there are no planet detections and very low or zero completeness, so an FGK result necessarily involves extrapolation from cooler stars.

Table 7.  Parameter Fits and η_⊕ with 68% Confidence Limits for Model 1 from Equation (5) Computed Using the Population Model from the Poisson Likelihood Method Applied to the GK and FGK Stellar Populations

	With Uncertainty		Without Uncertainty
	Based on GK Stars	Based on FGK Stars	Based on GK Stars	Based on FGK Stars
	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound	Low Bound–High Bound
	Model 1
F0	$+{1.01}_{-0.4}^{+0.82}$ – $+{1.23}_{-0.52}^{+1.12}$	$+{1.26}_{-0.55}^{+1.2}$ – $+{2.36}_{-1.2}^{+2.94}$	$+{0.89}_{-0.33}^{+0.64}$ – $+{1.12}_{-0.44}^{+0.83}$	$+{1.24}_{-0.51}^{+1.03}$ – $+{2.05}_{-0.98}^{+2.14}$
α	$-{1.00}_{-0.9}^{+1.0}$ – $-{1.05}_{-0.9}^{+1.0}$	$-{1.03}_{-0.87}^{+0.95}$ – $-{1.16}_{-0.85}^{+0.94}$	$-{0.80}_{-0.83}^{+0.93}$ – $-{0.90}_{-0.78}^{+0.87}$	$-{1.00}_{-0.77}^{+0.84}$ – $-{0.96}_{-0.78}^{+0.87}$
β	$-{0.90}_{-0.32}^{+0.34}$ – $-{1.13}_{-0.36}^{+0.37}$	$-{0.82}_{-0.3}^{+0.32}$ – $-{1.23}_{-0.34}^{+0.36}$	$-{0.84}_{-0.3}^{+0.33}$ – $-{1.08}_{-0.34}^{+0.36}$	$-{0.82}_{-0.27}^{+0.3}$ – $-{1.21}_{-0.32}^{+0.34}$
γ	$-{3.00}_{-1.72}^{+1.75}$ – $-{2.25}_{-1.86}^{+1.9}$	$-{2.67}_{-1.58}^{+1.57}$ – $-{1.07}_{-1.77}^{+1.8}$	$-{2.76}_{-1.62}^{+1.61}$ – $-{1.98}_{-1.69}^{+1.75}$	$-{2.67}_{-1.43}^{+1.43}$ – $-{1.15}_{-1.66}^{+1.64}$
${\eta }_{\oplus }^{{\rm{C}}}$	${0.34}_{-0.20}^{+0.48}$ – ${0.46}_{-0.28}^{+0.71}$	${0.36}_{-0.20}^{+0.46}$ – ${0.63}_{-0.37}^{+0.92}$	${0.29}_{-0.16}^{+0.37}$ – ${0.41}_{-0.23}^{+0.53}$	${0.35}_{-0.19}^{+0.40}$ – ${0.53}_{-0.30}^{+0.68}$
${\eta }_{\oplus ,\mathrm{GK}}^{{\rm{C}}}$	${0.32}_{-0.18}^{+0.41}$ – ${0.40}_{-0.23}^{+0.54}$	${0.32}_{-0.18}^{+0.38}$ – ${0.48}_{-0.27}^{+0.64}$	${0.26}_{-0.15}^{+0.31}$ – ${0.35}_{-0.19}^{+0.40}$	${0.32}_{-0.16}^{+0.33}$ – ${0.41}_{-0.22}^{+0.49}$
${\eta }_{\oplus ,\mathrm{FGK}}^{{\rm{C}}}$	${0.35}_{-0.21}^{+0.54}$ – ${0.47}_{-0.29}^{+0.81}$	${0.37}_{-0.21}^{+0.53}$ – ${0.63}_{-0.39}^{+1.18}$	${0.30}_{-0.17}^{+0.43}$ – ${0.42}_{-0.24}^{+0.64}$	${0.36}_{-0.20}^{+0.46}$ – ${0.53}_{-0.31}^{+0.88}$
${\eta }_{\oplus }^{{\rm{O}}}$	${0.52}_{-0.30}^{+0.72}$ – ${0.68}_{-0.41}^{+1.01}$	${0.56}_{-0.31}^{+0.70}$ – ${0.92}_{-0.54}^{+1.29}$	${0.45}_{-0.25}^{+0.56}$ – ${0.61}_{-0.34}^{+0.77}$	${0.55}_{-0.29}^{+0.60}$ – ${0.77}_{-0.43}^{+0.96}$
${\eta }_{\oplus ,\mathrm{GK}}^{{\rm{O}}}$	${0.48}_{-0.27}^{+0.59}$ – ${0.59}_{-0.33}^{+0.76}$	${0.50}_{-0.27}^{+0.56}$ – ${0.70}_{-0.39}^{+0.90}$	${0.41}_{-0.22}^{+0.46}$ – ${0.51}_{-0.27}^{+0.57}$	${0.49}_{-0.25}^{+0.49}$ – ${0.59}_{-0.32}^{+0.68}$
${\eta }_{\oplus ,\mathrm{FGK}}^{{\rm{O}}}$	${0.54}_{-0.32}^{+0.81}$ – ${0.69}_{-0.42}^{+1.17}$	${0.57}_{-0.33}^{+0.81}$ – ${0.91}_{-0.55}^{+1.70}$	${0.46}_{-0.27}^{+0.64}$ – ${0.62}_{-0.35}^{+0.93}$	${0.57}_{-0.30}^{+0.70}$ – ${0.77}_{-0.45}^{+1.28}$

Notes. The low and high bounds correspond to the constant- and zero-completeness extrapolation of Section 3.3.2. The superscripts C and O on η_⊕ refer to the conservative and optimistic HZs. η_⊕,GK is the HZ occurrence for GK stars (3900 K ≤ T_eff ≤ 6000 K), and η_⊕,FGK is the HZ occurrence for FGK stars (3900 K ≤ T_eff ≤ 7300 K).

Download table as:  ASCIITypeset image

While this study takes care to incorporate detection and vetting completeness and importantly both reliability and observational uncertainty, there are still unresolved issues. We summarize these issues here, each of which can motivate future improvements to our occurrence rate model and methodology.

Power-law assumption. Products of power laws in terms of radius and period are commonly adopted for planet population models in occurrence rate studies, but there is growing evidence that calls their suitability into question. For instance, improvements to stellar radius measurements have revealed that the radius distribution for small, close-in planets is bimodal, rather than a smooth or broken power law (Fulton et al. 2017), which has also been observed in K2 data (Hardegree-Ullman et al. 2020). Power laws are not capable of describing such non-monotonic populations. Looking at the bottom panels of Figure 10 (without uncertainty), we can see some data points in the radius and instellation flux distributions do not lie along our inferred power laws. However, using input uncertainties (top panels of Figure 10) washes out this structure, making it more difficult to discern the failure or success of a power-law model as a descriptor of the data. There is also strong evidence that populations are not well described by products of power laws in terms of radius and period (Petigura et al. 2018; Lopez & Rice 2018) for orbital periods of <100 days. Therefore a product of power laws such as Equation (5) in terms of radius and instellation flux is unlikely to be a good description of planet populations at high instellation. At low instellation of the HZ, however, the observed planet candidate population does not indicate any obvious structure (see Figure 1). Given that our domain of analysis is plagued by few detections, low completeness, and low reliability, more observations are likely needed to determine more appropriate population models. Therefore, because most of our planet population have radii larger than 1.5 R_⊕, those larger planets are likely driving the population model and may be causing bias in the model in the smaller-planet regime due to our simple product power laws in Equation (5).

Planetary multiplicity. Zink et al. (2019) pointed out that when short-period planets are detected in the Kepler pipeline, data near their transits are removed for subsequent searches, which can suppress the detection of longer-period planets around the same star. They found that for planets with periods greater than 200 days detection completeness can be suppressed by over 15% on average. Our stellar population has removed stars for which more than 30% of the data has been removed due to transit detection via the dutycycle_post stellar property from the DR25 stellar property table (for details see Bryson et al. 2020a). We do not attempt to quantify the extent to which this stellar cut mitigates the impact identified in Zink et al. (2019), nor do we account for this effect in our analysis.

Stellar multiplicity contamination. Several authors (Ciardi et al. 2015; Furlan et al. 2017; Furlan & Howell 2017, 2020) have shown that undetected stellar multiplicity can impact occurrence rate studies in at least two ways. Stellar multiplicity can reveal planet candidates to be false positives, reducing the planet population, and valid planet candidates in the presence of unknown stellar multiplicity will have incorrect planet radii due to flux dilution. They estimate that these effects can have an overall impact at the 20% level. Stellar multiplicity can also bias the parent stellar sample because unaccounted-for flux dilution will bias the completeness estimates. Our analysis does not take possible stellar multiplicity into account. However stellar multiplicity has been shown to be associated with poor-quality metrics, specifically the BIN flag of Berger et al. (2018) and the Gaia RUWE metric (Lindegren 2018). For example, A. Kraus et al. (2021, in preparation) found that few Kepler target stars with RUWE > 1.2 are single stars. As described in Section 3.1, we remove stars from our parent stellar population that have been identified as likely binaries in Berger et al. (2018) or have RUWE > 1.2, which is expected to remove many stars with undetected stellar multiplicity (for details see Bryson et al. 2020a). We do not attempt to quantify the impact of undetected stellar multiplicity for our stellar population after this cut.

Planet radius and HZ limits. There are several stellar, planetary, and climate model–dependent factors that could reduce the occurrence rates calculated in this work. It is quite possible that uncertainties in the stellar radii may alter the planet radii, moving some rocky planets into the mini-Neptune regime of >1.5 R_⊕. Or, it is possible that the upper limit of 1.5 R_⊕ is an overestimate of the rocky planet limit, and the rocky-to-gaseous transition may lie lower than 1.5 R_⊕. Although, as pointed out in Section 2 following Otegi et al. (2020), the rocky regime can extend to as large as 2.5 R_⊕, many of these large–radius regime planets are highly irradiated ones, so they may not be relevant to HZ rocky planets.

The HZ limits themselves may be uncertain, as they are model and atmospheric composition dependent. Several studies in recent years have calculated HZ limits with various assumptions (see Kopparapu et al. 2019 for a review). In particular, the inner edge of the HZ could extend further in, closer to the star, due to slow rotation of the planet (Yang et al. 2014; Kopparapu et al. 2016; Way et al. 2016), and the outer edge of the HZ may shrink due to "limit cycling," a process where the planet near the outer edge of the HZ around FG stars may undergo cycles of globally glaciated and unglaciated periods with no long-term stable climate state (Kadoya & Tajika 2014, 2015; Menou 2015; Haqq-Misra et al. 2016). Consequently, the number of planets truly in the HZ remains uncertain.

Our computation of η_⊕ has large uncertainties, with the 68% credible interval spanning factors of 2 (see Tables 3–5). The 99% credible intervals in Table 6 span two orders of magnitude. In Section 5.3 we discuss how comparing occurrence rates with and without input uncertainties in Tables 3 and 4 indicates that these large uncertainties are present before considering the impact of uncertainties in the input data. We also observe in Table 5 that the uncertainties are considerably smaller for planets larger than those contributing to our η_⊕. We conclude that while input uncertainties make a contribution, the dominant cause of our large uncertainties is Poisson uncertainty due to the very small number of HZ planets smaller than 1.5 R_⊕ in very low completeness regions of the DR25 planet catalog (see Figure 1). Our uncertainties may be close to a noise floor induced by the small number of small HZ planets resulting from the low completeness.

These large Poisson-driven uncertainties are unlikely to be reduced by resolving the issues discussed in Section 5.10. Only by increasing the small-planet catalog's completeness, resulting in a larger small-planet HZ population, can these uncertainties be reduced.

There are two ways in which a well-characterized catalog with more small planets can be produced:

• 1.  
  Develop improved planet-vetting metrics that produce a catalog that is both more complete and more reliable than the DR25 catalog. There are several opportunities for such improved metrics, discussed in Bryson et al. (2020a), such as more fully exploiting pixel-level data and existing instrumental flags that can give more accurate reliability characterization than that given using DR25 products. This approach requires new vetting metrics. Bryson et al. (2020b) have shown that varying the DR25 Robovetter thresholds does not significantly change occurrence rates or their uncertainties once completeness and reliability are taken into account. In Appendix D we show that such changes in Robovetter metrics also do not significantly change the occurrence rates we find in this paper.
• 2.  
  Obtain more data with a quality similar to Kepler's, likely through more space-based observations. In Section 1 we describe how the unrealized Kepler extended mission, doubling the amount of data relative to DR25, was expected to significantly increase the yield of small planets in the HZ. An additional 4 years of data from observing the same stars as Kepler with similar photometric precision would be sufficient. Eight years of observation on a different stellar population would also suffice. As of this writing, plans for space-based missions such as TESS and PLATO do not include such long stares on a single field. For example, PLATO is currently planning no more than 3 years of continuous observation of a single field.⁵⁵