A CONTINUUM OF PLANET FORMATION BETWEEN 1 AND 4 EARTH RADII

The probability that a giant planet orbits a star is a steeply rising function of the host star's metallicity (e.g., Santos et al. 2004; Fischer & Valenti 2005). This observation is the key piece of evidence that the giant planets identified by the radial velocity and transit techniques form through core accretion and not through gravitational instability. This observation is perhaps the most important constraint placed on models of planet formation since the discovery of the first exoplanets.

The connection between stellar metallicity and the presence of small planets is less clear. The Neptune-mass planets discovered by radial velocity surveys do not appear to preferentially orbit metal-rich FGK stars (e.g., Sousa et al. 2008; Mayor et al. 2011). While Kepler has discovered a large number of small exoplanet candidates (planets from here), it has not yet settled the issue. Schlaufman & Laughlin (2011) showed that while the giant planets discovered by Kepler orbit metal-rich stars, the small planets discovered around F and G stars did not appear to prefer metal-rich stars. This observation was later confirmed by Buchhave et al. (2012).

Recently, Buchhave et al. (2014, B14 hereafter) argued that the observed distribution of metallicity in a sample of more than 400 Kepler planet host stars revealed 3 distinct clusters of exoplanets: terrestrial planets with planet radius R_p ≲ 1.7 R_⨁, "gas dwarf" planets with 1.7 R_⨁ ≲ R_p ≲ 3.9 R_⨁, and ice or gas giants with R_p ≳ 3.9 R_⨁. They suggested that these three populations formed via distinct planet formation channels.

To reach that conclusion, B14 repeatedly split their sample of planet host star metallicity measurements into small-planet and large-planet subsamples for different choices of the radius boundary dividing the two subsamples. They calculated the p-value from a two-sample Kolmogorov–Smirnov test on the two subsamples as a function of planet radius and identified local minima p-values. They saved the radii at which the local minima occurred. To account for measurement uncertainties, they repeated this process 10⁶ times, sampling the planet radius and host star metallicity from their uncertainty distributions on each iteration. They identified a distinct p-value minimum at R_p = 1.7 R_⨁ and argued that it represents a boundary between terrestrial and "gas dwarf" planets. This approach is inappropriate because it performs a large number of hypothesis tests on the same data set without correcting the test thresholds to account for the large number of tests. That strategy is known to produce a high false-discovery rate (e.g., Dunn 1959, 1961). Moreover, the B14 technique creates a sequence of p-values at many split points for data subject to measurement uncertainty. Consequently, before attaching any significance to features in that sequence of p-values, it is also critical to ensure that the p-values that result from the Monte Carlo simulation accurately represent the p-value measurement uncertainties that result from uncertainties in the input sample.

There are at least four more problems with the analysis presented in B14. First, B14 overlooked the effect of planet radius uncertainty due to transit depth uncertainty. Second, their approach used an asymptotically inconsistent estimator of the average p-value at each split point in the presence of observational uncertainty. Third, their analysis is subject to the multiple comparisons problem, which reduces the significance of their observation by a large amount. Fourth, while B14 assert that local minima in a plot of p-value as a function of split radius indicate transitions between distinct clusters of exoplanets, this is not necessarily so. I describe my sample selection in Section 2, I detail each issue with the B14 calculation in Section 3, I outline a more rigorous way to investigate the issue in Section 4, and I discuss the implications and my conclusion in Section 5.

I use the planet host star data from B14. Those data include ${{T}_{{\rm eff}}},{\rm log} g$, [M/H], M_*, R_*, and their associated uncertainties. B14 did not include in their planet radius uncertainties the effect of uncertainties in transit depth, even though transit depth uncertainties are more important than the host star radius uncertainties in 25% of the sample. As a result, I supplement the B14 data with the latest Kepler object of interest period and R_p/R_* estimates from the Kepler CasJobs database² hosted by the Mikulski Archive for Space Telescopes. I then recompute planet radii from the B14 stellar radii and the updated transit depths. Following B14, I remove from the sample all planets smaller than 3 R_⨁ subject to strong stellar irradiation (i.e., F_ν > 5 × 10⁵ J s⁻¹ m⁻²), as these planets may have lost a significant fraction of their initial atmospheres. I plot these data in Figure 1.

Zoom In Zoom Out Reset image size

3.1. An Asymptotically Inconsistent p-value Estimator

An asymptotically inconsistent estimator of a parameter does not converge to the true value of the parameter in the large-sample limit. One problem with the B14 analysis is that they used an asymptotically inconsistent estimator of the p-value averaged over planet radius measurement uncertainty in their Monte Carlo simulation. The p-value measurements depend on the planet radius measurements, which are subject to measurement uncertainty in the inferred stellar radii R_* and the measured ratios R_p/R_*. The true p-values in the absence of uncertainty cannot be measured directly. Instead, one can only measure p'

$$\begin{eqnarray}&&p^{\prime} =p+N\left( 0,{{\sigma }^{2}} \right),\end{eqnarray} \tag{ 1 }$$

where p is the true p-value and $N(0,{{\sigma }^{2}})$ is due to measurement uncertainties in R_* and R_p/R_*. Repeatedly calculating p' after perturbing each planet radius due to the uncertainties in R_* and R_p/R_* and averaging the result will provide an asymptotically consistent estimate of p by the central limit theorem

$$\begin{eqnarray}&&E\left[ p^{\prime} \right]=E\left[ p+N(0,{{\sigma }^{2}}) \right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&E\left[ p^{\prime} \right]=E\left[ p \right]+E\left[ N(0,{{\sigma }^{2}}) \right],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}\;p^{\prime} =\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}\;p+\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}\;N(0,{{\sigma }^{2}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\overline{p^{\prime} }=\bar{p}+0\Rightarrow \overline{p^{\prime} }=\bar{p}.\end{eqnarray} \tag{ 5 }$$

B14 never averaged the p-value produced for each split point over all iterations. Instead, after each iteration of their Monte Carlo simulation they identified the local p-value minima and saved them. After completing 10⁶ Monte Carlo iterations, they determined the mean radii at which local p-value minima occurred by averaging over the individual radii calculated on each iteration. In other words, they applied the nonlinear function f that takes a sequence of p-values and identifies the radii of local p-value minima before averaging over all iterations to identify the mean radii at which p-value minima occur. The central limit theorem does not apply in this case, as

$$\begin{eqnarray}&&f(p^{\prime} )=f\left( p+N(0,{{\sigma }^{2}}) \right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&E[f(p^{\prime} )]=E\left[ f\left( p+N(0,{{\sigma }^{2}}) \right) \right],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}f(p^{\prime} )=\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}f\left( p+N(0,{{\sigma }^{2}}) \right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\overline{f(p^{\prime} )}=\overline{f\left( p+N(0,{{\sigma }^{2}}) \right)}{\nRightarrow }\overline{f(p^{\prime} )}=\overline{f(p)}.\end{eqnarray} \tag{ 9 }$$

As a result, the p-values in their Figure 1 improperly account for measurement uncertainty and are asymptotically inconsistent with the true p-values absent measurement uncertainty.

To address that problem, I first generate 10⁵ realizations of each planet radius from the distributions that result from the propagation of measurement uncertainties in R_* and R_p/R_*. I split the metallicity data into small-planet and large-planet subsamples at 321 split points from 0.3 to 13.1 R_⨁ in steps of 0.04 R_⨁ and compute the p-value from a two-sample Kolmogorov–Smirnov test. I save the resulting p-value for each split point and repeat this process 10⁵ times. At the end of the calculation, I average the p-values for each split point. I plot the result in Figure 2. The apparent local minimum in the p-value distribution identified by B14 at R_p = 1.7 R_⨁ is not present.

Zoom In Zoom Out Reset image size

3.2. The Multiple Comparisons Problem

3.3. Do Local p-value Minima Indicate Distinct Exoplanet Regimes?

While B14 argue that local minima in a plot of split radius versus p-value indicate transitions between distinct exoplanet clusters, this is not always the case. To demonstrate this, I use the same Monte Carlo simulation described in Section 3.1. However, instead of using the observed metallicities, on each iteration I randomly sample the metallicities of stars hosting planets with R_p ≤ 1.7 R_⨁, 1.7 R_⨁ < R_p ≤ 3.9 R_⨁, and R_p ≥ 3.9 R_⨁ from their observed distributions. Those distributions are N(0.00, 0.20²), N(0.05, 0.19²), and N(0.18, 0.19²). I plot the result in Figure 3. Despite the fact that distinct metallicity distributions were imposed on each planet cluster, there is no local minimum in the p-value distribution at the boundary between the terrestrial and "gas dwarf" planets. The inability of the B14 technique to identify a metallicity boundary imposed by construction as a local p-value minimum implies that the technique is not sensitive to subtle features in the metallicity distribution.

Zoom In Zoom Out Reset image size

A better way to identify the number of subpopulations required by the planet radius and host star metallicity data is to compare statistical models with varying numbers of components, then identify the model that has the minimum number of parameters yet the maximum likelihood of producing the observed data. I consider two classes of models. First, I fit single-population linear models of the form

$$\begin{eqnarray}&&[M/H]={{a}_{0}}+\mathop{\sum }\limits_{j=1}^{m}{{a}_{j}}R_{p}^{j}+\epsilon ,\end{eqnarray} \tag{ 10 }$$

where is the standard uncertainty term in the regression equation. Second, I fit finite Gaussian mixture models with varying numbers of subpopulations of the form

$$\begin{eqnarray}&&{{\mathop{\prod }\limits_{i=1}^{n}}^{{}}}\mathop{\sum }\limits_{j=1}^{m}{{w}_{j}}{{N}_{j}}({{{\bf x}}_{i}}\mid {{{\boldsymbol{ \mu }} }_{j}},{{{\bf \Sigma }}_{j}}),\end{eqnarray} \tag{ 11 }$$

where x is the data, m is the number of components in the model, the w_j are weights such that $\mathop{\sum }\limits_{j=1}^{m}{{w}_{j}}=1,$ and each N_j is a two-dimensional Gaussian component of the overall density

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{j}}({\bf x}\mid {{{\boldsymbol{ \mu }} }_{j}},{{{\Sigma }}_{j}}) & = & \frac{1}{2\pi \mid {{{\bf \Sigma }}_{j}}{{\mid }^{1/2}}} \\ {} & {} & \times {\rm exp} \left\{ -\frac{1}{2}{{({{{\bf x}}_{i}}-{{{\boldsymbol{ \mu }} }_{j}})}^{T}}{\bf \Sigma }_{j}^{-1}({{{\bf x}}_{i}}-{{{\boldsymbol{ \mu }} }_{j}}) \right\}. \\ \end{array}\end{eqnarray} \tag{ 12 }$$

Here ${{{\boldsymbol{ \mu }} }_{j}}$ and ${{{\bf \Sigma }}_{j}}$ are the mean and covariance of each of the m components of the model. I fit the Gaussian mixture models using the mclust ³ package in R ⁴ (Fraley & Raftery 2002; Fraley et al. 2012; R Core Team 2014).

To account for the observational uncertainties, I use a Monte Carlo simulation. I sample the planet radii from the distributions that result from the propagation of measurement uncertainties in R_* and R_p/R_* and directly use the measured metallicities (since the uncertainty in [M/H] is already reflected in the uncertainty in R_*). On each iteration, I fit linear models of the form of Equation (10) for m = 1, 2, ..., 5 and Gaussian mixture models m = 1, 2, ..., 7. I choose both the best linear and Gaussian mixture models using the Bayesian information criterion (BIC; Schwarz 1978), then use the Akaike information criterion (Akaike 1974) to choose between the favored linear and Gaussian mixture models. I repeat this process 10³ times. In all cases, the best linear model is preferred. For the linear model, the m = 1 model is favored 76.4% of the time, the m = 2 model is favored 22.3% of the time, while a higher-order model is favored 1.3% of the time. While the Gaussian mixture model is disfavored relative to the linear model, the two-component model is the best of the mixture models: the two-component model is preferred on 92.8% of the iterations, while the three-component model is preferred 7.2% of the time. I plot representative examples of the BIC-selected models from one iteration of my Monte Carlo simulation in Figure 4.

Zoom In Zoom Out Reset image size

The performance of a large number of tests on the same data set without correcting the test thresholds, the use of an asymptotically inconsistent estimator of the p-value in the presence of measurement uncertainty, the oversight of the multiple comparison problem, or the inability of the B14 technique to identify an imposed metallicity effect as a p-value minimum are all sufficient reasons to be skeptical of the claimed transition at R_p = 1.7 R_⨁. The problems with the B14 analysis technique cannot be mitigated by examining a larger or independent data set—they are inherent in the analysis technique itself. Instead, the analysis in Section 4 shows that a smooth, one-component linear model is a better fit to the data than any multi-component model. If a multi-component model is used, then the two-component model is consistently a better fit than the three-component model. As a result, the planet radius and metallicity data for the Kepler F and G star planet hosts does not support the idea of multiple types of small planets. Instead, a continuum of planet sizes between 1 R_⨁ and 4 R_⨁ are likely formed independent of the amount of solids present.

While observational evidence suggests that most planets larger than about 2 R_⨁ have significant hydrogen atmospheres (e.g., Marcy et al. 2014), Kepler-10c is an exception with a radius of 2.35 R_⨁ and a density of 7.1 g cm⁻³ (Dumusque et al. 2014). Likewise, smaller planets probably have a wide range of atmospheric properties (e.g., Rogers 2014; Wolfgang & Lopez 2014). Moreover, a wide range of densities can be present even in the same system, with Kepler-36 the best example (Carter et al. 2012). For these reasons, near solar metallicity it does not seem likely that the final masses or compositions of small exoplanets are controlled primarily by the amount of solid material present in their parent protoplanetary disks.

I thank Lars Buchhave, Andy Casey, Bryce Croll, David W. Latham, Dimitar Sasselov, and Josh Winn. I am especially grateful to the referee, Eric Feigelson, for suggestions that substantially improved this paper. This research has made use of NASA's Astrophysics Data System Bibliographic Services. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-Hubble Space Telescope data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Support for this work was provided by the MIT Kavli Institute for Astrophysics and Space Research through a Kavli Postdoctoral Fellowship.

Facilities: Kepler