Exoplanet population inference and the abundance of Earth analogs from noisy, incomplete catalogs

Introduction

NASA’s Kepler mission has enabled the discovery of thousands of exoplanet candidates (Batalha et al. 2013; Burke et al. 2014). While many of these candidates have not been confirmed as bona fide planets, there is evidence that the false positive rate is low (Morton & Johnson 2011; Fressin et al. 2013), enabling conclusions about the population of planets based on the catalog of candidates. Many of these planets orbit Sun-like stars (Petigura et al. 2013b), where the definition of Sun-like is given in terms of the star’s temperature and surface gravity. Given these catalogs, it is interesting to ask what we can say about the population of exoplanets as a function of their physical parameters (period, radius, etc.). Observational constraints on the population can inform theories of planet formation and place probabilistic bounds on the abundance of Earth analogsFor our purposes, an “Earth analog” is an Earth-sized exoplanet orbiting a Sun-like star with a year-long period..

Petigura et al. (2013b) recently published an exoplanet population analysis based on an independent study of the Kepler light curves for 42,557 Sun-like stars. This study was especially novel because the authors developed their own planet search pipeline (TERRA; Petigura et al. 2013a) and determined the detection efficiency of their analysis empirically by injecting synthetic signals into real light curves measured by Kepler. The occurrence rate function determined by Petigura et al. (2013b) agrees qualitatively with previous studies of small planets orbiting Sun-like stars (Dong & Zhu 2013). In particular, both papers describe a “flattening” rate function (in logarithmic radius) for planets around Earth’s radius. Even though no Earth analogs were discovered in their search, Petigura et al. (2013b) used the small candidates that they did find to place an extrapolated constraint on the frequency of Earth-like exoplanets, assuming a flat occurrence rate density in logarithmic period.

A very important component of any study of exoplanet populations is the treatment of detection efficiency. Speaking qualitatively, in a transit survey, small planets with long periods are much harder to detect than large planets orbiting close to their star. This effect is degenerate with any inferences about the rate density and it can be hard to constrain quantitatively. In practice, there are three methods for taking this effect into account: (a) making conservative cuts on the candidates and assuming that the resulting catalog is complete (Catanzarite & Shao 2011; Traub 2012; Tremaine & Dong 2012), (b) asserting an analytic form for the detection efficiency as a function of approximate signal-to-noise (Youdin 2011; Howard et al. 2012; Dressing & Charbonneau 2013; Dong & Zhu 2013; Fressin et al. 2013; Morton & Swift 2013), and (c) determining the detection efficiency empirically by injecting synthetic signals into the raw data and testing recovery (Christiansen et al. 2013; Petigura et al. 2013a, b).

There are two qualitatively different methods that are commonly used to infer the occurrence rate density from a catalog and a detection efficiency specification. The first is an intuitive method that we will refer to as “inverse-detection-efficiency” and the second is based on the likelihood function of the catalog given a parametric rate density. The inverse-detection-efficiency method involves making a histogram of the objects in the catalog where each point is weighted by its inverse detection probability. This method is very popular in the literature (Howard et al. 2012; Dong & Zhu 2013; Dressing & Charbonneau 2013; Swift et al. 2013; Petigura et al. 2013b) even though it is not motivated probabilistically. The alternative likelihood method models the catalog as a Poisson realization of the observable rate density of exoplanets taking the survey detection efficiencies and transit probabilities into account. This technique has been used to constrain parametric models—a broken power law, for example—for the occurrence rate density (Tabachnik & Tremaine 2002; Youdin 2011; Dong & Zhu 2013). In this Article, we start from the likelihood method but model the rate density non-parametrically as a piecewise-constant step function. Using this formulation of the problem, we derive a generalization that takes observational uncertainties into account. In Appendix A, we show that the inverse-detection-efficiency method can be derived as a special case of the likelihood method in the limit of a smoothly varying completeness function.

In every previous study of exoplanet occurrence rates, the authors have assumed that the measurement uncertainties are negligible. This assumption is not justified because these uncertainties—especially on measurements (like exoplanet radius) that depend on the stellar parameters—can be large compared to the scales of interest. In this Article, we develop a flexible framework for probabilistic inference of exoplanet occurrence rate density that can be applied to incomplete catalogs with non-negligible observational uncertainties. Our method takes the form of a hierarchical probabilistic (Bayesian) inference. We generalize the method introduced by Hogg et al. (2010b) to account for survey detection efficiencies. We run tests on simulated datasets—comparing results with the standard techniques that neglect observational uncertainties—and apply our method to a real catalog of small planets transiting Sun-like stars (Petigura et al. 2013b).

For the purposes of this Article, we make some strong assumptions, although we argue that they are weaker than the implicit assumptions in previous studies. None of these assumptions is necessary for the validity of our general method but they do simplify the specific procedures we employ. We assume that

the candidates in the catalog are independent draws from an inhomogeneous Poisson process set by the censored occurrence rate density,

every candidate is a real exoplanet (there are no false positives),

the observational uncertainties on the physical parameters are non-negligible but known (the catalog provides probabilistic constraints on the parameters),

the detection efficiency of the pipeline is known, and

the TrueIn this Article, we use “True” to describe an observable (for example, the exoplanet occurrence rate density) that would be trivially measured in the limit of very high signal-to-noise data. We use “true” to describe a simulation quantity with a value exactly known to us. occurrence rate density is smoothWe give our definition of “smooth” in more detail below but our model is very flexible so this is not a strong restriction..

The first assumption—conditional independence of the candidates—is reasonable since the dataset that we consider explicitly includes only single transiting systems (Petigura et al. 2013b). The second assumption—neglecting false positives—is also strong and only weakly justified by estimates of low false positive rates in the Kepler data (Fressin et al. 2013; Morton & Johnson 2011). For this Article, we will neglect this issue and only comment on the effects but the prior distributions published by Fressin et al. (2013) could be directly applied in a generalization of our method.

We must emphasize one very important consequence of our assumptions. We assume that the catalog of exoplanet candidates is only missing planets with probabilities given by the empirical detection efficiency. In detail this must be false; the detection efficiency we use doesn’t take into account the fact that the catalog doesn’t include multiple transiting systems. A large fraction of the transiting planets discovered by the Kepler transit search pipeline are members of multiple transiting systems (see Lissauer et al. 2011, for example). Since Petigura et al. (2013b) only detected at most one planet per system, their catalog is actually a list of planet candidates without a more detectable companion. The global effects of this selection are not trivial and an in-depth discussion is beyond the scope of this Article but all of the results should be interpreted with this caveat in mind.

Conditioned on our assumptions and the choices made in the planet detection, vetting and characterization pipeline (Petigura et al. 2013a, b), we constrain the rate density of small exoplanets orbiting Sun-like stars. As part of this analysis we also place probabilistic constraints on the rate densityIn this Article, we use the word “rate” to indicate the dimensionless expectation value of a Poisson process and the words “rate density” to indicate a quantity that must be integrated over a finite bin in period and radius to deliver a rate. of Earth analogs $\Gamma_{\oplus}$, which we define as the expected number of planets per star per natural logarithmic bin in period and radius, evaluated at the period and radius of Earth

Since no Earth analogs have been detected, this constraint requires an extrapolation in both period and radius. Petigura et al. (2013b) performed this extrapolation by assuming that the period distribution of planets in a small bin in radius is flat, obtaining ${\Gamma_{\oplus}}\approx 0.12$. We relax this assumption and extrapolate only by assuming that the occurrence rate density is a smooth function of period and radius; we find lower values for $\Gamma_{\oplus}$. We enforce the smoothness constraint by applying a flexible Gaussian process regularization to the bin heights.

In the next Section, we summarize the likelihood method for exoplanet population inference and in Section 3, we describe how to include the effects of observational uncertainties. The technical term for this procedure is hierarchical inference and while a general discussion of this field is beyond the scope of this Article, in Section 3, we present the basic probabilistic question and derive a computationally tractable inference procedure. In Section 4, we summarize the technique and derive the key equation for our method: Equation (16). We test our method on synthetic catalogs in Sections 6 and 7. In Section 8, we use the catalog of planet candidates and the empirically determined detection efficiency from Petigura et al. (2013b) to measure the occurrence rate density of small planets with long orbital periods.

Sections 2 and 3 provide a general pedagogical introduction to the methods used in this Article. Readers looking to implement a population inference are directed to Appendix A if measurement uncertainties are negligible or Section 4 (especially Equation 16) for problems with non-negligible uncertainties. Readers interested in our results—the inferred population of exoplanets and Earth-analogs—can safely skip to Section 8 and continue to the discussion in Section 9.

The likelihood method

The first ingredient for any probabilistic inference is a likelihood function; a description of the probability of observing a specific dataset given a set of model parameters. In this particular project, the dataset is a catalog of exoplanet measurements and the model parameters are the values that set the shape and normalization of the occurrence rate density. Throughout this Article, we use the notation $\Gamma_{\boldsymbol{{\theta}}}({\boldsymbol{w}})$ for the occurrence rate density $\Gamma$—parameterized by the parameters $\boldsymbol{{\theta}}$—as a function of the physical parameters $\boldsymbol{w}$ (orbital period, planetary radius, etc.). In this framework, the occurrence rate density can be “parametric”—for example, a power law—or a “non-parametric” function—such as a histogram where the bin heights are the parameters $\boldsymbol{{\theta}}$.

We’ll model the catalog as a draw from the inhomogeneous Poisson process set by the observable rate density $\hat{\Gamma}_{\boldsymbol{{\theta}}}$. This leads to the previously known result (see Tabachnik & Tremaine 2002; Youdin 2011 for some of the examples from the exoplanet literature)

In this equation, the integral in the normalization term is the expected number of observable exoplanets in the sample.

The main thing to note here is that $\hat{\Gamma}_{\boldsymbol{{\theta}}}$ is the rate density of exoplanets that you would expect to observe taking into account the geometric transit probability and any other detection efficiencies. In practice, we can model the observable rate density as

Finally, we model the rate density as a piecewise constant step function

where the parameters ${\theta}_{j}$ are the log step heights and the bins ${\Delta}_{j}$ are fixed a priori. In Appendix A, we use this parameterization and derive the analytic maximum likelihood solution for the step heights. This result is similar to and just as simple as the inverse-detection-efficiency method and it is guaranteed to provide a lower variance estimate of the rate density than the standard procedure.

One major benefit of expressing the problem of occurrence rate inference probabilistically is that it can now be formally extended to include the effects of observational uncertainties.

A brief introduction to hierarchical inference

The general question that we are trying to answer in this Article is: what constraints can we put on the occurrence rate density of exoplanets given all the light curves measured by Kepler? In the case of negligible measurement uncertainties, this is equivalent to optimizing Equation (2) but when this approximation is no longer valid, we must instead compute the marginalized likelihood

where $\{{\boldsymbol{x}}_{k}\}$ is the set of all light curves, one light curve ${\boldsymbol{x}}_{k}$ per target $k$, $\boldsymbol{{\theta}}$ is the vector of parameters describing the population occurrence rate density $\Gamma_{\boldsymbol{{\theta}}}({\boldsymbol{w}})$ and ${\boldsymbol{w}}_{k}$ is the vector of physical parameters describing the planetary system (orbital periods, radius ratios, stellar radius, etc.) around target $k$. In this equation, our only assumption is that the datasets depend on the rate density of exoplanets only through the catalog $\{{\boldsymbol{w}}_{k}\}$. In our case, this assumption qualitatively means that the signals found in the light curves depend only on the actual properties of the planet and star, and not on the distributions from which they are drawn. It is worth emphasizing that—as we will discuss further below—the catalog only provides probabilistic constraints on $\{{\boldsymbol{w}}_{k}\}$; not perfect delta-function measurements.

In other words, we treat the catalog as being a dimensionality reduction of the raw data with all the relevant information retained. In the context of Kepler, the catalog reduces the set of downloaded time series (approximately 70,000 data points for the typical Kepler target) to probabilistic constraints on a handful of physical parameters—$\boldsymbol{w}$ from above—like the orbital period and planetary radius. If we take this set of parameters $\{{\boldsymbol{w}}_{k}\}$ as sufficient statistics of the data then we can, in theory, compute Equation (10)—up to an unimportant constant—without ever looking at the raw data again! This is important because the high-dimensional integral in Equation (10) won’t generally have an analytic solution and each evaluation of the per-object likelihood $p({\boldsymbol{x}}_{k}\,|\,{\boldsymbol{w}}_{k})$ is expensive, making numerical methods intractable.

Instead, we will reuse the hard work that went into building the catalog. We must first notice that each entry in a catalog is a representation of the posterior probability

of the parameters ${\boldsymbol{w}}_{k}$ conditioned on the observations of that object ${\boldsymbol{x}}_{k}$. The notation $\boldsymbol{\alpha}$ is a reminder that the catalog was produced under a specific choice of a—probably “uninformative”—interim prior $p({\boldsymbol{w}}_{k}\,|\,{\boldsymbol{\alpha}})$. This prior was chosen by the author of the catalog and it is different from the likelihood $p({\boldsymbol{w}}_{k}\,|\,{\boldsymbol{{\theta}}})$ from Equation (2).

Now, we can use these posterior measurements to simplify Equation (10) to a form that can, in many common cases, be evaluated efficiently. To find this result, multiply the integrand in Equation (10) by

The data only enter this equation through the posterior constraints provided by the catalog $\{{\boldsymbol{w}}_{k}\}$! For our purposes, this is the definition of hierarchical inference.

The constraints in Equation (11) can always be—and often are—propagated as a list of $N$ samples $\{{\boldsymbol{w}}_{k}\}^{(n)}$ from the posterior

We can use these samples and the Monte Carlo integral approximation to estimate the marginalized likelihood from Equation (13)—up to an irrelevant constant—as

where the constant $Z_{\boldsymbol{\alpha}}=p(\{{\boldsymbol{x}}_{k}\}\,|\,{\boldsymbol{\alpha}})$ is not a function of the parameters $\boldsymbol{{\theta}}$. This is very efficient to compute as long as an evaluation of $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{{\theta}}})$ is not expensive. That being said, Equation (15) could be a high variance estimator of Equation (13), depending on the number of independent samples $N$ and the initial choice of $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{\alpha}})$. Additionally, the support of $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{{\theta}}})$ in $\{{\boldsymbol{w}}_{k}\}$ space is restricted to be narrower than that of $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{\alpha}})$. Besides this caveat, in the limit of infinite samples, the approximation in Equation (15) becomes exact. Equation (15) is the importance sampling approximation to the integral in Equation (13) where the trial density is the posterior probability for the catalog measurements.

A very simple example is the familiar procedure of making a histogram. If you model the function $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{{\theta}}})$ as a piecewise constant rate density—where the step heights are the parameters—and if the uncertainties on the catalog are negligible compared to the bin widths then the maximum marginalized likelihood solution for $\boldsymbol{{\theta}}$ is a histogram of the catalog entries. The case of non-negligible uncertainties is described by Hogg et al. (2010b) using a method similar to the one discussed here.

Model generalities

Now, we can substitute Equation (2) into Equation (13) and apply the importance sampling approximation (Equation 15) to derive the following expression for the marginalized likelihood

where the values $\{{{\boldsymbol{w}}_{k}}^{(n)}\}$ are samples drawn from the posterior probability

as described in the previous section. Equation (16) is the money equation for our method. It lets us efficiently compute the marginalized likelihood of the entire set of light curves for a particular occurrence rate density.

In this equation, we’re making the further assumption that the catalog treated the objects independently. This is a somewhat subtle point if we were to consider targets with more than one transiting planet—a point that we will return to below—but for the considerations of the dataset considered here, it is a justified simplification.

For the remainder of this Article, we model the rate density as a two-dimensional histogram with fixed logarithmic bins in period and radius. When we include observational uncertainties—using Equation (16)—the maximum likelihood result is no longer analytic. Therefore, if we want to compute the “best-fit” rate density, we can use a standard non-linear optimization algorithm.

In the regions of parameter space that we tend to care about, the completeness is low and there are only a few observations with large uncertainties. In this case, we’re especially interested in probabilistic constraints on the occurrence rate density; not just the best-fit model. To do this, we must apply a prior $p({\boldsymbol{{\theta}}})$ on the rate density parameters and generate samples from the posterior probability

There is a lot of flexibility in the choice of functional form of $p({\boldsymbol{{\theta}}})$. In the well-sampled parts of parameter space there are a lot of detected planets and the choice of prior makes little difference, but in the regions that we care about, the detection efficiency is low and applying a prior that captures our beliefs about the rate density is necessary. This will be especially important when we extrapolate the rate density function to the location of Earth—in Section 7—where no exoplanets have been found. Therefore, instead of using an uninformative prior, we want to use a prior that encourages the occurrence rate density to be “smooth” but it should be flexible enough to capture structure that is supported by the data. To achieve this, we model the logarithmic step heights as being drawn from a Gaussian process (Rasmussen & Williams 2006; Gibson et al. 2012; Ambikasaran et al. 2014). This model encodes our prior belief that, on the grid scale that we consider, the rate density should be smooth but it is otherwise very flexible about the form of the function.

Mathematically, the Gaussian process density is

where $\Sigma^{-1}$ is the diagonal matrix

The Gaussian process model for the step heights given in Equation (19) is very flexible but the results will depend on the values of the hyperparameters $\mu$ and $\boldsymbol{{\lambda}}$. Therefore, instead of fixing these parameters to specific values, we add another level to our hierarchical probabilistic model and marginalize over this choice. In other words, we apply priors—uniform in the logarithm—on $\mu$ and $\boldsymbol{{\lambda}}$, and sample from the joint posterior

Strictly speaking, in this model, $p({\boldsymbol{{\theta}}}\,|\,{\mu},\,{\boldsymbol{{\lambda}}})$ can’t really be called a “prior” anymore and the constraints on the step heights are no longer independent.

There is an efficient algorithm called elliptical slice sampling (ESS; Murray et al. 2010; Murray & Prescott Adams 2010) for sampling the step heights $\boldsymbol{{\theta}}$ from the density in Equation (24). In practice, for problems with this specific structure, ESS outperforms more traditional MCMC methods commonly employed in astrophysics (e.g., Foreman-Mackey et al. 2012). Our implementation is adapted from Jo Bovy’s BSD licensed ESS codehttps://github.com/jobovy/bovy_mcmc/blob/master/bovy_mcmc/elliptical_slice.py. To simultaneously marginalize over the hyperparameter choice, we use the Metropolis–Hastings update from Algorithm 1 in Murray & Prescott Adams (2010). We tune the Metropolis–Hastings proposal by hand until we get an acceptance fraction of $\sim 0.2-0.4$ for the hyperparameters.

For all the results below, we run a Markov chain with $10^{6}$ steps for the heights and update the hyperparameters every 10 steps. We only keep the final $2\times 10^{5}$ steps and discard the earlier samples as burn-in. By estimating the empirical integrated autocorrelation time of the chain (Goodman & Weare 2010), we find that the resulting chain has $\gtrsim 4000$ independent posterior samples. These samples provide an approximation to the marginalized probability distribution for $\boldsymbol{{\theta}}$.

Data and completeness function

Using an independent exoplanet search and characterization pipeline, Petigura et al. (2013b) published a catalog of 603 planet candidates orbiting stars in their “Sun-like” sample of Kepler targets. For each candidate, Petigura et al. (2013b) used Markov chain Monte Carlo to sample the posterior probability density for the radius ratio, transit duration, and impact parameter assuming uninformative uniform priors. They then incorporated the uncertainties in the stellar radius and published constraints on the physical radii of their candidates. Given this data reduction and since we don’t have access to the individual posterior constraints on radius ratio and stellar radius, we can’t directly compute the importance weights $p(\{{\boldsymbol{w}}_{k}\}\,|\,{\boldsymbol{\alpha}})$ needed for Equation (15). For the rest of this Article, we’ll make the simplifying assumption that these weights are constant in log-period and log-radius but the results don’t seem to be sensitive to this specific choice.

Petigura et al. (2013b) did not publish or share posterior samples of their measurements of the physical parameter (Equation 14). They did publish a list of periods, radii and radius uncertainties based on their analysis. Assuming that there is no measurement uncertainty on the period measurement and that the radius posterior is Gaussian in linear radius (with a standard deviation given by the published uncertainty), we draw 512 samples for ${\boldsymbol{w}}_{k}$ and use these as an approximation to the posterior probability function.

A huge benefit of this dataset is that Erik Petigura and collaborators published a rigorous analysis of the empirical end-to-end completeness of their transit search pipeline. Instead of choosing a functional form for the detection efficiency of the pipeline as a function of the parameters of interest, Petigura et al. (2013b) injected synthetic signals of known period and radius into the raw aperture photometry and determined the empirical recovery after the full analysis.

We use all the injected samples from Petigura et al. (2013b) to compute the mean (marginalized) detection efficiency in bins of $\ln P$ and $\ln R$. In each bin, this efficiency is simply the fraction of recovered injections. For the purposes of this Article, we neglect the counting uncertainties introduced by the finite number of samples used to estimate the completeness. The largest injected signal had a radius of $16\,R_{\oplus}$ but, because of the measurement uncertainties on the radii, we need to model the distribution at larger radii. To do this, we approximate the survey completeness for $R>16\,R_{\oplus}$ as 1.

Given our domain knowledge of how detection efficiency depends on the physical parameters, the intuitive choice would be to measure the survey completeness in radius ratio or signal-to-noise instead of period and radius. It is also likely that a change of coordinates would yield a higher precision result. That being said, it is still correct to measure the completeness in period and radius, and there are a few practical reasons for our choice. The main argument is that since the radius uncertainties are dominated by uncertainties in the stellar parameters, it is not possible to use the published catalog (Petigura et al. 2013b) to compute constraints on radius ratios. In the future, this problem would be solved by publishing a representation of the full posterior density function for each object in the catalog. In this case, the most useful data product would be posterior samples for each target’s radius ratio and stellar radius.

The detection efficiency also depends on the geometric transit probability $R_{\star}/a$. Since we are modeling the distribution in the period–radius plane, we need to compute the transit probability marginalized over stellar radius and mass. This marginalized distribution scales only with the period of the orbit as $\propto P^{-2/3}$. In theory, this marginalization should be over the True distribution of these parameters in the selected stellar catalog but we’ll approximate it by the empirical distribution; a reasonable simplification given the size of the dataset. At a period of 10 daysThis period is chosen arbitrarily because the power law only needs to be normalized at one point., the median transit probability in the selected sample of stars is $5.061\%$ so we model the transit probabilityWe are using the letter $Q$ to indicate probabilities since we are already using $P$ to mean period. as a function of period as

Implicit in the expression for the transit probability in Equation (25) is the assumption that all of the planets are on circular orbits. Recently, Kipping (2014) demonstrated that when eccentric orbits are included, our given value is an underestimate by about 10%. This effect will propagate directly to our inferred rate densities. Even though the degeneracy is not exact—due to our choice of priors on the rate density parameters—it is not a bad approximation to assume that it is and scale the results down by your preferred factor. The right thing to do would be to marginalize over this effect directly during inference but that exercise is beyond the scope of the current Article. To complicate matters, the detection probability of a transit is also a non-trivial function of the duration. To account for this effect, so non-circular orbits should also be injected when measuring the survey completeness.

Validation using synthetic catalogs

In order to get a feeling for the constraints provided by our method and to explore any biases introduced by ignoring the observational uncertainties, we start by “observing” two synthetic catalogs from qualitatively different known occurrence rate density functions. For each of these simulations, we take the completeness function computed by Petigura et al. (2013b) as given. In general, Equation (2) can be sampled using a procedure called thinning (Lewis & Shedler 1979) but for our purposes, we’ll simply consider a piecewise constant rate density evaluated on a fine grid in log-period and log-radius. For this discrete function, the generative procedure is simple;

distribute $K_{i}$ catalog entries in the cell randomly.

We then choose fractional observational uncertainties on the radii from the Petigura et al. (2013b) catalog and apply them to the true catalog as Gaussian noise.

We generate synthetic catalogs from two qualitatively different rate density functions. Both distributions are generated by a separable model

but fit using the full general model. The first catalog—Catalog A—is generated assuming a smooth occurrence surface where both distributions are broken power laws. The second—Catalog B—is designed to be exactly the distribution inferred by Petigura et al. (2013b) in the range that they considered and then smoothly extrapolated outside that range. The catalogs generated from these two models are shown in LABEL:smooth-results and LABEL:simulation-results, respectively and the data are available onlinehttp://dx.doi.org/10.5281/zenodo.11507.

For each catalog, we directly apply both the inverse-detection-efficiency procedure as implemented by Petigura et al. 2013bOur implementation reproduces their results when applied to the published catalog. and our probabilistic method, marginalizing over the hyperparameters of the Gaussian process regularization. Figure 1 and LABEL:simulation-results show the results of this analysis in both cases. In particular, the side panels compare the marginalized occurrence rate density in period and radius to the true functions that were used to generate the catalogs. Figure 1 shows that even if the True rate density is a smooth function, the density inferred by the inverse-detection-efficiency method can appear to have sharp features. In this first example—where the true distribution is well described by our Gaussian process model—the probabilistic inference of the occurrence rate density is both more precise and accurate.

In the second example, the true rate density includes a sharp feature chosen to reproduce the result published by Petigura et al. (2013b). In this case, LABEL:simulation-results shows that the probabilistic constraints on the rate density are less precise but more accurate than results using the inverse-detection-efficiency method. This effect is most apparent in the parts of parameter space where the detection efficiency is low—long period and small radius.

When applied to either simulated catalog, the inverse-detection-efficiency method gives a high-variance estimate of the true occurrence rate density. One effect of this variance is that the inferred distribution will appear to have more small-scale structure than the true underlying distribution.

Extrapolation to Earth

As well as inferring the occurrence distribution of exoplanets, this dataset can also be used to constrain the rate density of Earth analogs. Explicitly, we constrain the occurrence rate density of exoplanets orbiting “Sun-like” starsIn this Article, we adopt the Petigura et al. (2013b) sample of G-stars as our definition of “Sun-like”., evaluated at the location of Earth:

That is, $\Gamma_{\oplus}$ is the rate density of exoplanets around a Sun-like star (expected number of planets per star per natural logarithm of period per natural logarithm of radius), evaluated at the period and radius of Earth.

In Equation (27), we use the symbol $\Gamma$ instead of the more commonly used $\eta$ since we define “Earth analog” in terms of measurable quantities with no mention of habitability or composition. This might seem unsatisfying but the composition of an exoplanet is notoriously difficult to measure even with large uncertainty and any definition of habitability is still extremely subjective. With this in mind, we stick to the observable definition for this Article.

Since no Earth analogs have been found, any constraints on this density must be extrapolated from the existing observations. This is generally done by assuming a functional form for the occurrence rate density, constraining it using the observed candidates and extrapolating. All published extrapolations are based on rigid models of the occurrence rate density (for example, a power law) fit to the catalog and evaluated at the location of Earth (Catanzarite & Shao 2011; Traub 2012). Petigura et al. (2013b) used their catalog of planet candidates to constrain the rate of Earth analogs in a specific period–radius bin assuming an extremely rigid model: flat in logarithmic period. These results are all sensitive to the choice of extrapolation function and the specific definition of “Earth analog”.

We weaken the assumptions necessary for extrapolation by only assuming that the distribution is smooth using the Gaussian process regularization described in Section 4. Under this model, the occurrence rate density at periods and radii where no objects have been detected will be constrained—with large uncertainty—by the heights of nearby bins. Therefore, even though there are no candidates that qualify as Earth analogs, we simply fit our model of the occurrence rate density in a large enough region of parameter space (including Earth) and compute the posterior constraints on $\Gamma_{\oplus}$. This works because the Gaussian process regularization actually captures our prior beliefs about the shape of the rate density function. This model—and any other extrapolation—will, of course, break down if there is an unmeasured sharp feature in the occurrence rate density near the location of Earth but our method is the most conservative extrapolation technique published to date.

Figures 3 and 4 compare our results and the results of the Petigura et al. (2013b) extrapolation procedure when applied to the synthetic catalogs. Since these catalogs were simulated from a known population model, we know the true value of $\Gamma_{\oplus}$ and it is indicated in the figures with a vertical gray line. In both cases, our method returns a less precise but more accurate result for the rate density and the error bars given by the functional extrapolation are overly optimistic. One major effect that leads to this bias is that the period distribution is not flat. Restricting the result to only include uniform models is equivalent to applying an extremely informative prior that doesn’t have enough freedom to capture the complexity of the problem. As a result, the posterior constraints on $\Gamma_{\oplus}$ are dominated by this prior choice and the resulting uncertainties are much smaller than they should be.

Results from real data

Having developed this probabilistic framework for exoplanet population inferences and demonstrating that it produces reasonable results when applied to simulated datasets, we now turn to real data. As described in Section 5, we will use the catalog of small exoplanet candidates orbiting Sun-like stars published by Petigura et al. (2013b). This is a great test case because those authors empirically measured the detection efficiency of their pipeline as a function of the parameters of interest.

We directly applied our method to the Petigura et al. (2013b) sample and generated MCMC samples from the posterior probability for the occurrence rate density step heights, marginalizing over the hyperparameters of the Gaussian process model. The resulting MCMC chain is available onlinehttp://dx.doi.org/10.5281/zenodo.11507.

Figure 5 shows posterior samples from the inferred occurrence rate density as a function of period and radius conditioned on the catalog. The marginalized distributions are qualitatively consistent with the occurrence rate density measured using the inverse-detection-efficiency method with larger uncertainties.

The period distribution integrated over various radius ranges is shown in LABEL:period. In agreement with Dong & Zhu (2013), we find that the period distribution of large planets ($R>8\,R_{\oplus}$) is inconsistent with the distribution of smaller planets. The rate density of large planets appears to monotonically increase as a function of log period while the distribution for small planets seems to turn over at a relatively short period (around 50 days) and decrease for longer periods.

The equivalent results for the radius distribution are shown in Figures 7 and 8. Figure 7 shows the log-radius occurrence rate density integrated over various logarithmic bins in period. The distributions in each period bin are qualitatively consistent; the rate density is dominated by small planets (around two Earth radii) with potential “features” near $R\sim 3R_{\oplus}$ and $R\sim 10R_{\oplus}$. These features appear in every period bin. They were also detected—using a completely different dataset and technique—by Dong & Zhu (2013) and a similar result is visible in the occurrence rate determined by Fressin et al. (2013, their Figure 7) at low signal-to-noise. Figure 8 shows the same result but presented as a function of linear radius. In these coordinates, the rate density in a single bin is no longer uniform; instead, scales as inverse radius.

Our constraint on the rate density of Earth analogs (as defined in Section 7) is in tension—even though our result has large fractional uncertainty—with the result from Petigura et al. (2013b). This is shown in LABEL:real-rate where we compare the marginalized posterior probability function for $\Gamma_{\oplus}$ to the published value and uncertainty. Quantitatively, we find that the rate density of Earth analogs is

Although they are mainly nuisance parameters, we also obtain posterior constraints on the hyperparameters $\mu$ and $\boldsymbol{{\lambda}}$. In particular, the constraints on the length scales in $\ln P$ and $\ln R$ are ${\lambda}_{P}=3.65\pm 1.03$ and ${\lambda}_{R}=0.65\pm 0.12$ respectively. Both of these scales are larger than a bin in their respective dimension. For completeness we also find the following constraints on the other hyperparameters

The MCMC chains used to compute these values is available onlinehttp://dx.doi.org/10.5281/zenodo.11507.

Comparison with previous work

Our inferred rate density of Earth analogs (Equation 29) is not consistent with previously published results. In particular, our result is completely inconsistent with the earlier result based on exactly the same dataset (Petigura et al. 2013b). This inconsistency is due to the different assumptions made and the detailed cause merits some investigation. The two key differences between our analysis and previous work are (a) the form of the extrapolation function, and (b) the presence of measurement uncertainties on the planet radii.

To make their estimate of $\Gamma_{\oplus}$, Petigura et al. (2013b) asserted a flat distribution in logarithmic period for small planets. Our results suggest that the data do not support this assumption (see LABEL:period). We find that the data require a decreasing period distribution in the relevant range. A similar result was also found by Dong & Zhu (2013) and it is apparent in Figure 2 of Petigura et al. (2013b).

With the large error bars, this result is consistent with both results (see LABEL:comparison where this value is labeled “linear extrapolation”) but it does not fully account for the discrepancy.

To examine the effects of measurement uncertainties, we repeat our analysis with the error bars on the radii artificially set to zero, keeping everything else the same. This analysis (labeled “uncertainties ignored” in LABEL:comparison) gives the result

This result is relatively more precise and higher than our final result and consistent with the value obtained with linear extrapolation. This confirms the hypothesis that the discrepancy between our result and the previously published values is the combined result of both of our key generalizations.

Discussion

We have developed a hierarchical probabilistic framework for inferring the population of exoplanets based on noisy incomplete catalogs. This method incorporates systematic treatment of observational uncertainties and detection efficiency. One major benefit of this framework is that it provides the best possible probabilistic measurements of the population under the assumptions listed in Section 1 and repeated below. After demonstrating the validity of our method on two qualitatively different synthetic exoplanet catalogs, we run our inference on a published catalog of small exoplanet candidates orbiting Sun-like stars (Petigura et al. 2013b) to determine the occurrence rate density these planets as a function of period and radius. We extrapolate this measurement to the location of Earth and constrain the rate density of Earth analogs with large error bars. In order to perform this extrapolation, we don’t assume a specific functional form for the rate density. Instead, we only assume that it is a smooth function of logarithmic period and radius.

The occurrence rate density function that we infer is qualitatively consistent with previously published results using different inference techniques (Dong & Zhu 2013; Fressin et al. 2013; Petigura et al. 2013b). In particular, we find (see LABEL:radius) previously recorded features in the radius distribution around $R\sim 3\,R_{\oplus}$ and $R\sim 10\,R_{\oplus}$, although not at high signal-to-noise. We find that the period distributions for planets in different radius bins are different, in qualitative agreement with previous results (Dong & Zhu 2013). Figure 6 shows that larger planets tend to be on longer periods than smaller planets.

Our extrapolation of the rate density to the location of Earth is more general and conservative than any previously published method. We find a rate density of Earth analogs that is inconsistent with the result published by Petigura et al. (2013b). This discrepancy can be attributed to both the rigidity of the assumptions about the period distribution and the effects of non-negligible measurement uncertainties. Our extrapolation is also less confident than previous measurements. Again, this difference is due to the fact that we allow a much more flexible extrapolation function. This is another illustration that, against the standard data analysis folklore, the correct use of flexible models is conservative.

In contrast to previous work, we don’t define “Earth analog” in terms of habitability or composition. Instead, we advocate for a definition in terms of more directly observable quantities (in this case, period and radius). Furthermore, we define $\Gamma_{\oplus}$ as a rate density (per star per logarithmic period per logarithmic radius) so that its value doesn’t depend on choices about the “Earth-like” bin.

In our analysis we make a few simplifying assumptions. Every assumption has an effect on the results and could be relaxed as an extension of this project. For completeness, we list and discuss the effects of our assumptions below.

Conditional independence We assume that every object in the catalog is a conditionally independent draw from the observable occurrence rate density. This is a bad assumption when applying this method to a different catalog where multiple transiting systems are included. In practice, the best first step towards relaxing this assumption is probably to follow Tremaine & Dong (2012) and assume that the mutual inclination distribution is the only source of conditional dependence between planets. For this Article, the assumption of conditional independence is justified because the dataset explicitly includes only systems with a single transiting exoplanet.

False positives In our inferences, we assume that all of the candidates in the catalog are True exoplanets. The rate of false positives in the Kepler catalog has been shown to be low but not negligible (Morton & Johnson 2011; Fressin et al. 2013). Since some of the objects in the catalog are probably false positives, our inferences about the occurrence rate density are biased high but without explicitly including a model of false positives, it’s hard to say in detail what effect this would have on the distributions. In an extension of this work, we could incorporate the effects of false positives by switching to a mixture model (see Hogg et al. 2010a, for example) where each object is modeled as a mixture of True exoplanet and false positive. In this mixture model, the false positives would be represented using prior distributions similar to those used by Morton (2012) or Fressin et al. (2013).

Known observational uncertainties To apply the importance sampling approximation to the published catalog, we assume that the measurement uncertainties are known and, in this case, Gaussian. The assumption of normally distributed uncertainties could be relaxed given a sampling representation of the posterior probability function for the physical parameters (period, radius, etc.). There is recent evidence that the stellar radii of Kepler targets might, on average, be underestimated (Bastien et al. 2014), introducing another source of noise. It is possible to relax the noise model and include effects like this but inference would be substantially more computationally expensive.

Given empirical detection efficiency Petigura et al. (2013b) determined the end-to-end detection efficiency of their planet detection pipeline as a function of True period and radius by injecting synthetic signals into real light curves and testing recovery. We used these simulations as an exact representation of the detection efficiency of the catalog but there are several missing components. The biggest effect is probably the fact that this formulation doesn’t include the selection of only the most detectable signal in each light curve. This bias will be largest in the parts of parameter space where the baseline detection efficiency is lowest: at long periods and small radius. As a result, our inferences (and the results from Petigura et al. 2013b) about the occurrence rate of small planets on long periods is probably underestimated relative to Truth. In detail there is another limitation due to the fact that the stellar parameters are only known noisily and the transit light curve only constrains the radius ratio. This means that the marginalized detection efficiency should be measured as a function of radius ratio and the interpretation in terms of True radius is only approximately correct. Given the size of the dataset and the number of injection simulations, this effect should be small.

Smooth rate function Throughout our analysis, we make the prior assumption that the occurrence rate density is a smooth function of logarithmic period and radius. This model is useful because it allows us to make probabilistically justified inferences about the exoplanet population in regions of parameter space with low detection efficiency. The assumption that the rate density should be smooth is intuitive but there is no theoretical indication that it must be true at all scales. That being said, the Gaussian process regularization that we use to enforce smoothness is flexible enough to capture substantial departures from smooth if they were supported by the data.

Our assumptions are severe but we believe that this is the most conservative population inference method currently on the market.

where the uncertainties are only on the expectation value and don’t include the Poisson sampling variance. This is an exciting result because it means that, if we can improve the sensitivity of exoplanet search pipelines to small planets orbiting on long periods, then we should find some Earth analogs in the existing data. Furthermore, because of the treatment of multiple transiting systems in the catalog, the True expected number of transiting Earth-like exoplanets orbiting Sun-like stars is almost certainly larger than the values in Equation (34)!

Some of the caveats on the results in this paper are due to assumptions made for computational simplicity but a much more robust study would be possible given a complete representation of the posterior probability function for the physical parameters in the catalog. The use of MCMC to fit models to observations is becoming standard practice in astronomy and the results in many catalogs (including Petigura et al. 2013b) are given as statistics computed on posterior samplings. For the sake of hierarchical inferences like the method presented here, it would be very useful if the authors of upcoming catalogs also published samples from these distributions along with the value of their prior function evaluated at each sample. In this spirit, we have released the results of this paper as posterior samplingshttp://dx.doi.org/10.5281/zenodo.11507 for the occurrence rate density function.

All of the code used in this project is available from http://github.com/dfm/exopop under the MIT open-source software license. This code (plus some dependencies) can be run to re-generate all of the figures and results in this Article; this version of the paper was generated with git commit d56324d (2014-08-28).

APPENDIX

Appendix A Inverse-detection-efficiency

One huge benefit of the inverse-detection-efficiency procedure is its simplicity. Therefore, it’s worth noting that there is a probabilistically justified procedure that will always provide less biased results while being only marginally more complicated.

To motivate this derivation, let’s start by considering the following pathological example: a single bin where the completeness sharply drops from one to zero halfway across the bin. If we observe $K$ objects in this bin, we would have observed about $2K$ objects in a complete sample. If we apply the inverse-detection-efficiency procedure to this dataset, each sample will get unit weight because they are all found in the part of the bin where the completeness is one. Therefore, we would underestimate the true rate in the bin by half. It’s clear in this specific case that giving the points a weight of two would give a better solution and we’ll derive the general result below.

If we model the occurrence rate density as a histogram with $J$ fixed bin volumes ${\Delta}_{j}$ (Equation 9) then Equation (2) becomes

where the indicator function $\mathbf{1}[\cdot]$ is one if $\cdot$ is true and zero otherwise. Taking the gradient of this function with respect to $\boldsymbol{{\theta}}$ and setting it equal to zero, we find the maximum likelihood result

where $K_{j}$ is the number of objects that fall within the bin $j$. We estimate the uncertainty $\delta{\theta}_{j}$ on this value by examining the curvature of the log-likelihood function near the maximum and find

In our pathological example from above, the integral of the completeness function over the bin is $1/2$, giving each sample the expected weight of $2$. In more realistic cases, where the completeness function varies smoothly, the inverse-detection-efficiency result will begin to agree with Equation (A2) but the severity of this bias will be very problem dependent. Therefore, if you have a dataset with negligible observational uncertainties, we recommend that you always apply Equation (A2) instead of the standard inverse-detection-efficiency procedure. As the uncertainties become more significant, there is no longer an analytic result and the method derived in this Article is necessary.

References