The Exoplanet Census: A General Method, Applied to Kepler
Abstract
We develop a general method to fit the planetary distribution function (PLDF) to exoplanet survey data. This maximum likelihood method accommodates more than one planet per star and any number of planet or target star properties. Application to Kepler data relies on estimates of the efficiency of discovering transits around Solar type stars by Howard et al. (2011). These estimates are shown to agree with theoretical predictions for an ideal transit survey. Using announced Kepler planet candidates, we fit the PLDF as a joint powerlaw in planet radius, down to , and orbital period, up to 50 days. The estimated number of planets per star in this sample is —, where the broad range covers systematic uncertainties in the detection efficiency. To analyze trends in the PLDF we consider four planet samples, divided between shorter and longer periods at 7 days and between large and small radii at 3 . At longer periods, the size distribution of the small planets, with index steepens to for the larger planet sample. For shorter periods, the opposite is seen: smaller planets follow a steep powerlaw, that is much shallower, at large radii. The observed deficit of intermediatesized planets at the shortest periods may arise from the evaporation and sublimation of Neptune and Saturnlike planets. If the trend and explanation hold, it would be spectacular observational confirmation of the core accretion and migration hypotheses, and allow refinement of these theories.
Subject headings:
Methods: statistical — Planetary Systems — Planets and satellites: detection — Planets and satellites: dynamical evolution and stability — Planets and satellites: formation — Stars: statistics1. Introduction
Individual exoplanet discoveries highlight the extraordinary diversity of worlds in the Solar neighborhood (Mayor & Queloz, 1995; Marcy & Butler, 1996; Charbonneau et al., 2009; Lissauer et al., 2011a). For an accurate census of the planet population, the statistical analysis of large samples of exoplanets is required (Cumming et al., 2008; Howard et al., 2010). Trends revealed by the data provide powerful constraints on dynamical theories of planet formation (Goldreich et al., 2004; Kenyon & Bromley, 2006; Youdin, 2010).
The correlation of giant planets with host star metallicity is perhaps the most interesting trend revealed by radial velocity surveys (Gonzalez, 1997; Fischer & Valenti, 2005; Johnson et al., 2010). This trend has been a powerful guide to identifying the mechanisms responsible for planetesimal formation (Youdin & Shu, 2002; Youdin & Goodman, 2005; Johansen et al., 2009), see Chiang & Youdin (2010) for a review.
The Kepler transit survey is currently revolutionizing the field of exolanets from space (Borucki et al., 2010). Borucki et al. (2011, hereafter BKep11) released 1,235 partially vetted planet “candidates.” Morton & Johnson (2011) estimate the false positive rate as around the brighter stars that are considered here. For brevity and statistical purposes, we mostly refer to the candidates as “planets,” though significant followup work remains with only 18 candidates currently confirmed.
Howard et al. (2011, hereafter HKep11) presented a detailed statistical analysis of the BKep11 planets around bright, Solartype stars. Accounting for detection efficiencies, HKep11 report planet occurrence rates of planets per star for radii and orbital periods . The purpose of this paper is to apply different statistical methods to the Kepler dataset, making use of the estimates of detection efficiencies — whose importance is on par with the detections themselves — provided by HKep11.
Our method is based on the Tabachnik & Tremaine (2002, hereafter TT02) technique to analyze the mass and period distributions of planets from radial velocity surveys. TT02 emphasized the importance of determining planet mass and period distributions simultaneously. A simultaneous fit is crucial because radial velocity detection thresholds depend on both the mass (times the sine of inclination) and period. Transit surveys have a similar, though weaker, coupling between planet radius and orbital period in the efficiency of discovering transits.
A primary advantage of the TT02 technique is that it naturally accommodates more that one planet per star. Planet hosting is not treated as a binary proposition. This distinction is not just a technical nicety, because the number of planets per star is order unity, at least. We show that Kepler data already imply about one planet per Solar type star within AU and with a radius bigger than . This number is certain to grow with time as longer periods and smaller sizes are probed.
High multiplicity rates in the Kepler sample (Latham et al., 2011) means that the fraction of stars that host a planetary system (at most one, of course) can be significantly smaller than the number of planets per star (NPPS). We treat the planet population (including known multiples) as a whole, and ignore multiplicity issues. Ragozzine & Holman (2010) discuss the value of multitransiting systems in constraining mutual inclinations. Using these constraints, Kepler data show that the number of planets per planetary system is at least 2 — 3 (Lissauer et al., 2011b).
This paper is organized as follows. Section 2 describes the detection efficiencies relevant to the Kepler transit survey, including an analytic fit to discovery efficiencies reported by HKep11 in §2.2. Section 3 describes our method for analyzing exoplanet survey data. Section 4 applies the method by fitting Kepler data to a joint powerlaw in radius and period. Differences between the shorter and longer period planets (in §4.2) and also between the smaller and larger planets (in §4.3) are analyzed. Section 5 gives nonparametric estimates of planet occurrence, and comments speculatively on the consequences of rising period distributions. We conclude with a summary and discussion in §6. Appendix A gives an analytic solution for powerlaw fits to data from an idealized transit survey.
2. Selection Effects For Transit Surveys
A robust statistical analysis must account for relevant selection effects. We quantify selection effects as detection efficiencies, , which give the ratio of detections to actual planets. Individual efficiencies multiply to give the net detection efficiency.
The three main selection effects for transit surveys are: (i) the transit probability for a planetary orbit to cross our line of sight to the star, ; (ii) the discovery efficiency of the survey in finding transiting planets, ; and (iii) the rate of false positive events that mimic a planet transit, . For false positives, the relevant , the only efficiency that can exceed unity. As mentioned already, the false positive rate is estimated to be low enough that we ignore it for simplicity.
For most analyses, an average detection efficiency is insufficient. The dependence of the efficiencies on relevant properties of planets and target stars is required. We consider how the efficiencies vary with planet radius, , and orbital period, . See §3.3 for a general discussion of which parameters should be included.
The left panel of Fig. 1 shows the Kepler discovery efficiency, calculated as described below (in section 2.2). The discovery efficiency is nearly unity over a signifiant range of planetary radii and periods, thanks to the high photometric precision of Kepler. There is a sharp drop in discovery efficiency for sufficiently small planets whose transit depths compete with noise in the photometric data. At shorter periods, the signal from smaller planets can rise above the noise because more transits are seen. The right panel of Fig. 1 shows how the geometic transit probability reduces the probability of finding long period planets.
To connect with the extensive work on radial velocity (RV) surveys, we note that RV upper limits can be expressed as a detection efficiency for a given star. Most simply, for planet parameters above the detection threshold and otherwise. This step function can be smoothed to account for confidence intervals.
2.1. Transit Probability
The probability that a planet, on a randomly oriented circular orbit of semimajor axis , transits it’s host star, with radius and mean density , is
(1)  
for .
Eccentricity tends to increase the detectability of planets at fixed (Burke, 2008).
For simplicity, we approximate the mean stellar density (or more specifically the mean ) as Solar. The densities of planet hosts can be estimated from precise light curve parameters, though corrections for the eccentricity apply (Tingley et al., 2011). Furthermore, planet hosts have a lower average density than the target star population, on general. This bias arises because lower density stars have higher transit probabilities. For our purposes, uncertainties in the average stellar density modestly affects the magnitude (but not the shape) of the inferred planet distribution.
One might (incorrectly) expect multiplanet systems to require special values of . When mutual orbital inclinations are small, finding one planet does increase the odds that others will be found. However since the planetary systems as a whole are randomly oriented, low mutual inclinations also make it easier to miss an entire planetary system. The mutual inclination of planets within systems has no effect on the overall detection rate (aside from sampling noise, as always).
2.2. The Kepler Discovery Efficiency
Precisely quantifying the discovery efficiency of a survey is quite difficult and is ongoing work for the Kepler mission. For our study we rely on the estimates of discovery efficiency in HKep11. HKep11 quantified the discovery efficiency for a subsample of bright solartype stars with relatively high . The sample consists of stars that satisfy the effective temperature, surface gravity and Kepler magnitude cuts:
{subeqnarray}
& 4100 K≤T_eff ≤6100 K&
&4.0 ≤logg /(cm s^2) ≤4.9&
&K_P ≤15 mag& .
Section 2.3 discusses the planet candidates in this sample.
The data for is presented in Figure 4 of HKep11. Here the discovery efficiency is reported for planets with and on a loguniform grid. HKep11 report , the number of stars around which a planet with that radius and period could be detected, in the bottom left of each cell. The related values of are plotted in Fig. 2.
For bins with no detected planets, HKep11 do not report a discovery efficiency. This omission is not because planet detections are required to estimate efficiencies. Rather HKep11 applied efficiencies to the detections only, which differs from our approach of applying the efficiencies across all parameter space. The missing efficiencies in Figure 4 of HKep11 can be readily obtained by interpolation. Most of the empty bins correspond to small and large where . Thus the lack of detections in these bins carries heavy statistical significance.
Instead of merely interpolating the reported Kepler efficiencies, we are motivated by their smooth variation to find an analytic fit. We find an excellent joint powerlaw fit, not to itself, but to the related
(2) 
the ratio of discoverable to nondiscoverable planets. Since ranges from zero to arbitrarily large values, it is simpler to fit with a powerlaw. The restriction to is automatically satisfied.
Fig. 2 shows that the broken powerlaw is an excellent fit to reported discovery efficiencies. The powerlaw fit is
(3) 
below a break at
(4) 
The efficiency at the break is quite high, (and remarkably constant, probably reflecting the high S/N threshold). Note that drops to 50% ( in Fig. 2) at , rising roughly from 1 to 2 over the periods considered. When discovery efficiencies are small, , the powerlaw in Equation (3) applies directly to .
While the scalings for are of little practical concern for our study, we report the fit above the break for completeness:
(5) 
This fit describes a small fraction of noisy stars, but does not affect our results because in this regime.
The relevant fit in Equation (3) agrees well with theoretical estimates. For S/N limited detections, the combined efficiency should scale as
(6) 
see e.g. Equation (7) of Gaudi (2007). The small limit of Equation (3) gives , rather good agreement.
The efficiencies summarized here will likely be improved by more detailed studies that include extraction of simulated transits from the actual Kepler data analysis pipeline. Accurate estimates of the Kepler discovery efficiency for all relevant parameters — especially stellar — are the most crucial ingredient for determining the frequency of Earthlike planets.
2.3. The Planetary Sample
Tables 1 and 2 of BKep11 list the properties of 1,235 planet candidates and their host stars. We must restrict our attention to the planet candidates around host stars with known discovery efficiencies, . Thus we only consider host stars within the Solar sample of HKep11, set by Equation (2.2). We further restrict attention to detections with signaltonoise (S/N) 10 and days, as these are also conditions for the validity of . Applying these filters reduces the number of planet candidates from 1,235 to 566.
Our “full” planet sample considers
{subeqnarray}
0.5 days¡ P ¡ 50 days
0.5 R_⊕¡ R ¡ 20 R_⊕ .
Adding these radius and minimum period cuts only removes four planets (from 566 to 562). Instead of using the BKep11 reported values for , we compute from the reported values for the host star radius, , and .
We now discuss relevant differences between our planet sample and that of HKep11. HKep11 conservatively restricted attention to , because the discovery efficiency is high for these planets (except at longer periods). We extend our analysis down to , which includes essentially all small planets in the Solar sample. This extension involves (a) trusting the HKep11 reported efficiencies down to the level to which they were reported, and (b) extrapolating the efficiencies down to . As shown above, both the smooth variation of the reported efficiencies and their agreement with theory give us confidence in making the extrapolation. We also report results for samples of larger planets that are unaffected by this extrapolation.
HKep11 defined S/N slightly differently than the publicly reported BKep11 S/N values on which we rely. HKep11 define S/N based on only one quarter of data (Q2, Kepler’s first full quarter). By contrast BKep11 report S/N values up to Q5, so that roughly four times more data has been collected. With perfect noise scaling, the BKep11 S/N values should be roughly twice as high as those used by HKep11. Table 2 of HKep11 presents noise scalings for 4 planets with radii and day periods. The ideal scaling roughly holds for three objects, but S/N saturated and only rose modestly for the fourth. It’s unclear what the general noise scaling is, and how it varies with size and period.
We proceed conservatively by performing all our analyses on both a S/N and a S/N planet sample (as defined by BKep 11). These cases cover the complete range of possibilities between saturated noise and perfect scaling. Fortunately the adopted S/N threshold has a modest effect on quantitative results (most importantly, the number of planets per star) but no effect on the interesting trends we discuss.
We briefly note a discrepancy between the planet counts reported by BKep11 and HKep11. HKep11 report a sample of 438 planet candidates (of the 1,235 released by BKep11) that (a) have hosts that satisfy the stellar parameter cuts of Equation (2.2), (b) have planet parameters days and and (c) satisfy the HKep11 definition of S/N 10. If we apply cuts (a) and (b) only to the BKep11 tables, but allow all S/N, there are 378 candidates. That number increases to 393 using the reported instead of calculating it as . Thus there are at least 60 (or 45) planets that HKep11 say should appear in the BKep11 tables, but they do not.
Applying any S/N threshold to the BKep11 tables can only make the discrepancy larger. We made no choices that would cause us to miss planets in the BKep11 tables. Cuts were applied inclusively (including any values that fall on a boundary) and all the planets in any system were counted. We conclude that a discrepancy exists between the planet and/or stellar parameters used by HKep11 and reported by BKep11. This discrepancy does not affect our results as long as the estimates of the discovery efficiencies in HKep11 are sufficiently accurate.
3. A General Method for Exoplanet Statistics
This section describes a general method for estimating the planetary distribution function (PLDF) from the results of a survey with quantified detection efficiencies. The purpose of developing this formalism is twofold.
First, this method can investigate a PLDF that depends on properties of both the planet detections and the target stars. Further, the PLDF can take an arbitrary functional form. While the exoplanet community is unlikely to agree on a single statistical methodology anytime soon, at least such an approach is possible.
Second, this formalism allows the inclusion of all detected planets — including those in multiplanet systems. This inclusion is possible because we treat the PLDF as the average number of planets per star (NPPS) not the fraction of stars with planets (FSWP). Following TT02, our technique treats planet occurrence as a Poisson process, i.e. a series of independent random events. This technique naturally gives the NPPS.
By contrast, many statistical analyses treat planet hosting as a binary proposition, i.e. a star either does or does not have a detected planet. The use of binomial statistics naturally gives the FSWP. Including multiple planets per system in this type of analysis is formally incorrect, though discrepancies are small if the multiplicity rate is low. For more discussion of these points see Cumming et al. (2008), who rigorously analyze the FSWP in radial velocity surveys by only including the most detectable planet around any star.
For transit surveys, the FSWP is more difficult to determine and is not addressed here. The main issue is applying the transit efficiency correction. Should the missed planets that do not transit be assigned evenly among all stars, or into high order multiplanet systems? As mentioned in the introduction, multitransiting systems are a powerful constraint, as are transit timing variations and comparison to RV surveys (Ragozzine & Holman, 2010).
Before developing the general method for fits to a parameterized PLDF in Section 3.2, we consider nonparametric estimates of the NPPS in Section 3.1. Parametric fits also give an estimate of the planet occurrence, whose quality depends on the appropriateness of the assumed functional form. The main reason to consider parameterized fits is to identify and assess trends in the data — and if feeling bold, to extrapolate. Section 3.3 discusses the consequences of fitting some parameters and ignoring others.
3.1. NonParametric Estimates of Planet Occurrence
For a survey of stars, we divide planet detections into bins indexed by , with planet detections in a given bin. If the average detection efficiency per bin is , then the best estimate of the NPPS in each bin is
(7) 
the number of detected planets per star divided by the efficiency.
The total planet occurrence, , depends on bin size when is not constant. For arbitrary small bins, there is only one planet per bin, and the efficiency is evaluated for the parameters of each planet and its host. This small bin limit is the only truly assumptionfree way to estimate the planet fraction, and suggests that binning of data is never required.
A maximum likelihood analysis using Poisson statistics can reproduce the estimates of given by Equation (7), and develop intuition for the general method. The number of expected detections in all bins is
(8) 
The likelihood function for a Poisson process is
(9) 
which represents the product of the probabilities of each individual detection (the term in square brackets) times the probability of finding no additional planets in any bin (the exponential).
Any constant multiplicative factors can be ignored in the likelihood function, because they do not affect the maximization of likelihood. (Thus factorials do not appear in Equation (9).) Here constant means independent of the PLDF, i.e. the values. Remarkably, this freedom allows us to ignore the efficiencies that characterize the detections, but not in . The likelihood thus simplifies to
(10) 
The efficiencies now appear only in .
The maximum likelihood values of are the roots of , which reproduces the expected result of Equation (7). It is easier (analytically and numerically) to maximize .
A fundamental feature of the method is that efficiencies are used to predict the number of expected planets and are not directly applied to the detections.
3.2. Fitting the Planetary Distribution Function
We now describe a general method to fit a parameterized PLDF to unbinned data. We consider a PLDF that depends on planet properties, (an component vector) and stellar properties (with components). Usually the elements of these vectors are logarithms of measured quantities such as the orbital period or stellar mass.
The probability that a star with properties has a planet with properties that lie in a volume
(11) 
Defined this way, the integrated represents the NPPS.
We express the differential distribution,
(12) 
as an amplitude, , times a shape function, . The shape parameters, , can describe the behavior of individual planetary or stellar properties as well as any correlations. In the simplest case, the values of are powerlaw exponents.
The total number of planets around stars (indexed by ) is
{subeqnarray}
N_tot &=& C ∑_j^N_∗ ∫g(,_j;) d
&=& N_∗C∫F_∗() g(,;) dd .
where integration covers a specified volume of the phase space of . In Equation (3.2) the sum over stars is replaced by an integral over the stellar distribution function . We proceed with the integral notation to minimize indices, but direct summation over all stars is always possible (and preferred when feasible).
Consider a survey with a net detection efficiency
(13) 
that is a product of efficiencies from individual selection effects (indexed by ).
The number of expected planet detections is
{subeqnarray}
N_exp&=& N_∗C ∫η(,) F_∗() g(,;) dd ,
&≡& N_∗C F
The shape integral, ,
weights the shape function over both the stellar distribution and the efficiencies.
Maximizing Likelihood
Consider a survey that detects planets around stars. We define the likelihood of the data analogously to Equation (10) as
(14) 
which represents the probabilities of each individual detection, labelled by , times the probability that no other planets were detected, given by the exponential factor.
Applying Equations (11) and (12), we again eliminate unnecessary constants, here the phase space volume, to define the likelihood function as
(15) 
where refers to the star that hosts planet detection . We let represent the shape function evaluated for the properties of detection and express the log likelihood as
(16) 
The best fit normalization, , is found by maximizing the likelihood as the root of :
(17) 
using Equation (13).
Following TT02, we use this constraint to eliminate from the likelihood:
(18)  
where the constant term in curly brackets can again be ignored.
The best fit values of maximize as the roots of , which are
(19) 
These (the number of parameters in ) equations generally couple to each other and require numerical solution, but see §A. Equation (19) weights the shape function by detection efficiencies on the left hand side and over individual detections on the right hand side.
The basic steps in obtaining a best fit solution are as follows. First, select a form of the PLDF to fit the data and express in the form of Equation (12). Second, solve Equation (19) for the shape parameters as just described. Third, calculate the normalization via Equation (17) with evaluated using the best fit shape parameters. The best fit solution is now complete and can be used (e.g.) to estimate the NPPS as using Equation (3.2). Fourth, estimate the errors on the fit via likelihood contours given by Equation (16), as described in §3.2.3 below.
Interpretation
This formalism has a remarkably straightforward interpretation. The normalization matches the numbers of detected and expected planets, , as seen by comparing Equations (13) and (17).
The shape parameters similarly match the expected and detected averages of planet and host star observables. Consider the case of powerlaw distributions with a shape function
(20) 
Here and are the logarithms of quantities being fit to a powerlaw, and distinguishes the powerlaw exponents for planetary and stellar parameters. With this powerlaw shape function Equation (19) gives
(21) 
where
{subeqnarray}
{ ⟨⟩_F, ⟨⟩_F } &≡& 1F ∫{, } ηF_∗g dd
{ ⟨⟩_obs , ⟨⟩_obs } &≡& 1 Npl ∑_i^N_pl { _i , _i } .
Nonpowerlaw shape functions can investigate more detailed properties of the detections than simply the average value of observables.
Errors
To quantify the uncertainty in the shape parameters, we compare contours of the likehood to a Gaussian distribution. The maximum likehood, , is given by Equation (18) with the best fit . The 1 uncertainties are defined by the contour. The general contours follow . The uncertainty in the normalization is just the range of values taken by within the error ellipse of the .
The 1D errors on an individual shape parameter (a component of ) are usually given assuming other parameters are held at their best fit values. When there is strong covariance, and a very elongated error ellipse (not the case in this study), 1D errors are not very meaningful.
The errors on describe the precision of a fit, not the quality. A large planet sample precisely defines the average powerlaw even if the actual PLDF deviates strongly from a powerlaw. The quality of fit can be determined by KolmogovSmirnoff (KS) tests.
This derivation ignores uncertainties in the detection efficiencies and the planetary and stellar parameters. Here, we describe how to include these errors but do not include them in our analysis. Uncertainties associated with an individual detection, , affect the value of the shape function, , that appear in the likelihood analysis. To include these uncertainties, should be a weighted integral over the uncertainties
(22) 
where is an appropriately normalized Gaussian distribution (for instance) centered on the best fit values of the planet detection and its host star. Ignoring uncertainties is equivalent to setting the error distributions, , to functions.
Uncertainties in the target star properties affect the number of expected planets via the shape integral . Errors can be included similarly to Equation (22) as
{subeqnarray}
F =1 N∗∑_j^N_∗ ∫ϕ_j(’  _j) η(,’)g(,’) dd’
≈∫ϕ(’  )F_∗() η(,’)g(,’) dd’ d .
The two forms above show how errors could be applied to each star , or (more approximately) how average errors could be assigned to the stellar distribution function. The integration over both and in Equation (3.2.3) is potentially confusing. The integration over is over the stellar distribution and replaces the sum over individual stars in Equation (3.2.3), while the integral over covers the range of uncertainties allowed by the error distributions .
Estimates of the total number of planets should also include target star uncertainties, by similarly modifying Equation (3.2) (i.e. setting in Equation (3.2.3)).
Systematic uncertainties in the detection efficiencies are more difficult to quantify. Assigning detection efficiencies are the most important (and dangerous) part of any statistical analysis. It is safest to perform multiple analyses that cover a range of (hopefully reasonable and nearly complete) possibilities for detection efficiencies. That approach is the one taken here.
3.3. Which Parameters to Study
In any analysis one must choose to study some subset of planetary and stellar parameters and to ignore others. For a robust fit, the PLDF should include all parameters that cause the detection efficiency to change significantly across the parameter space studied. Ignoring a parameter in the PLDF amounts to averaging the PLDF over the relevant parameter. However as expressed in Equation (13), the correct way to predict the yield of a survey is to weight the PLDF by the detection efficiencies before averaging.
A safer way to ignore a parameter (relevant to the detection efficiency) is to first fit it and then marginalize over it. With this caveat made explicit we now show how the general formalism treats the case of a PLDF that depends only on planetary parameters, , or only on stellar parameters, , (most famously stellar metallicity).
Planetary Parameters Only
For a shape function, , that is independent of stellar properties, the shape integral simplifies to
(23) 
where the stellaraveraged efficiencies are
(24) 
The fit procedure is simplified in this case. The stellar averaging is done a single time, and not for every choice of during optimization.
This case is the one most relevant to this work. Specifically in Section 4 we consider a shape function
(25) 
This powerlaw distribution in planetary radius and orbital period, identifies the planetary parameters as
(26) 
with the help of a reference radius and period , and defines the shape parameters as powerlaws. Appendix A gives an analytic solution for these best fit powerlaw exponents for an ideal transit survey with unity discovery efficiency and no (unvetted) false positives.
Stellar Parameters Only
We now consider fitting a PLDF that depends on stellar parameters only. As cautioned at the beginning of this section, it would be a bad idea to do this for a transit survey since the detection efficiencies always vary with orbital period, and also vary with planetary radius unless small planets with a discovery efficiency below unity are ignored. We consider this case for other survey types and for completeness.
Since the PLDF is independent of by assumption, we can integrate the differential distribution of Equation (12) over a volume, , of planetary parameter space. Defining gives
(27) 
which represents the NPPS with characteristics . The shape integral simplifies slightly to
(28)  
(29) 
by averaging the efficiencies over planetary parameters.
The best fit solutions are still given by Equations (17) and (19), and can be replaced with (no approximation needed) to remove all references to the differential distribution.
4. Powerlaw Fits to Kepler Data
We now use the Kepler planet candidates announced by BKep11 to constrain the underlying PLDF (planetary distribution function) as a powerlaw,
(30) 
in planet radius, , and orbital period . Equation (30) represents the NPPS (number of planets per star) per logarithmic interval in and . The actual PLDF can deviate from a powerlaw (and the data show that it does). Nevertheless, powerlaws are a common and useful approximation that captures basic trends in the data. Identifying deviations from powerlaw behavior reveals the relevant scales of planet formation and migration, and we show that Kepler data identify these trends extraordinarily well.
We study the full planet sample around Solar type stars in §4.1. In §4.2, we compare the distributions of shorter and longer period planets. Further comparing results for larger and smaller planets in §4.3 shows that the sample has a deficit of intermediatesized planets at the shortest orbital periods. For all these correlations, the S/N threshold has limited impact (§4.4).
4.1. The Full Sample
Our most complete sample of planets covers and . (See §2.3 for details.)
Fig. 3 plots planet counts of this full sample. The planet counts are binned by () in the left (right, respectively) panel. Raw counts are divided by the number of stars surveyed and the logarithmic bin size.
(31) 
where brackets indicate marginalization over and , respectively. Aside from sampling noise, these values are independent of the survey size or choice of bin width.
The best fit PLDF is “projected” into observational space in Fig. 3. This projection requires weighting by the detection efficiencies of transits, , and discovery, , before marginalizing over (or ), as in Equation (31). The comparison of the fits to histograms is somewhat misleading. The actual analysis uses unbinned data and fits and simultaneously, not separately to the size and period distributions. However the comparison of uncorrected counts to a projected PLDF does mimic the fitting procedure, which applies the efficiencies to the “theory” (i.e. powerlaw) not the detections.
We now explain how the powerlaw slopes are altered by the projection into observational space. Though the best fit period distribution in the right panel of Fig. 3 has , the curve in observational space has , smaller by about one. A decrease of precisely 2/3 is due to the transit efficiency. The remainder (here ) is due to the period dependence of and depends on the amount of small planets in the sample. From Equation (3), this correction could be as large as , the period exponent of at low efficiencies. Though not easily visible, there is modest curvature to the projected period powerlaw due to lower discovery efficiencies at longer periods.
The discovery efficiency has a much more dramatic effect on the projected size distribution, whose curve is evident in the left panel of Fig. 3. There are no free parameters to adjust either the position or angle of the break, because is fixed. Above , the Kepler discovery efficiency is unity for all periods considered. Thus the slope of the projected size distribution at large radii matches the bestfit . The break starts at , because drops below for these sizes, as shown in Fig. 2. At small , the powerlaw slope of the projected radial PLDF is , as fixed by the radial dependence of the discovery efficiency in Equation (3).
The thickness of the fits shown in Fig. 3 (and other figures) indicate deviations from maximum likelihood. Fig. 6 shows the corresponding error ellipse in the and parameters (the filled black oval is for the full planet sample). The 1 errors on and reported in Table 1 are onedimensional errors holding the other parameter fixed. As explained in Section 3, the normalization is not independently varied during the fit procedure, but follows from the values of , and the planettostar ratio. The errors on indicate the range of values taken inside the 1 error ellipse of and .
Our best fit size distribution has . This value indicates that smaller planets are much more abundant, but that larger planets contain most of the mass. At constant planet density, would give small planets most of the mass. However larger planets tend to have lower density, . Take as an example.
While planet counts drop towards small sizes, so does the projection of the bestfit powerlaw, as described above. The good qualitative agreement between the counts and the curve (in the left panel of Fig. 3) means that the drop in planet counts is mostly due to selection effects. The actual size distribution continues to rise towards smaller planets.
The best fit period law, , indicates that planets are more closely packed further from the star, for logarithmic intervals in or semimajor axis . For the planet density per linear interval in would increase.
However Fig. 3 shows that a single powerlaw is a qualitatively poor fit to the data due to a sharp drop in planet counts below days. Planet counts flatten at longer periods. As explained above, this means that the actual period distribution continues to rise towards longer periods. This break in the period distribution motivates our division of the planet sample below.
4.2. Short Vs. Long Periods
Since the period distribution deviates strongly from a simple powelaw, we divide the planets into short and long period samples. We use days as the dividing line. The choice gives us comparable numbers of planets in each sample and also gives an adequate range of periods over which to measure a powerlaw slope. Planet counts in different cuts are summarized in the column of Table 1.
Fig. 4 shows the planet counts and the results of independent powerlaw fits to the short and long period data. These fits are not broken powerlaws, as the values are not required to match at the period boundary. Two powerlaws vastly improve the qualitative fit to the data. More complicated functions could more precisely describe the location and nature of the period break, but are not considered here.
The radius distribution in the left panel includes overlapping histograms for the short and long period samples. The longer period (“slow”) sample has a steeper radius distribution than the shorter period (“fast”). The difference of has about 2.5 significance. By eye, the difference in the distribution is clearly driven by the flatness of the fast sample between —. To clarify this behavior we now examine the differences between small and large radius planet samples.
4.3. Quadrants
We subdivide the planets into four samples by splitting both the shorter and longer period populations (still defined relative to days) into “small” and “big” samples relative to . The results of the four independent powerlaw fits are plotted against planet counts in Fig. 5. Fig. 6 plots error ellipses to compare the slope changes and their statistical significance. The error ellipses of the subsamples are much larger than the full sample. Smaller number counts and the decreased leverage of the reduced range in and (over which to define a slope) both play a role.
The period distributions are only modestly affected by this split. At both short and long periods, the best fit of the small and big samples agree within the statistical uncertainties. Nevertheless some qualitatively interesting features appear in the right panel of Fig. 5. The small planet counts begin their decline below about 7 days (as noted by HKep11). The large planet sample remains flat towards shorter periods — with evidence of a peak at days. This peak is responsible for the flat appearance of the overall period distribution down to 3 days, as seen in the right panels of Figs. 3 and 4.
The behavior of the size distributions is even more intriguing. For big planets, the best fit radius law changes by , with almost significance. Since the discovery efficiency of Kepler is near unity for these large planets, no obvious systematic uncertainties should be at work. For small planets, the change in the size distributions, . While of more modest significance, this change is of the opposite sign.
Fig. 7 plots the best fit size distributions for each quadrant, with the full planet sample shown for comparison. At longer periods the size distribution flattens towards smaller radii, i.e. is less negative for the small planets. For the shortest period planets the behavior is precisely opposite. The size distribution is quite shallow at large sizes and steepens for the small radius sample.
The implications of this behavior are intriguing. The processes that form and/or migrate planets inside days clearly disfavors planets of a few to several . A plausible scenario (but lacking in specifics) follows. At short periods, only the most massive giants — true Jovians with can efficiently retain their atmospheres. Lesser giants are stripped of their atmospheres, including ices that sublimate from the surface. The demotion of ice and gas giants to small rocky cores can simultaneously explain the flatter size distribution at large sizes and the steeper distribution at small sizes.
If this hypothesis is born out by further observations and theoretical study, it would provide yet another pillar of support for the coreaccretion hypothesis and for migration as the source of short period planets.
4.4. Systematic Uncertainty
To test the quality of the fits, we consider an increased S/N threshold. Increasing this threshold preferentially removes smaller planets from the sample. In principle the appropriate S/N threshold is set by the analysis of the discovery efficiency, . As discussed in §2.3, increasing the S/N threshold is a conservative way to check how the assumed impacts the results.
Fig. 3 includes dotted histograms that correspond to the higher S/N sample. Planet counts are noticeably lower below and across all periods.
Table 1 includes fits to the higher S/N planet samples. All of the qualitative trends discussed in this section persist in the higher S/N sample and with only mostly reduced significance. Indeed, the interesting trends we identified for small planets was a steepening of the size distribution at small periods and a flattening at longer periods. It would be rather difficult for systematic uncertainties to conspire to produce both trends.
Despite this robustness, is crucial for the statistical analysis of small planets. The final rows of Table 1 shows that setting (but still including the transit probability) gives very discrepant results. These discrepancies are expected due to the low discovery efficiency below , as shown in Fig. 2.
5. Planet Occurrence
To derive the number of planets per star, NPPS, we consider our nonparametric analysis. For and days, we find planets per star, consistent with HKep11.
When we consider smaller planets the NPPS is roughly one. For , S/N thresholds of 10, 15 and 20 give , 1.05 and 0.72 planets per star. Since the S/N = 20 threshold is conservative, it seems likely .
The likely existence of more than one planet per Solar type star, especially within 50 days, is remarkable. However it does not imply that the Sun is exceptional for lacking a planet so close. For there are 2—3 planets per planetary system (Lissauer et al., 2011b). Multiplicity rates will almost certainly increase when planets down to are considered.
It is inherently speculative to extrapolate into unobserved regions of parameter space. Nevertheless, all the powerlaw fits of the previous section have distributions that rise with period. To emphasize the implications of this rise, we calculate the extrapolated density of Earthsized planets at 1 AU. The final column of Table 1 expresses the normalization to the fits as
(32) 
the NPPS in a logarithmic interval centered on an Earth radius and a year.
The most relevant fit is to the small planet, long period sample, which gives . Integrating this powerlaw gives on average Earthlike planets with periods below a year. If the distribution function of planets does not turn over shortly beyond periods of 50 days, the implications for habitable Earths around Sunlike stars is truly staggering.
6. Concluding Discussion
We develop a general technique for the analysis of exoplanet survey data. The merits of the approach, a generalization of the TT02 maximum likelihood analysis, include the following:

Data are analyzed without binning in order to preserve the statistical significance of each detection.

Describing the number of planets per star (NPPS) allows multiplanet systems to be included.

Selection effects are used to calculate the number of expected detections. They are not used to enhance (or diminish) the number of actual detections, which affects error estimates.

The overall normalization is analytically removed from the likelihood function. With one less parameter in the numerical optimization, the fit is simpler and possibly more robust.

Different surveys can be jointly analyzed, as shown by TT02. Even different survey methods that probe different regions of parameter space can be combined if selection effects are well characterized.
We apply the technique to Kepler planet candidates released by BKep11. We focus on planets orbiting Solar type stars for which HKep11 has quantified the Kepler discovery efficiency. Fits to these efficiences (see Fig. 2) reveal a remarkably smooth variation with planet radius and orbital period that agrees with analytic predictions for transit surveys. This regular behavior gives us confidence to analyze detections down to .
From this analysis, the best estimate of the number of planets per star bigger than and with periods below is — . Uncertainty in the discovery efficiencies dominate these estimates. Since planet occurrence rises towards both small sizes and longer periods, the only question is when (not if) we can claim with certainty that there is more than one planet per star.
The shape of the radius and period distributions is even more informative than absolute number counts. The most notable trend is a difference between the size distributions of shorter and longer period planets (below and above 7 days), plotted in Fig. 7. While the shorter period planets are less abundant overall, their relative deficit is most pronounced around . We interpret this finding as arising from the preferential evaparation of ice and lower mass gas giants that migrated too close to their host star. In this picture, the steepness of the size distribution of small planets at short periods arises from the remnant cores of these stripped giants.
If this interpretation holds, it would add support to the core accretion hypothesis (Pollack et al., 1996), though this theory is not seriously challenged at masses below and possibly much higher (Kratter et al., 2010). In particular, cores would have to be present by the end of migration. This constraint applies to models where the cores sediment over time (Helled & Schubert, 2008).
We do not here attempt to constrain the mode of migration with these size trends. Possibilities include smooth inward migration through the disk, a theory that was highly developed even before the discovery of extrasolar planets (Lin & Papaloizou, 1986; Ward, 1986; Artymowicz, 1993). Planetplanet scattering (Chatterjee et al., 2008) and Kozai oscillations (Wu & Murray, 2003) have also been proposed as mechaisms to deliver planets close to their hosts. Recently Wu & Lithwick (2010), proposed a new mechanism — secular chaos — and compared all the proposed mechanisms to observations of hot Jupiters. These theories and others should now be compared to Kepler data for low mass planets as well.
Appendix A Analytic Solutions for Powerlaw Fits
This appendix derives an analytic solution to the best fit powerlaws of the planetary radius and orbital period. The transit probability is included, but all other detection efficiencies are set to unity, as appropriate for a perfect transit survey. Since numerical fits are more flexible, this derivation is intended to provide insight and to provide a check on numerical solutions (when only the transit probability is included as a selection effect).
Consider the detection of planets in a transit survey. We fit the PLDF to the powerlaw of Equation (25) (equivalent to Equation (30) used in our main analysis of Kepler data). The logarithms of radius and period are defined as the parameters and as in Equation (26). We define the reference planet radius, , and orbital period so the (here rectangular) domain of the survey is centered on and all detections (indexed by ) fall within
(A1) 
The transit efficiency of Equation (1) can be written as
(A2) 
and we ignore other selection effects.
The number of expected transits for such a perfect survey is with
(A3) 
given by Equation (23) (which averages the general Equation (13) over stellar parameters).
The bestfitting and follow from Equation (19) as
(A4) 
where the obsevational mean , and similarly for . The relevant transcendental function, is monotonic in with as and for . The closed form solutions of Equation (A4) have several features that relate to more complete analyses:

The powerlaws increase monotonically with the average value of the (log of the) observables, independent of how the values are distributed around that mean. Higher order shape functions can describe more detailed features in the data.

The period powerlaw in the PLDF is larger than the period law describing the detections. The transit efficiency, causes this increase. In general any powerlaw selection effect adjusts from the observed to the underlying distribution this way.

The overall magnitude of the detection efficiency does not affect powerlaws, or shape parameters in general. For instance, uncertainty in due to uncertainty in stellar densities does not affect the shape function. The overall planet occurrence of course does depend on the magnitude of detection efficiencies, via the normalization given by Equation (17).

Solutions for the best fit powerlaws decouple. In general, the solution for two shape parameters decouple if the shape parameters (and the physical parameters they describe) are separable in the shape function and the relevant physical parameters are separable in the detection efficiencies. This separability condition is not met for out full treatment of the Kepler data. While the efficiency ratio of Equation (3) is separable, the efficiency itself, , is not. This leads to the modest covariance between and seen in Fig. 6.
Planet Sample  Fit Properties  

Tag  Radii  Periods 

planets 

[] 
[days]  (corrected)  per star  (radius law)  (period law)  
Full  0.5 — 20  0.5 — 50  562  2.38  
Full/  2 — 20  0.5 — 50  372  0.21  
(Short vs. Long Period: “Fast” vs. “Slow”)  
Fast  0.5 — 20  0.5 — 7  211  0.18  E2  
Slow  0.5 — 20  7 — 50  351  2.26  
Fast/  2 — 20  0.5 — 7  111  0.02  E2  
Slow/  2 — 20  7 — 50  261  0.17  
Small  0.5 — 3  0.5 — 50  389  2.30  
Big  3 — 20  0.5 — 50  173  0.08  
(Quadrants)  
Small/Fast  0.5 — 3  0.5 — 7  148  0.25  E2  
Big/Fast  3 — 20  0.5 — 7  63  0.01  )E2  
Small/Slow  0.5 — 3  7 — 50  241  1.31  
Big/Slow  3 — 20  7 — 50  110  0.06  
(S/N20)  
Full/sn20  0.5 — 20  0.5 — 50  453  1.37  
SmFa/sn20  0.5 — 3  0.5 — 7  105  0.11  E2  
BiFa/sn20  3 — 20  0.5 — 7  62  0.01  E2  
SmSl/sn20  0.5 — 3  7 — 50  178  0.62  
BiSl/sn20  3 — 20  7 — 50  108  0.06  
(TRANSIT PROBABLILITY ONLY, Perfect Discovery Efficiency, )  
Full/tp  0.5 — 20  0.5 — 50  562  0.24  
SmFa/tp  0.5 — 3  0.5 — 7  148  0.03  
BiFa/tp  3 — 20  0.5 — 7  63  0.01  E2  
SmSl/tp  0.5 — 3  7 — 50  241  0.12  
BiSl/tp  3 — 20  7 — 50  110  0.06 
Footnotes
 slugcomment: Draft Modified March 6, 2018
 For an instructive explanation see http://oklo.org/2011/03/23/theeccentricitydistribution/
 This step overcomes the rounding of the reported to intervals, and is merely for convenience in defining boundaries in an assumptionfree way.
 Nondetections however can enter via upper limits that define the efficiency, see §2.
 For brevity the exponent in the phase space volume is dropped, i.e. here and throughout. We keep this exponent in the denominator of derivatives to be explicit that these are not gradients, as in Equation (11).
 E.g. divided by for a bin between and and similarly for period bins.
 Appropriate if 1 planets have and 10 planets have .
 :Number of Kepler planet candidates with S/N 10 (or 20 where indicated) as reported in BKep11 .
 See §5 for more accurate nonparametric fits to the planet occurrence.
 Planet radii are computed form the products of and given in Tables 1 and 2 of BKep11.
References
 Artymowicz, P. 1993, ApJ, 419, 166
 Borucki, W. J., et al. 2010, Science, 327, 977
 —. 2011, ArXiv eprints
 Burke, C. J. 2008, ApJ, 679, 1566
 Charbonneau, D., et al. 2009, Nature, 462, 891
 Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580
 Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
 Cumming, A., Butler, R. P., Marcy, G. W., Vogt, S. S., Wright, J. T., & Fischer, D. A. 2008, PASP, 120, 531
 Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102
 Gaudi, S. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 366, Transiting Extrapolar Planets Workshop, ed. C. Afonso, D. Weldrake, & T. Henning, 273–+
 Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549
 Gonzalez, G. 1997, MNRAS, 285, 403
 Helled, R., & Schubert, G. 2008, Icarus, 198, 156
 Howard, A. W., et al. 2010, Science, 330, 653
 —. 2011, ArXiv eprints
 Johansen, A., Youdin, A., & Mac Low, M. 2009, ApJ, 704, L75
 Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905
 Kenyon, S. J., & Bromley, B. C. 2006, AJ, 131, 1837
 Kratter, K. M., MurrayClay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375
 Latham, D. W., et al. 2011, ArXiv eprints
 Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846
 Lissauer, J. J., et al. 2011a, Nature, 470, 53
 —. 2011b, ArXiv eprints
 Marcy, G. W., & Butler, R. P. 1996, ApJ, 464, L147+
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355
 Moorhead, A. V., et al. 2011, ArXiv eprints
 Morton, T. D., & Johnson, J. A. 2011, ArXiv eprints
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62
 Ragozzine, D., & Holman, M. J. 2010, ArXiv eprints
 Tabachnik, S., & Tremaine, S. 2002, MNRAS, 335, 151
 Tingley, B., Bonomo, A. S., & Deeg, H. J. 2011, ApJ, 726, 112
 Ward, W. R. 1986, Icarus, 67, 164
 Wu, Y., & Lithwick, Y. 2010, ArXiv eprints
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605
 Youdin, A. N. 2010, in EAS Publications Series, Vol. 41, EAS Publications Series, ed. T. Montmerle, D. Ehrenreich, & A.M. Lagrange, 187–207
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
 Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494