A Distinct Population of Small Planets: Sub-Earths

By injecting simulated transit signals into real light curves and then recovering them using the same search pipeline, one can determine the so-called "pipeline efficiency" as a function of planet size, orbital period, and stellar properties. The Kepler team has published multiple studies of pipeline completeness (e.g., Christiansen et al. 2016; Thompson et al. 2018). The most recent ones are those by Burke & Catanzarite (2017a) and Christiansen et al. (2020) for the DR25 pipeline. The Christiansen et al. (2020) work is based on the so-called pixel-level transit injection tests. Limited by the computational capacities, each Kepler star only receives one planet injection, and the resulting completeness is an average over the entire star sample. In contrast, the much simpler flux-level tests (Burke & Catanzarite 2017a) can inject many different signals into a given star, and one can survey a range of planet parameters (radius and period) on a per-star basis, to obtain the detection efficiency, DE(p, r, θ_*), as a function of the orbital period, planet radius, and host-star properties. Such a model captures the target-to-target variations and is recommended for studies involving a small subset of Kepler targets, as is the case here. We will adopt this model and use the KeplerPorts python code (Burke & Catanzarite 2017b) to generate the pipeline efficiency.

In detail, KeplerPorts calculate the detection efficiency, DE(P, R_p , θ_*), based on the so-called multiple-event statistics, $\mathrm{MES}(P,{R}_{p},\bar{\theta })$ (Jenkins 2002; Christiansen et al. 2012), with $\bar{\theta }$ being stellar parameters (R_*, $\mathrm{log}g$ and T_eff). This statistic is closely related to the transit signal-to-noise ratio. The detection efficiency is also empirically found to correlate strongly with the so-called CDPP slope, the temporal behavior of the stellar noise. This dependence is taken into account by the Burke & Catanzarite (2017b) model, based on injection experiments. Furthermore, the method of "MES Smearing" is applied to transits with nonzero impact parameters.

The resultant detection efficiencies are then averaged over the N_* = 13,297 GK dwarfs in our sample. For the radius and period ranges of interests, these efficiencies are plotted in the left-hand panel of Figure 1. One observes that, inward of 16 days, transiting planets larger than 1.3R_⊕ can be detected fairly completely (DE ≥ 90%), and the efficiency drops below ∼50% for planets smaller than 0.9R_⊕. The right-hand panel presents the overall efficiency, η(p, r), one that also includes the geometric transit probability. This is calculated using Equation (1):

$$\begin{eqnarray}&&\eta (P,{R}_{p})=\displaystyle \frac{1}{{N}_{* }}\displaystyle \sum _{i=1}^{{N}_{* }}{\mathrm{DE}}_{i}(P,{R}_{p})\times {P}_{\mathrm{geo},i}(P,{R}_{p}),\end{eqnarray} \tag{ 1 }$$

where P_geo,i = 2R_*,i /a, and we adopt an estimate for the star mass from Berger et al. (2020b).

Zoom In Zoom Out Reset image size

Binary stars and some astrophysical noise can masquerade as transit signals. So, every TCE needs to be vetted before qualifying as a planet candidate. This can introduce a reduction in completeness as some genuine planets may be vetted out as false positives. This vetting efficiency is presented in Figure 2 for our ranges of interest, produced using the robotic vetting code Robovetter (Coughlin 2017; Thompson et al. 2018). For the short-period planets we are interested in, the vetting efficiency has an average of 93%, with no obvious dependence on planet radius or period. So, in this study, we simply assume this to be 100% for the entire space.

Zoom In Zoom Out Reset image size

Morton et al. (2016) assigned an FPP for each planet candidate using an automated validation procedure. This uses a full range of observational constraints to calculate the probabilities of various hypothetical scenarios (genuine planets, brown dwarfs, eclipsing binaries, etc). We have discarded all false positives, but retained candidates in both the "confirmed" and "candidate" categories, with the proviso that FPP < 0.1. We follow convention and define a reliability index =(1 − FPP). These are shown in Figure 3 by the colors of the dots. Because reliability is high in our sample, we simply ignore this distinction and assume all candidates are genuine.

Zoom In Zoom Out Reset image size

One straightforward, nonparametric way of obtaining the occurrence rate is called "IDEM" (Howard et al. 2012; Petigura et al. 2013). Assuming the observed population is drawn from an underlying population weighted by the detection efficiency, one can retrieve the differential occurrence rate by dividing the number of planets detected in the bin ${\rm{\Delta }}\mathrm{ln}P{\rm{\Delta }}\mathrm{ln}{R}_{P}$, n(P, R_P ), by the detection efficiency,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}f}{d\mathrm{ln}{Pd}\mathrm{ln}{R}_{P}}=\displaystyle \frac{1}{{\rm{\Delta }}\mathrm{ln}P{\rm{\Delta }}\mathrm{ln}{R}_{P}}\times \displaystyle \frac{1}{{N}_{* }}\,\times \displaystyle \frac{n(P,{R}_{P})}{\eta (P,{R}_{P})},\end{eqnarray} \tag{ 2 }$$

where η(P, R) is the detection efficiency, suitably averaged over the relevant grid, and N_* is the stellar sample size. We then sum over the period bins to obtain a size distribution

$$\begin{eqnarray}&&\displaystyle \frac{{df}}{d\mathrm{ln}{R}_{P}}=\sum \displaystyle \frac{{d}^{2}f}{d\mathrm{ln}{Pd}\mathrm{ln}{R}_{P}}{\rm{\Delta }}\mathrm{ln}P.\end{eqnarray} \tag{ 3 }$$

The dominant source of error is Poisson error, as opposed to measurement errors. To obtain an estimate on the uncertainty, we note that because we are dealing with small candidate numbers in each bin (some bins even contain zero detection), the confidence interval cannot be estimated using the usual $1/\sqrt{n}$ approach. Instead, we adopt Table 6 of Feldman & Cousins (1998) to obtain the 95% confidence limits when presented with a detection n(P, R_P ) ≥ 0. When summing over different period bins, we assume that the confidence limits from each bin are independent and adds the errors quadratically. The uncertainty range will then increase with bin number. This, we believe, yields a maximal uncertainty range.

This approach departs from previous IDEM studies (e.g., Howard et al. 2012; Petigura et al. 2013; Kunimoto & Matthews 2020) in providing a more conservative error estimate. To obtain more realistic and smaller uncertainties, one needs to adopt a model-based approach (e.g., Bayesian; see below) that accounts for correlated information among bins, while paying a price for extra assumptions.

The result of this analysis is presented in tandem with the Bayesian model in Figure 5.

4.2.1. Model Ansatz

To conduct a Bayesian analysis, we need to make an array of assumptions for the underlying distribution.

First, we will depart from tradition and consider the distribution in planet mass and period. We use the planet core radius (R_c ) as a proxy for its mass (roughly, ${M}_{p}\propto {R}_{c}^{4}$; Fortney et al. 2007). Because mini-Neptunes and super-Earths are of the same origin and therefore likely share the same mass distribution, it seems sensible to combine their signals together by working on planet mass. Moreover, this enables us to separate the distribution in period–radius space (see below).

While R_c ≈ R_P for super-Earths, mini-Neptunes need some explanation. According to the theory of photoevaporation, the radii of the two populations are spaced by a ratio that is nearly 2. We adopt the number 1.7 taken from Owen & Wu (2017) and Wu (2019). This ratio is the result of fitting the observed distributions with evaporation theory. It is also the ratio between the observed super-Earth peak (∼1.4 R_⊕) and mini-Neptune peak (∼2.4 R_⊕) in our sample. So in our study, for any planets lying above 2R_⊕ (the so-called "Fulton" gap), we simply divide their observed radii by 1.7 to obtain their core radii. When correcting for detection completeness, we note that while super-Earths are smaller and are in principle more incomplete compared to mini-Neptunes, the difference is minor and can be ignored (Figure 3).

One naturally wonders if this procedure produces artificial features. We show that folding Neptunes into super-Earths, while enhancing the super-Earth peak, is not critical to our conclusions on the existence and characters of sub-Earths.

Next, we assume that the planet distribution in the P–R_c plane is separable,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}f}{d\mathrm{ln}P\,d\mathrm{ln}{R}_{c}}=\displaystyle \frac{{{df}}_{P}}{d\mathrm{ln}P}\times \displaystyle \frac{{{df}}_{R}}{d\mathrm{ln}{R}_{c}}.\end{eqnarray} \tag{ 4 }$$

Such a separability assumption, often adopted in the literature (e.g., Silburt et al. 2015; Mulders et al. 2018; Neil & Rogers 2020), is problematic when describing the distribution of super-Earths and/or mini-Neptunes individually, because their relative proportion depends on distances from the star. It is, however, reasonable when describing the joint distribution. Studies of the planet-mass distribution (e.g., Wu 2019) have shown that it is relatively distance dependent. This is another reason behind our adoption of R_c as the relevant parameter.

In the following, we will experiment with different forms of f_R , while using the same broken power law for f_P ,

$$\begin{eqnarray}\displaystyle \frac{{{df}}_{P}}{d\mathrm{ln}P}=\left\{\begin{array}{ll}{\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{\alpha }, & \mathrm{if}\ P\leqslant {P}_{b}\\ {\left(\displaystyle \frac{{P}_{b}}{1\mathrm{day}}\right)}^{\alpha -\beta }{\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{\beta }, & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where P_b , α, and β are free parameters to be determined from data. This form is motivated by previous studies of super-Earths/mini-Neptunes (e.g., Silburt et al. 2015; Mulders et al. 2018; Neil & Rogers 2020). We assume that any new population of planets also satisfies the same rules. We relax this assumption in Model 4.

For the planet-mass (or core-radius) distributions, we adopt the following three models, each with a different physical significance.

Model 1: one Gaussian. Here, we assume that nature produces only one single population of planets, for core radius from 0.6 ⊕ to 2 ⊕ (or from ∼0.2 to 16M_⊕, if these planets all share the terrestrial composition). We describe this using a log-normal distribution in core radius,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{df}}_{R}}{d\mathrm{ln}{R}_{c}} & = & A\times N(\mathrm{ln}{R}_{c},{\mu }_{0},{\sigma }_{0})\\ & = & A\times \displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{0}^{2}}}\exp \left[-\displaystyle \frac{{\left(\mathrm{ln}{R}_{c}-\mathrm{ln}{\mu }_{0}\right)}^{2}}{2{\sigma }_{0}^{2}}\right].\end{array}\end{eqnarray} \tag{ 6 }$$

We adopt broad linear priors for both μ₀ and σ₀: μ₀ ∈ [0.6, 2]R_⊕ and σ₀ ∈ [0, 1].

Model 2: two Gaussians. To allow for a separate population of planets, we assume the core-radius distribution is a superposition of two Gaussians,

$$\begin{eqnarray}&&\displaystyle \frac{{{df}}_{R}}{d\mathrm{ln}{R}_{C}}=B\times N(\mathrm{ln}{R}_{c},{\mu }_{1},{\sigma }_{1})+C\times N(\mathrm{ln}{R}_{c},{\mu }_{2},{\sigma }_{2}),\end{eqnarray} \tag{ 7 }$$

where the first term describes the super-Earth/Neptune population, and we adopt linear priors of μ₁ ∈ [1.3, 1.6]R_⊕ and σ₁ ∈ [0, 0.5]. These narrower priors are intuited from both the above IDEM analysis and previous studies (e.g., Hsu et al. 2019; Neil & Rogers 2020). For the second term, employed to describe a new lower-mass population, we assign broad linear priors of μ₂ ∈ [0.6, 1.2] R_⊕ and σ₂ ∈ [0, 1].

Model 3: power-law+Gaussian. While super-Earth/Neptunes appear to be well described by (narrow) Gaussians, this does not appear so for the new population of planets we are after. The above two-Gaussian results, as well as our IDEM analysis, suggest that they may have a broad radius range, compatible with a power-law distribution. We attempt the following power-law/Gaussian combination,

$$\begin{eqnarray}&&\displaystyle \frac{{{df}}_{R}}{d\mathrm{ln}{R}_{C}}=B\times N(\mathrm{ln}{R}_{c},{\mu }_{1},{\sigma }_{1})+D\times {\left(\displaystyle \frac{{R}_{C}}{{R}_{\oplus }}\right)}^{\gamma }.\end{eqnarray} \tag{ 8 }$$

We select a linear prior of γ ∈ [ − 3, 1]. This is not restrictive, as it allows for both ascending and descending trends toward small sizes. Moreover, neither this model nor the two-Gaussian model presumes the presence of a gap-like feature.

Model 4: two periods. The previous three models adopt the same period distribution for any planet population. Here we relax this assumption by allowing for a separate period distribution for each population (based on model 3),

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{d}^{2}f}{d\mathrm{ln}{Pd}\mathrm{ln}{R}_{C}} & = & B\times N(\mathrm{ln}{R}_{c},{\mu }_{1},{\sigma }_{1})\times {\left(\displaystyle \frac{{{df}}_{P}}{d\mathrm{ln}P}\right)}_{I}\\ & & +D\times {\left(\displaystyle \frac{{R}_{c}}{{R}_{\oplus }}\right)}^{\gamma }\times \left(\displaystyle \frac{{{df}}_{P}}{d\mathrm{ln}P}\right).\end{array}\end{eqnarray} \tag{ 9 }$$

Here both period distributions have the same form as Equation (5) but they can have different parameters. Moreover, because super-Earths/Neptunes extend to well beyond 16 days and their period distribution has been well determined (see, e.g., Neil & Rogers 2020), we can adopt narrower priors on their parameters based on previous studies (see Table 2).

Table 2. The Best-fit Value of Parameters for Each Occurrence Rate Density Model with 68% Credible Interval and the Corresponding Akaike Information Criteria

Parameter	1 Gaussian	2 Gaussian	Power-Gaussian	2 Period	Prior
A	${0.008}_{-0.002}^{+0.003}$	⋯	⋯	⋯	[0, 0.05]
μ0	${1.31}_{-0.07}^{+0.13}$	⋯	⋯	⋯	[0.6, 2]
σ0	${0.50}_{-0.08}^{+0.17}$	⋯	⋯	⋯	[0, 1]
B	⋯	${0.0029}_{-0.0010}^{+0.0013}$	${0.0027}_{-0.0008}^{+0.0010}$	${0.0031}_{-0.0012}^{+0.0014}$	[0, 0.01]
μ1	⋯	${1.44}_{-0.03}^{+0.02}$	${1.43}_{-0.02}^{+0.02}$	${1.43}_{-0.02}^{+0.02}$	[1.3, 1.6]
σ1	⋯	${0.14}_{-0.02}^{+0.03}$	${0.14}_{-0.02}^{+0.02}$	${0.13}_{-0.02}^{+0.02}$	[0, 0.5]
C	⋯	${0.005}_{-0.003}^{+0.003}$	⋯	⋯	[0, 0.01]
μ2	⋯	${0.74}_{-0.09}^{+0.12}$	⋯	⋯	[0.6, 1.2]
σ2	⋯	${0.6}_{-0.3}^{+0.3}$	⋯	⋯	[0, 1]
D	⋯	⋯	${0.0028}_{-0.0008}^{+0.0010}$	${0.0030}_{-0.0009}^{+0.0015}$	[0, 0.01]
γ	⋯	⋯	$-{1.2}_{-0.7}^{+0.5}$	$-{0.9}_{-0.6}^{+0.4}$	[−3, 1]
α	${2.4}_{-0.264}^{+0.293}$	${2.3}_{-0.2}^{+0.3}$	${2.3}_{-0.3}^{+0.3}$	${2.5}_{-0.4}^{+0.3}$	[1, 3]
β	${0.72}_{-0.18}^{+0.16}$	${0.73}_{-0.19}^{+0.18}$	${0.73}_{-0.19}^{+0.18}$	$-{0.0}_{-0.5}^{+0.4}$	[−1, 1]
Pb	${4.1}_{-0.5}^{+0.5}$	${4.2}_{-0.5}^{+0.5}$	${4.2}_{-0.6}^{+0.5}$	${4.4}_{-0.6}^{+0.6}$	[2, 10]
αI	⋯	⋯	⋯	${1.9}_{-0.2}^{+0.2}$	[1.5, 3]
βI	⋯	⋯	⋯	${0.3}_{-0.3}^{+0.2}$	[−0.5, 0.5]
Pb,I	⋯	⋯	⋯	${8.8}_{-1.3}^{+0.9}$	[5, 10]
AIC	2841	2802	2799	2801
relative likelihood	∼10−10	0.2	1	0.4

Note. Best models are the power-Gaussian as it has the lowest AIC. The last row is the likelihoods of other models relative to the power-Gaussian one.

Download table as:  ASCIITypeset image

4.2.2. Bayesian Operations

Here, we apply Bayesian inference to determine the parameters of the maximum posterior probability.

Given the observed properties of a planet sample ("PCs"), the posterior probability of a hypothesis θ (with its attendant parameters) can be calculated using the Bayes' theorem as

$$\begin{eqnarray}&&P({\boldsymbol{\theta }}| \{\mathrm{PCs}\})=P(\{\mathrm{PCs}\}| {\boldsymbol{\theta }})\pi ({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 10 }$$

where P({PCs}∣ θ ) is the likelihood and π( θ ) the prior. We use exclusively flat priors.

As the planets are discrete points in the period–radius space, we treat them as a Poisson point process, following earlier works (e.g., Youdin 2011; Burke et al. 2015; Bryson et al. 2020). Defining the occurrence rate within a bin ${\rm{\Delta }}\mathrm{ln}P\,{\rm{\Delta }}\mathrm{ln}{R}_{C}$ as

$$\begin{eqnarray}&&{\rm{\Delta }}f({\boldsymbol{\theta }})=\displaystyle \frac{{d}^{2}f}{d\mathrm{ln}{Pd}\mathrm{ln}{R}_{C}}\,{\rm{\Delta }}\mathrm{ln}P\,{\rm{\Delta }}\mathrm{ln}{R}_{C}\,,\end{eqnarray} \tag{ 11 }$$

where θ stands for an ansatz and its attendant parameters, we can write the total expected number of planet detection as

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\exp } & = & {N}_{* }\displaystyle \int \displaystyle \int \displaystyle \frac{{d}^{2}f}{d\mathrm{ln}{Pd}\mathrm{ln}{R}_{C}}\times \eta (P,{R}_{C})\,d\mathrm{ln}{Pd}\mathrm{ln}{R}_{C}\\ & = & {N}_{* }\displaystyle \sum _{\mathrm{bins}}{\rm{\Delta }}f({\boldsymbol{\theta }})\,\times \eta (P,{R}_{C}),\end{array}\end{eqnarray} \tag{ 12 }$$

where η(P, R_c ) is the averaged detection completeness within the bin (Section 3). We choose a fine grid, with 20 period bins (1–16 days) and 30 core-radius bins (0.6–2R_⊕).

The likelihood function is then

$$\begin{eqnarray}&&P(\{\mathrm{PCs}\}| {\boldsymbol{\theta }})=\exp (-{N}_{\exp })\displaystyle \prod _{i=1}^{{n}_{p}}{\rm{\Delta }}f({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 13 }$$

where the product is carried over each detected candidate.

In detail, we first use the package scipy.optimize to obtain values that maximize the likelihood function. We then sample the posteriors using the Markov Chain Monte Carlo (MCMC) algorithm provided by the emcee package (Foreman-Mackey et al. 2013). We use 50 walkers, 1000 "burn-in" steps, and 5000 MCMC steps to find convergence, where the walkers of "burn-in" are initiated in a narrow Gaussian whose mean is the maximum likelihood solution. The parameters corresponding to the median of posterior distribution are shown in Table 2, with Figure 4 illustrating the corresponding size distributions for our best models.

Zoom In Zoom Out Reset image size

To compare the relative merits of various models, we employ the Akaike information criterion (AIC), a measure of information loss. This takes account of both the quality of the fit and the simplicity of the model to avoid both under-fitting and over-fitting. Let the latter be quantified by the number of free parameters, n _θ , we write

$$\begin{eqnarray}&&\mathrm{AIC}=2{n}_{{\boldsymbol{\theta }}}-2\mathrm{ln}\left[P(\{\mathrm{PCs}\}| {\boldsymbol{\theta }})\right].\end{eqnarray} \tag{ 14 }$$

The best model has the lowest AIC score, and the relative merit (likelihood) between models is quantified as e^(−ΔAIC/2). So, any model that has an AIC score higher by ∼10 relative to the best model is clearly rejected.

Here, we describe the results of our Bayesian inferences and compare them against those from IDEM. Their AICs and relative likelihoods are presented in Table 2.

First, among the models we tested, the Bayesian analysis strongly rejects Model 1 (single population) and requires the presence of at least one other population of planets, in addition to super-Earths. Model 1 is extremely unlikely, with a relative likelihood of $\exp (-40/2)\sim {10}^{-10}$ compared to the two-population models.

Second, the two-period model (Model 4), while allowing for more flexibility in the period distributions, performs slightly worse than the models that insist on a single period distribution. In other words, there is no statistical evidence suggesting that the two populations (super-Earths and sub-Earths) have different period distributions within the limited range that we probe (1–16 days). In fact, the obtained period distributions all feature a power-law rise to a few days, and a roughly logarithmically flat distribution afterwards, analogous to earlier studies using more extended samples (Silburt et al. 2015; Mulders et al. 2018; Neil & Rogers 2020).

Third, Models 2 and 3 are statistically indistinguishable (Model 3 is slightly preferred). Their size and period distributions also look similar (Figure 5). This is because the Gaussian distribution of sub-Earths, in the two-Gaussian model, is very broad and is compatible with a power law. From now on, we choose to focus on Model 3, the so-called power-Gaussian model.

Zoom In Zoom Out Reset image size

In Model 3, super-Earths are centered around 1.4R_⊕ with a narrow width of σ₁ ∼ 0.13. In planet mass, these correspond to a centroid of ∼4M_⊕ and an FWHM of a factor of 3.2. Given that we are focusing on host stars of lower masses (median mass 0.89M_⊙), this centroid falls below that obtained by Wu (2019), M_p ≈ 8M_⊕(M_*/M_⊙). The reason for the discrepancy is unclear for now. The FWHM does agrees with that in Wu (2019). The super-Earth population appears to peak narrowly in core mass, a mystery of the formation.

What is of most interest to us is the size distribution of sub-Earths. The power-Gaussian model yields a size index $\gamma =-{1.2}_{-0.7}^{+0.5}$, or the number of sub-Earths most likely rises toward small sizes. This is in drastic contrast to super-Earths, which fall off steeply toward small sizes. In addition, not all objects bigger than Earth belongs to the super-Earth group. Even at the super-Earth peak (∼1.4 R_⊕), sub-Earths make a non-negligible contribution (∼24%). Lastly, when we integrate all planets inward of 16 days and with core sizes from 0.6 to 2R_⊕, we find that super-Earths and sub-Earths are comparable in number (21% ± 4% versus 26% ± 4%).

There is a possibility that the negative γ value is driven not by sub-Earths below 1R_⊕, but by larger planets. Following the advice of the reviewer, we experiment with a more sophisticated model that includes a Gaussian (for super-Earths) and a broken power law. This introduces two new parameters, the break radius, and a new power-law index. We find a break radius of 0.77R_⊕, with the power index above this size of −1.8 ± 0.8, even steeper than our original result ($-{1.2}_{-0.7}^{+0.5}$). Due to the lack of objects below this radius, the index below is very uncertain (${0.7}_{-1.7}^{+1.0}$). This allows us to conclude that the rising power law we find is natural to the data and is not an artifact of our adopted model.

To compare the parameter-free IDEM results against the Bayesian results, we first convert the observed planet radii to core radii, as detailed above. We then interpolate the IDEM results between bins, using a third-degree spline smoothing to obtain the desired core-size distribution, ${df}/d\mathrm{ln}{R}_{C}$. The comparison is shown in Figure 5. They generally agree with each other. And as expected, our treatment of uncertainties in the IDEM approach yields larger error bars compared to those from forward models.

Our analysis is performed on a small set of "clean" planet samples (298 in total). Here, we vary the sample criteria and observe how the results for the power-Gaussian model change. The priors remain the same as before. For each experiment, we remove (or add) one constraint on the property of the stars or PCS, keeping other cuts unchanged.

We also experiment with alternative priors and completeness treatment. All experiments yield similar results compared to our default case (see Figure 6). We comment on each briefly below.

• 1.  
  No stellar radius cut: we experiment with employing all stars, regardless of radii. We exclude stars that satisfy the following crude cuts,

  $$\begin{eqnarray}\left\{\begin{array}{ll}\tfrac{{R}_{* }}{{M}_{* }^{0.9}}\gt 1.3, & \mathrm{if}\ {M}_{* }\lt {M}_{\odot }\\ \tfrac{{R}_{* }}{{M}_{* }^{1.25}}\gt 1.3, & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 15 }$$

  as they are most likely evolved. There are then 48,642 stars and 794 PCs. The stellar sample has expanded by a factor of 3.7, while the planet sample only a factor of 2.7, leading to a slightly lower occurrence (Figure 6). Moreover, the super-Earth peak is broadened, due to the known correlation between planet mass and stellar mass. The sub-Earth slope is γ ≈ −1.1.
• 2.  
  No K_p cut: while dimmer stars return fewer super-Earths and sub-Earths, the completeness correction largely accommodates this correctly. There is no perceptible change to the occurrence rates.
• 3.  
  No FPP cut: retaining all objects, regardless of their FPP values, adds 36 new candidates to our list, but produces no change to the results. The newly added candidates do not make the gap features more conspicuous, presumably because they are more likely to be potential false positives. However, we believe that unconfirmed PCs with a low FPP value are reliable enough and the sub-Earth distribution in our default case is likely genuine.
• 4.  
  Alternative priors: we have chosen broad and flat priors previously. Here, we follow the insights from Wu (2019) and use a narrower prior on super-Earths: μ₁ is now normally distributed with a mean at 1.43R_⊕ and a variance of 0.05R_⊕, and σ₁ is normally distributed with a mean of 0.1 and a variance of 0.1. The results remain the same.
• 5.  
  Alternative completeness: while we have adopted Burke & Catanzarite (2017a) to obtain the detection completeness, there exists an alternative approach (Christiansen et al. 2020). In this method, the completeness of a detection is expressed as a function of its MES (multiple-event-statistics) value, as well as some parameters a, b, c that are empirically obtained from pixel-level injection experiments,

  $$\begin{eqnarray}&&\eta (\mathrm{MES},a,b,c)=\displaystyle \frac{c}{{b}^{a}{\rm{\Gamma }}(a)}{\int }_{0}^{\mathrm{MES}}{t}^{a-1}{e}^{-\tfrac{t}{b}}{dt}.\end{eqnarray} \tag{ 16 }$$

  Because pixel-level injections are computationally expensive, typically one injection per star is possible. These empirical parameters are therefore the best fits for the entire Kepler stellar sample. This type of correction is used by, e.g., Neil & Rogers (2020) and Hsu et al. (2019). Here, we employ KeplerPORTS to map radius and period into MES for every planet host in our sample, and then adopt a = 29.41, b = 0.284, and c = 0.891 from Neil & Rogers (2020). Running the Bayesian analysis returns results that do not differ significantly from our original ones (see Figure 6).
• 6.  
  Reliability correction: as some planet candidates are less reliable (a higher false-positive probability) than 100%, we experiment with down-weighting them in the analysis and repeat the IDEM analysis. This makes negligible difference, presumably because all candidates are above 90% reliable, with most close to 100%.

Zoom In Zoom Out Reset image size

There are two recent works that employ updated completeness corrections and have included sub-Earths into their modeling efforts (Hsu et al. 2019; Neil & Rogers 2020). We compare our results with theirs here for the period range from 1 to 16 days. The Neil & Rogers (2020) work infers distributions in planet masses. We use simple scaling relations of $M/{M}_{\oplus }={({R}_{C}/{R}_{\oplus })}^{4}$ for R_C > 1.2R_⊕ and $M/{M}_{\oplus }={({R}_{C}/{R}_{\oplus })}^{3}$ for R_C ≤ 1.2R_⊕ to translate their mass distribution into the core-radius distribution. The other study (Hsu et al. 2019) infers the distribution in planet transit radius. We convert from these to the core radius as laid out in Section 4.2 and interpolate between bins where necessary. The comparisons are displayed in Figure 7, and we briefly describe any similarities or differences here.

Zoom In Zoom Out Reset image size

Similar to our study here, both works show a dropping off of super-Earths toward small sizes, and a new population below about 1R_⊕.

In the parameter-free, grid-based inversion model of Hsu et al. (2019), one observes a radius gap between super-Earths (including Neptunes) and sub-Earths. However, because their stellar sample covers a wide range of spectral types (mostly FGK dwarfs with no cuts in stellar radius), their super-Earth peak is broader and lower in amplitude, resembling closely our "no-radius-cut" case in Figure 6. By focusing on a narrow range of spectral types, we end up with a smaller planet sample but can detect more prominent features.

The forward-modeling approach of Neil & Rogers (2020) is very similar in spirit to our Bayesian models. By parameterizing the size distributions as a number of distinct Gaussians, they find that a new population of small planets (which they term "intrinsically rocky") is strongly demanded by the data. However, their overall size distribution does not show a clear inflection point, let alone a "gap" as we find here. We compare their results with our two-Gaussian model in detail. Their best-fit super-Earth peak lies at R_C = 1.53R_⊕ (slightly larger than our 1.44R_⊕), but with a much broader dispersion of σ = 0.63 (compared to our 0.14 ± 0.03). We attribute this latter difference to their wide range of stellar properties. They include almost all dwarfs, and because super-Earth properties are known to correlate with stellar properties, this tends to broaden the super-Earth peak. More significantly, their sub-Earths are described by a log-normal distribution with a peak at size R_C ≈ 0.96R_⊕ (compared to our 0.74) and a radius dispersion of σ ∼ 0.4 (same as our 0.4). The fact that their sub-Earths are more massive than ours, and that their super-Earths are more broadly distributed, give rise to a gap-free transition from super-Earths to sub-Earths (Figure 7), as well as a dropping occurrence toward very small sizes. These are qualitatively different from our results and those in Hsu et al. (2019). We believe that our choice of parameter space (narrower while preserving almost all sub-Earth candidates) is better suited for resolving any possible feature in the size distribution. A recent investigation, using their original model but with more relaxed assumptions, leads to better agreements with our results (A. Neil, 2020 private communication).

We have selected an optimum range to study sub-Earths. Our choices of 1–16 days and bright hosts minimize the uncertainties due to extreme completeness correction, without sacrificing many candidates. We have high confidence in our results because the corrected model actually does not look dramatically different from the raw data (Figure 8).

Zoom In Zoom Out Reset image size

Another support for our model comes from an intriguing class of objects, the so-called ultra-short-period planets (USP for short), planets that have orbital periods shortward of one day (e.g., Rappaport et al. 2013; Sanchis-Ojeda et al. 2014; Winn et al. 2018). At these short periods, Kepler can be sensitive to much smaller planets than otherwise possible. The origin for these boiling worlds is still under debate (e.g., Lee & Chiang 2017; Petrovich et al. 2019; Pu & Lai 2019; Millholland & Spalding 2020), though scenarios involving dynamical migration by companion planets have received some observational support (Sanchis-Ojeda et al. 2014; Adams et al. 2020). If so, their mass distribution should be similar to those farther away. Here, we exploit this advantage to shed light on our model. ²

Figure 8 compares the observed size distribution of USPs (<1 day) against that of our default sample (1–16 days) and against our preferred model. These 127 USPs are discovered by the Kepler main mission, are smaller than 2R_⊕, and include both "confirmed" and "candidate" dispositions. Their host magnitudes are mostly between 14 and 16 mag. We have also updated their radii using Gaia stellar parameters (Berger et al. 2020b).

Compared to our planet sample, the USPs show an excess of sub-Earths, thanks to their better completeness in these small sizes. More importantly, one can see a near-perfect agreement between the uncorrected USPs and our power-Gaussian model. If anything, the USPs suggest that the sub-Earth rise may be steeper than our current model, after completeness correction is properly applied.

Other than supporting our model, this result also suggests that USPs are indeed formed similarly as planets farther away.

In our sample, there are 71 PCs with radii smaller than 1.1R_⊕. These most likely belong to the sub-Earth population, according to our power-Gaussian model. Among them, 33 sub-Earths have at least one (though frequently more than one) transiting sibling planet. These multiple-planet systems are shown in Figure 9. We see that sub-Earths appear to happily coexist with larger planets (super-Earths, mini-Neptunes, even Jovians), with no substantial differences in both their period distribution and their stellar mass distribution.

Zoom In Zoom Out Reset image size

The fact that about half of all sub-Earths have transiting siblings may indicate that the majority of them prefer to reside in systems with other larger planets.