Examining the Radius Valley: a Machine Learning Approach

TL;DR

This study confirms the existence of the radius valley in exoplanets using a large dataset and characterizes it as a negative power law slope, supporting photoevaporation theories over gas-poor formation models.

Contribution

It provides a comprehensive analysis of the radius valley with a larger sample, confirming its power-law nature and quantifying its slope across multiple methods.

Findings

Radius valley confirmed with larger sample

Slope of the valley is approximately -0.319

Results support photoevaporation as the main process

Abstract

The ''radius valley" is a relative dearth of planets between two potential populations of exoplanets, super-Earths and mini-Neptunes. This feature appears in examining the distribution of planetary radii, but has only ever been characterized on small samples. The valley could be a result of photoevaporation, which has been predicted in numerous theoretical models, or a result of other processes. Here, we investigate the relationship between planetary radius and orbital period through 2-dimensional kernel density estimator and various clustering methods, using all known super-Earths ($R<4.0R_E$). With our larger sample, we confirm the radius valley and characterize it as a power law. Using a variety of methods, we find a range of slopes that are consistent with each other and distinctly negative. We average over these results and find the slope to be $m=-0.319^{+0.088}_{-0.116}$. We repeat our analysis on samples from previous studies. For all methods we use, the resulting line has a negative slope, which is consistent with models of photoevaporation and core-powered mass loss but inconsistent with planets forming in a gas-poor disk.

Tables

Table 1. Table 1: The number of planets and average (median) uncertainty in planetary radius for planets in each of the data sets we use. For data set ”all,” we pull all exoplanets from the exoplanet archive (April 5, 2019). We only use planets with published uncertainties in both planetary radius and in orbital period. We also limit all samples to super-Earths ( R < 4.0 ​ R E 𝑅 4.0 subscript 𝑅 𝐸 R<4.0R_{E} ) with orbital periods P < 50 𝑃 50 P<50 days.

Data set	Sample Size	($R_E$)	Ref
All	2066	0.252	–
F17	766	0.19	Fulton et al. (2017)
V18	81	0.0615	Van Eylen et al. (2018)
F18	1334	0.089	Fulton & Petigura (2018)

Table 2. Table 2: Resulting slope measurements and their uncertainties for our different methods and data sets. Here, the ”all” data set refers to all known exoplanets with orbital period and planetary radius estimates and uncertainties and P < 50 𝑃 50 P<50 days and R P < 4.0 ​ R E subscript 𝑅 𝑃 4.0 subscript 𝑅 𝐸 R_{P}<4.0R_{E} , σ r subscript 𝜎 𝑟 \sigma_{r} refers to data sets where we accounted for uncertainties in the planetary radii by drawing 1000 planets with replacement from the data set and drawing each planet’s radius from a normal distribution centered on the nominal value with standard deviation of the uncertainty in the measurement, and P < 25 𝑃 25 P<25 refers to adding a cut to orbital period at 25 days to help mitigate the effects of incompleteness. The other samples refer to data sets used in F17: Fulton et al. ( 2017 ) , V18: Van Eylen et al. ( 2018 ) , and F18: Fulton & Petigura ( 2018 ) . See Section  2 for a description of each data set. ”DBSCAN” is density-based clustering, ”HIER” is hierarchical clustering cut at two branches, ”KDE” is 2D kernel density estimation, ”CLUST” is the results from all clustering methods (excluding KDE), and ”ALL M” is the results from all methods. See Section  3 for a description and comparison of all methods. All measurements from all methods and datasets are consistent with a negative slope. We take our bestfit slope to be the result of ”all” data set and all methods, in bold in the table above, m = − 0.319 − 0.116 + 0.088 𝑚 subscript superscript 0.319 0.088 0.116 m=-0.319^{+0.088}_{-0.116} .

Data set	Num. points	Kmeans	Kmedians	DBSCAN	HIER	KDE	CLUST	ALL M
All	2066	$-{0.329}_{-0.039}^{+0.037}$	$-{0.306}_{-0.031}^{+0.029}$	$-{0.324}_{-0.087}^{+0.085}$	$-{0.253}_{-0.307}^{+0.127}$	$-{0.356}_{-0.144}^{+0.146}$	$-{0.316}_{-0.084}^{+0.077}$	-0.319 ${}_{-0.116}{}^{}{}_{}{}^{+0.088}$
All, ${\sigma }_r$	1000	$-{0.590}_{-0.331}^{+0.181}$	$-{0.651}_{-0.282}^{+0.160}$	$-{0.658}_{-0.587}^{+0.332}$	$-{0.366}_{-0.636}^{+0.193}$	$-{0.333}_{-0.185}^{+0.166}$	$-{0.608}_{-0.371}^{+0.291}$	$-{0.531}_{-0.401}^{+0.265}$
All P	1855	$-{0.398}_{-0.050}^{+0.042}$	$-{0.366}_{-0.038}^{+0.036}$	$-{0.530}_{-0.177}^{+0.120}$	$-{0.221}_{-0.373}^{+0.104}$	$-{0.403}_{-0.059}^{+0.090}$	$-{0.393}_{-0.159}^{+0.104}$	$-{0.397}_{-0.123}^{+0.096}$
All P, ${\sigma }_r$	1000	$-{0.793}_{-0.243}^{+0.304}$	$-{0.745}_{-0.155}^{+0.191}$	$-{0.725}_{-0.276}^{+0.198}$	$-{0.603}_{-0.684}^{+0.325}$	$-{0.340}_{-0.106}^{+0.074}$	$-{0.738}_{-0.288}^{+0.281}$	$-{0.640}_{-0.352}^{+0.293}$
F17	766	$-{0.525}_{-0.074}^{+0.059}$	$-{0.450}_{-0.066}^{+0.042}$	$-{0.480}_{-0.143}^{+0.197}$	$-{0.783}_{-0.237}^{+0.378}$	$-{0.420}_{-0.020}^{+0.019}$	$-{0.508}_{-0.227}^{+0.116}$	$-{0.474}_{-0.183}^{+0.076}$
F17, ${\sigma }_r$	1000	$-{0.379}_{-0.078}^{+0.058}$	$-{0.348}_{-0.045}^{+0.032}$	$-{0.146}_{-0.068}^{+0.056}$	$-{0.236}_{-0.607}^{+0.172}$	$-{0.378}_{-0.032}^{+0.023}$	$-{0.326}_{-0.120}^{+0.202}$	$-{0.348}_{-0.086}^{+0.205}$
V18	81	$-{0.561}_{-0.196}^{+0.144}$	$-{0.553}_{-0.025}^{+0.025}$	$-{0.369}_{-0.120}^{+0.066}$	$-{0.272}_{-0.018}^{+0.018}$	$-{0.315}_{-0.042}^{+0.026}$	$-{0.476}_{-0.109}^{+0.204}$	$-{0.376}_{-0.194}^{+0.100}$
V18, ${\sigma }_r$	1000	$-{0.565}_{-0.039}^{+0.039}$	$-{0.570}_{-0.026}^{+0.031}$	$-{0.362}_{-0.073}^{+0.061}$	$-{0.226}_{-0.430}^{+0.021}$	$-{0.328}_{-0.015}^{+0.013}$	$-{0.535}_{-0.066}^{+0.292}$	$-{0.433}_{-0.160}^{+0.141}$
F18	1334	$-{0.370}_{-0.058}^{+0.038}$	$-{0.334}_{-0.030}^{+0.035}$	$-{0.303}_{-0.101}^{+0.086}$	$-{0.267}_{-0.325}^{+0.167}$	$-{0.372}_{-0.038}^{+0.031}$	$-{0.339}_{-0.088}^{+0.117}$	$-{0.351}_{-0.066}^{+0.099}$
F18, ${\sigma }_r$	1000	$-{0.417}_{-0.061}^{+0.081}$	$-{0.301}_{-0.060}^{+0.032}$	$-{0.376}_{-0.439}^{+0.112}$	$-{0.177}_{-0.328}^{+0.107}$	$-{0.345}_{-0.047}^{+0.025}$	$-{0.346}_{-0.147}^{+0.117}$	$-{0.346}_{-0.127}^{+0.089}$

Equations

$$\text{log}{R}_p=m\cdot \text{log}P+b$$

$$\hat{f}(x;H)=n^{-1}\sum _{i=1}^n{K}_H(x-X_i)$$

$$d_{(ij)k}={\alpha }_i{d}_{ik}+{\alpha }_j{d}_{jk}+\beta d_{ij}+\gamma |d_{ik}-d_{jk}|$$

$$\begin{array}{c}\text{arg min }\\ \bar{S}\end{array}\sum _{i=1}^K\sum _{\bar{x}\in S_i}||\bar{x}-{\mu }_i|{|}^2=\begin{array}{c}\text{arg min}\\ \bar{S}\end{array}\sum _{i=1}^K|S_i|\text{ Var }{S}_i$$

Full text

††pubyear: 2019††pagerange: Examining the Radius Valley: a Machine Learning Approach–References

1 Introduction

With the launch of Kepler in 2009 and the following years of analysis, the astronomical community witnessed the number of known exoplanets skyrocket from a couple hundred to a couple thousand (Lissauer et al., 2011, 2014; Fabrycky et al., 2014). Newfound discoveries such as compact systems (Lissauer et al., 2011; Muirhead et al., 2015), resonant chains (Lithwick & Wu, 2012; Fabrycky et al., 2014), and low density planets (e.g., Welsh et al., 2015; Weiss & Marcy, 2014; Ford, 2017) led to numerous studies of the exoplanet population as a whole, as well as countless studies focused on just one or a few systems. Many of these studies aimed to better characterize exoplanets and estimate the planet parameters to gain insight into planet formation.

Several theoretical models predict that the size of planets can be strongly influenced by the host stars. In particular, formation models find that atmospheric erosion, which mostly affects planets near their star, will result in a so-called ”photoevaporation valley,” a gap in the distribution of planetary radii around $2R_{E}$ (Owen & Wu, 2013; Jin et al., 2014; Lopez & Fortney, 2014; Chen & Rogers, 2016; Lopez & Rice, 2016; Owen & Wu, 2017). This valley is established by the boundary between planets that are massive enough to keep their gas envelopes and planets whose atmospheres are stripped away. The specific slope of the valley depends on various factors, such as the composition of the planets and the physics of the evaporation (e.g., Owen & Wu, 2017), but the overall sign of the slope is negative. In particular, the separation appears linear in a log-log plot of planetary radius vs. orbital period:

[TABLE]

where $P$ is the orbital period, $R_{p}$ is the planetary radius, $m$ is the slope in log-log space, and $b$ is the y-intercept in log-log space. In linear space, this translates to a power law $R_{p}=e^{b}\cdot P^{m}$.

Actually observing this valley can be challenging, as uncertainties in stellar parameters can obfuscate its presence and render characterization of the valley, to aid in understanding planet formation, quite difficult. Fulton et al. (2017) first observed the gap in distribution of planetary radii using precise radius measurements from the California-Kepler Survey and a sample of 2025 exoplanets. They note the gap as a dearth of planets with radii between 1.5 and 2.0 $R_{E}$ and suggest that this limit essentially separates super-Earths from mini-Neptunes, meaning that it determines whether a planet has a substantial H/He envelope. Using Gaia parallaxes, Fulton & Petigura (2018) released an updated catalog of $\sim$1000 planets with improved precision in the planet radius measurements ($\sim$5%). They found that the mass of the stellar host and the planet’s orbital distance are both key factors that sculpt the distribution of planetary radii, concluding that the gap does indeed have planets in it and is more a dearth of planets around $R_{p}\sim 2R_{E}$. Berger et al. (2018) also released an updated catalog with stellar radii derived from a combination of Gaia DR2 parallaxes and the DR25 Kepler Stellar Properties Catalog. They confirm a gap in the distribution of radii with their updated values, noting that the gap is mostly limited to large incident fluxes ($>200F_{E}$). These three studies did not attempt to constrain the valley as a function of orbital period.

Soon after Fulton et al. (2017) released their study, Van Eylen et al. (2018) explored the location and shape of the radius valley using 117 planets with accurate parameters determined from asteroseismology. Using a linear model (in log-log space) and an MCMC algorithm, they first characterize the valley by fitting the absence of data, recovering a slope of $m=-0.10\pm 0.03$. They then use support vector machines to determine the hyperplane of maximum separation between the planets above and below the valley, resulting in a slope of $m=-0.09^{+0.02}_{0.04}$. They conclude that the negative slope they measured is consistent with models of photoevaporation.

There are a few possible explanations for this radius valley aside from mass loss by photoevaporation, including sculpting by giant impacts (e.g., Liu et al., 2015; Schlichting et al., 2015; Inamdar & Schlichting, 2016), core-powered mass loss (e.g., Ginzburg et al., 2017; Gupta & Schlichting, 2018), and planets forming late in a gas-poor disk (e.g., Lee et al., 2014; Lee & Chiang, 2016). It is unclear whether impacts alone can cause a gap in the distribution of planetary radii since impact erosion is a highly stochastic process, but atmospheric heating caused by an impact can cause the envelope to expand, which would make it more susceptible to photoevaporation. Core-powered mass loss results in a negative slope in the valley and can sufficiently explain the gap in the distribution of planetary radii (Ginzburg et al., 2018; Gupta & Schlichting, 2018). Lopez & Rice (2016) explore the idea that super-Earths near their stars are actually a separate population of planets that developed late in a gas-poor disk and never accreted significant envelopes, instead of accreting atmospheres that are then later stripped away. They found that the transition radius between this population and other planets depends on orbital period $R\propto P^{m}$, where $m$ is between 0.07 and 0.10. Therefore, forming in a gas-poor disk would result in the slope of a radius valley that is positive.

To date, there has been no observational study that looks to characterize the radius valley, in both radius and orbital period, without suffering from a small sample size. Here, we aim to characterize the valley by measuring its slope and to identify the dominant mechanism that is responsible for sculpting the radii of exoplanets. We investigate the existence and slope of the radius valley using all known exoplanets with orbital period and planetary radius measurements of $P<50,~{}R_{p}<4.0R_{E}$.

The outline of this paper is as follows: in Section 2, we describe the data used in this study. In Section 3, we describe and explain the theories behind the methods and tests used to examine the radius valley. Following this, we show and discuss the results from our analysis in Section 4 and Section 5 before concluding in Section 6.

2 Data

We use all known exoplanets in our analyses, pulling data from the exoplanetarchive111exoplanetarchive.ipac.caltech.edu, data pulled on April 5, 2019. Although exoplanetarchive keeps track of all published parameters for planets, we use only the most recent parameters for the study presented in this paper. We include planets from all detection methods, although Kepler contributed the most data, resulting in 3933 known exoplanets. We only use planets with published values and uncertainties for radius and orbital period and limit our sample to $P<50$ days and $R_{p}<4.0R_{E}$, resulting in 2066 planets. The Solar System was not included in this study.

We perform our analysis on other data sets that have been used to study this phenomenon in the past, particularly those of Fulton et al. (2017), Van Eylen et al. (2018), and Fulton & Petigura (2018). We refer to these data sets as ”F17,” ”V18,” and ”F18” from here on. We present information about each data set, including number of planets and average uncertainty in planetary radius, in Table 1.

3 Methods and Theory

We analyze the exoplanetary population using 2-D kernel density estimator, hierarchical clustering, $K$-means clustering, $K$-medians clustering, and density-based clustering. The methods used here and the theories behind them are briefly described below. We perform all of our analysis in the R language (R Core Team, 2017).

3.1 Kernel Density Estimator

A kernel density estimation is a nonparametric estimation of probability density functions. The kernel density estimate is defined by:

[TABLE]

where $x=(x_{1},x_{2})^{T}$ and $X_{i}=(X_{i1},X_{i2})^{T}$ for $i$=1,2,…,n a bivariate random sample drawn from a density $f$, $K(x)$ is the kernel which is a symmetric probability density function, and H is the bandwidth matrix which is symmetric and positive-definite. Although the choice of the kernel is not very important, in this analysis we use an axis-aligned bivariate normal kernel, evaluated on a square grid. A KDE is dependent on the distribution of data given to it, so our results will be sensitive to incompleteness and effects from biases. We also must choose a bandwidth for both orbital period and planetary radius; the results of the KDE are not very dependent on the bandwidths, but bandwidths that are too large will result in a uni-modal distribution and bandwidths that are too small will result in multi-modal distributions. We choose bandwidths that lead to bi-modal distributions so that we are able to characterize a line between the modes. In our analysis, we use the function kde2d from CRAN package MASS (Venables & Ripley, 2002). For a more in-depth description of KDEs and their applications, we encourage the reader to explore Silverman (2018).

3.2 Hierarchical Clustering

Hierarchical clustering is a method of cluster analysis that builds a hierarchy of clusters. Originally, each object receives its own cluster; the algorithm computes the distance between clusters via the Lance-Williams dissimilarity formula and combines the two most similar clusters. The Lance-Williams algorithms are a family of agglomerative clustering algorithms which are represented by a recursive formula. After two clusters are merged, the updated cluster distance is computed recursively:

[TABLE]

where $i$ and $j$ indicate the clusters being merged, $d$ is the cluster distance and $\alpha_{i}$, $\alpha_{j}$, $\beta$, and $\gamma$ are parameters. We use complete linkage method to find similar clusters, where the distance between two clusters is the distance between the two elements (one in each cluster) that are farthest away from each other. This process continues until all objects have been sorted into one cluster, and a decision tree has been formed. Simply put:

1. Assign every point to its own cluster.

2. Compute the distance between each cluster.

3. Merge the two closest clusters.

4. Update the distances between the new cluster and

the original clusters.

5. Repeat steps 3 and 4 until only a single cluster

remains.

We cut the hierarchical tree at two branches so that we have two populations (super-Earths and mini-Neptunes). One caveat to this method is we cannot control along which access the resulting two clusters are separated. In our analysis, we use the function hclust from CRAN package MASS (Venables & Ripley, 2002). We refer the reader to Johnson (1967) for a discussion on various hierarchical clustering schemes.

3.3 K-means Clustering

Originally from signal processing, $K$-means clustering is a method of vector quantization. We use the original algorithm by Hartigan & Wong (1979) which aims to separate the data into $K$ number of groups such that the sum of squares from all points to the cluster centers is minimized. Mathematically, the algorithm aims to find:

[TABLE]

where $\bar{x}=\{x_{1},x_{2},...,x_{n}\}$ is a data set, $\bar{S}=\{S_{1},S_{2},...,S_{k}\}$ is a set of clusters, and $\mu_{i}$ is the mean of points in $S_{i}$. Simply put:

1. Select $K$ points as the initial centeroids.

2. Assign all points to the closest centeroid (mean point

of cluster).

3. Recompute the centroid (mean) of each cluster.

4. Repeat steps 2 and 3 until the centroids no longer

change.

We choose $K=2$ for two clusters. In our analysis, we use the function kmeans from base R (R Core Team, 2017). We refer the reader to MacQueen et al. (1967) or Hartigan & Wong (1979) for a more thorough description of the algorithm.

3.4 K-medians Clustering

$K$-medians clustering is another method of vector quantization. We use the same algorithm as presented above in Section 3.3. The only difference between the two methods is $\mu_{i}$ is now the median of the points in each cluster $S_{i}$, so the centeroid is the median point of the cluster. We choose $K=2$ for two clusters. In our analysis, we use the function kGmedian from CRAN package Gmedian (Cardot, 2017). We encourage the reader to explore Mulvey & Crowder (1979) for a detailed description and analysis of this method.

3.5 Density-based Clustering

Density Based Spatial Clustering of Applications with Noise (DBSCAN) is a clustering technique that was designed to discover the clusters and the noise in a spatial database (Ester et al., 1996). The aim of the algorithm is for each point of a cluster in the ”neighborhood” of a given radius $\epsilon$ to contain at least a minimum number of points. In other words, the density in the ”neighborhood” has to exceed some threshold. The algorithm works as:

1. Start with an arbitrary point $p$ in the data and

retrieve all points that are ”density-reachable,”

i.e. the other point belongs in the ”neighborhood”

of the initially selected point $p$ and the

”neighborhood” contains the minimum number of

points.

2. Repeat step 1 until some points are density-reachable.

3. Calculate the distances between the resulting

clusters.

4. Merge clusters that are at least as dense as the

thinnest cluster if the two clusters are closer than

the given radius $\epsilon$.

5. Repeat steps 3 and 4 until no more mergers occur.

The choice of distance function that is used will determine the shape of the ”neighborhood,” but the resulting clusters can have any shape (Ester et al., 1996). We use the Manhattan distance in 2D space which results in a rectangular neighborhood. We randomly draw our $\epsilon$ and minimum number of points from uniform distributions such that the clustering results in two populations. DBSCAN specifically will account for noise in the data, where noise is simply the data points that do not belong to any clusters. This method is therefore less sensitive to outliers than other methods (Ester et al., 1996; Campello et al., 2013). It also means that the clusters are smaller (spatially, and contain fewer points) than those from the other methods which force all points into a cluster. In our analysis, we use the function dbscan from CRAN packages dbscan (Hahsler & Piekenbrock, 2018). We refer the reader to Ester et al. (1996) for the full mathematical definition and implementation of the algorithm.

3.6 A Comparison and Caveats of Clustering Techniques

We use a variety techniques so that our resulting slope is not dependent on the technique. However, many computer scientists and statisticians have been studying and comparing these methods for decades. The general consensus is that a variety of clustering techniques should be studied for each new data set, given that different methods are sensitive to different pit-falls (e.g., Ester et al., 1996; Chen et al., 2002; Mingoti & Lima, 2006; Sharma et al., 2012).

In general, hierarchical clustering or density-based clustering are considered to be superior when the number of clusters is unconstrained or poorly constrained. Partitioning methods such as $K$-means or $K$-medians, however, can be much faster and less computationally intensive (e.g., Chen et al., 2002; Sharma et al., 2012). For the problem presented here, we know that we want two populations (super-Earths and mini-Neptunes) so this superiority no longer applies. Instead, the hierarchical clustering introduces a bit of a problem in our slope determination: given that all of the data have been sorted into one cluster via a tree, and we simply cut this tree at two branches, the axis of separation between the two clusters can readily be along orbital period instead of planetary radius. Oftentimes, the two populations will be separated along $R_{p}\sim 2R_{E}$, but sometimes the populations will be separated by $P\sim 7$ days. This results in lines that are near-vertical with very negative slopes.

Chen et al. (2002) found that, in analyzing ES cell gene expression data, hierarchical clustering performed worse than $K$-means clustering. They attributed its poor performance to the algorithm’s ”greediness” since it is not possible to do any refinement or correction once two clusters have merged. Other studies have found similar results (e.g., Mingoti & Lima, 2006).

DBSCAN has been noted to be inefficient for high-dimensional data sets, but is often praised for its notion of noise and ability to find arbitrarily shaped clusters, including clusters that are completely surrounded by other clusters (Ester et al., 1996; Sharma et al., 2012; Campello et al., 2013). It has also been found to be at least 100 times more efficient than CLARANS, a $K$-medoid algorithm similar to $K$-means and $K$-medians (Ester et al., 1996).

Both $K$-means and $K$-medians clustering are typically faster than hierarchical clustering, especially for large data sets (Ester et al., 1996; Sharma et al., 2012), but do not work very well with non-globular clusters (Sharma et al., 2012). $K$-means clustering is also affected by the presence of a large amount of outliers $\sim$40% (Mingoti & Lima, 2006). The accuracy of $K$-means and $K$-medians clustering is also quite dependent on the choice of initial centeroids (Milligan, 1980). To help mitigate this dependence, our bootstrapping technique starts with different centeroids each time for both methods. Our clusters are also fairly globular in shape and do not contain so many outliers, so these pit-falls do not effect our results.

4 Results

Here we analyze the resulting slopes from our various methods and data sets. We report our bestfit slope and its uncertainty as the median and 16th and 84th percentiles from a combination of all of our bootstrap realizations from using all known exoplanets with period $P<50$ days and radius $R_{p}<4.0R_{E}$. We report the resulting slopes from our different data sets and methods in Table 2.

4.1 Recovering a Line from Between the Two Populations

We aim to recover the slope of the linear separation between the two populations to help determine the cause of this observed dichotomy. For the 2-D KDE, we perform a form of bootstrapping to determine the median and uncertainty in the slope. For each of our 200 bootstrap samples, we draw (with replacement) 2000 planets from the data set. We randomly select our bandwidths for orbital period and for planetary radius from U[2.0,15.0] days and U[0.2,0.4] $R_{E}$, respectively, for each realization. These bandwidths mark the limits for a bi-modal distribution. For each realization, we first find the line between the two peaks in the full period-radius distribution, where period is in days and planetary radius is in $R_{\oplus}$. We then take the line that is perpendicular to this line to be our separator.

For each of our clustering methods (hierarchical, $K$-means, $K$-medians, density-based), we fit a line separating the two categories using a linear discriminant analysis, a generalization of Fisher’s linear discriminant (Fisher, 1936) that aims to find a linear combination of features that separates two populations of objects. Although an in-depth derivation of this process is beyond the scope of this work, we refer the reader to Mika et al. (1999) for an excellent review. We again perform a form of bootstrapping where we draw (with replacement) 2000 planets from the data set 200 times. For the density-based clustering, we choose a bandwidth $\epsilon$ from U[0.315,0.362] for each realization which results in two clusters.

4.2 Resulting slopes and intercepts

We apply a two-dimensional kernel density estimator (2-D KDE) to the data with an axis-aligned bivariate normal kernel, evaluated on a square grid. We perform a bootstrapping variant, sampling our data and measuring the line between populations 200 times. For each sample, we draw with replacement 2000 planets from the data set (see Section 2 for a description of each set). We plot the KDE with the resulting separators in Figure 1. The resulting slope estimation for ”all” data is $m=-0.35^{+0.15}_{-0.14}$ where our uncertainties mark the 16th and 84th percentiles. We present the slope estimates and uncertainties for all methods and data sets in Table 2.

We next cluster our data using a variety of methods described in Section 3: hierarchical clustering, $K$-means clustering, $K$-medians clustering, and density-based clustering. For each of our clustering methods, we repeat the bootstraping variant as mentioned above, selecting 2000 planets with replacement from the data set. We then cluster the data into two clusters for the $K$-means and $K$-medians clustering. For the hierarchical clustering, we cut the tree at two branches since hierarchical clustering clusters the data into an entire tree of classification. This sometimes led to the clustering being solely period-based, resulting in very large and negative slopes. For each clustering method and each bootstrapped realization, we fit a line between the two clusters, as described above in Section 4.1, and take the median slope to be our bestfit and the 16th and 84th percentiles as our uncertainties in the slope. The resulting slope estimation between the two populations for ”all” data is $m=-0.316^{+0.077}_{-0.084}$ for all clustering methods with a resulting y-intercept of $b=1.26^{+0.28}_{-0.17}$. We show a realization of clustering with all fitted lines plotted on top in Figure 2.

We present the slope estimates and uncertainties for all methods and data sets in Table 2. We take our bestfit slope to be the result of ”all” data set and all methods $m=-0.319^{+0.088}_{-0.116}$. All methods and data sets result in consistent and distinctly negative slopes, which are consistent with both models of photoevaporation (e.g., Owen & Wu, 2017) and core-powered mass loss (e.g., Gupta & Schlichting, 2018), but inconsistent with planets forming late in a gas-poor disk.

5 Discussion

We examine the distributions of planetary radii for each of the data sets in an attempt to observe the gap noted by previous studies (e.g., Fulton et al., 2017; Van Eylen et al., 2018). First, we create samples of equal size from each data set by drawing with replacement 1000 times. We also incorporate uncertainties in the planetary radii by drawing each planet radius from a normal distribution that is centered on its nominal value with a standard deviation of the published uncertainty. Our sample of all exoplanets with $P<50$ days and $R_{P}<4.0$ has a much larger average uncertainty for planetary radii than the other samples (see Table 1 for a full comparison of our data sets). In Figure 3, we plot a Kernel Density Estimation (KDE) for the radii distributions for each of our re-sampled data sets. We use a KDE to minimize dependence on bin width, although a KDE does require a bandwidth. We note that, once a KDE is used and uncertainties are accounted for, the gap in distribution of planetary radii is barely distinguishable for data set F18, and not at all distinguishable for the F17 and ”all” data sets. The gap is still prominent in data set V18, which has few planets (81) and the smallest radii uncertainties ($<\sigma_{R}>=0.06R_{E}$).

Accounting for the uncertainty in planetary radii obscures not only the gap in the distribution of radii, but also the radius valley as a function of period. When we cluster the ”all” data set, the slope of the line between the two populations grows from $m=-0.319^{+0.088}_{-0.116}$ using only nominal values for the radii to $m=-0.531^{+0.265}_{-0.401}$ when we account for uncertainties. This large increase in slope is due to the clustering algorithms failing. Instead of clustering along $R_{p}\sim 2R_{E}$, they tend to cluster along $P\sim 7$ days, which results in near vertical lines. This effect becomes even more apparent when we only use planets with $P<25$ days, as can be seen in Table 2.

We repeat our analysis with a cut in orbital period of $P<25$ days to help mitigate the effects of incompleteness due to observational biases. The resulting slope measurements are consistent with the measurements from the full sample, although larger and less precise (see Table 2 for all slope measurements).

The different methods that we employ result in different slopes. In Figure 4, we plot the KDE of the resulting slopes from our analysis on all exoplanets using the nominal values for radius (i.e., not accounting for uncertainties). Slopes larger than $m\sim-0.5$ are a result of the clustering methods failing, which is most apparent in the hierarchical clustering. This method creates an entire tree of classification and we only cut it at two branches. For the lines with steep slopes, a cut at three branches did result in a cluster with $R_{p}<2.5R_{E}$ and a cluster with $R_{p}>2.5R_{E}$, but we only include the results of the two-branch cut. $K$-means and $K$-medians clustering result in much higher precision than the other methods (for the ”all” data set as well as most of the other data sets). We take the median and uncertainties from all slopes from all methods as our bestfit slope $m=-0.319^{+0.088}_{-0.116}$.

Our measured slopes are much steeper than those reported in Van Eylen et al. (2018) and from theoretical models of photoevaporation and core-powered mass loss, although they are still consistent within $3\sigma$. Van Eylen et al. (2018) fit a line to their absence of data and then used support vector machines to determine the hyperplane of maximum separation between the planets above and below the valley, resulting in slope of $m=-0.10\pm 0.03$ and $m=-0.09^{+0.02}_{0.04}$, respectively. Given that the other data sets have planets occupying the radius valley, we cannot employ their methods on other data sets. Our methods (discussed in full in Section 3) result in a steeper slope when even we use the smaller asteroseismic data set from Van Eylen et al. (2018). This difference is due to the different classification techniques; as can be seen in comparing our classification (Figure 2) with the classification from support vector machines (Fig. 7 in Van Eylen et al., 2018), the methods employed herein tend to classify planets with periods $P>10$ days as part of the mini-Neptune population.

There is a dearth of mini-Neptunes at short orbital periods ($P<3$ days), that is predicted by photoevaporation models (e.g., Owen & Wu, 2013) and noted by many authors (e.g., Lundkvist et al., 2018). This absence of detected planets makes a slope difficult to constrain. Coupled with the observational bias of super-Earths far from their hosts being difficult to detect, this dearth of planets can confuse the clustering and lead to positive slopes. How much this effects the resulting slope is difficult to truly explore, but we do note that the clustering algorithms rarely ($<0.1\%$) result in a positive slope. As mentioned above, we partially investigate this by repeating our analysis with an orbital period cut of $P<25$ days. With this cut, we find slopes that are consistent with the slopes from using all data.

Given that the measured slopes for all data sets and methods are consistent and distinctly negative, we find that the radius valley is consistent with both models of photoevaporation (e.g., Owen & Wu, 2017) and core-powered mass loss (e.g., Gupta & Schlichting, 2018), but inconsistent with planets forming late in a gas-poor disk (e.g., Lee & Chiang, 2016). Core-powered mass loss, however, would produce a radius valley that is independent of incident flux, and is therefore independent of stellar type. Recently, Fulton & Petigura (2018) found that the radius gap shifts to higher incident stellar fluxes around higher mass stars, showing that the gap is dependent on stellar type. This is consistent with predictions from photoevaporation models (Owen & Wu, 2017). We note that these mechanisms are not mutually exclusive. It is plausible that both photoevaporation and core-powered mass loss, as well as other mechanisms, are responsible for the radius valley simultaneously or at different times during planetary formation.

The precise slope of the radius valley depends on the planet formation model and the composition of the planets. We compare our resulting slope with different models of photoevaporation from Owen & Wu (2017) in Figure 5. The photoevaporation models use the maximum radius at the bottom of the valley, but since the observational valley is more of a dearth of planets than a complete absence of planets, we do not have a ”bottom” of valley to directly compare to. Instead, we plot our resulting slope from all methods for our ”all” data set with a $1\sigma$ confidence interval, which is an overestimation for what we are comparing to.

Previous studies have found that the location of the radius valley is more consistent with iron-rich cores than with icy cores, assuming that the valley is caused primarily by photoevaporation (Owen & Wu, 2017; Jin & Mordasini, 2018; Van Eylen et al., 2018). We find that our measurements are consistent at $1\sigma$ with both the more complex models that include recombination and x-ray evaporation and with the steeper slopes predicted for pure energy-limited evaporation for both iron-rich and icy cores up until an orbital period of $P\sim 2.5$ days (although our measurements are consistent with all four models from Owen & Wu (2017) at $3\sigma$). Our slopes are more consistent with icy cores than iron-rich cores with orbital periods shorter than $P\sim 2.5$ days, possibly indicating a need for including planetary migration after the planets form.

6 Conclusion

We investigate the relationship between planetary radius and orbital period, confirming the “radius valley" that is suggested to separate super-Earths from mini-Neptunes using 2-dimensional kernel density estimator and various clustering techniques. Unlike previous studies, we use a large sample of exoplanets. We use all known exoplanets, but limit our sample to planets with orbital period and planetary radius measurements and $P<50$ days and $R_{P}<4.0R_{E}$, resulting in 2066 planets. We repeat our analysis for a sample with planetary radii drawn from normal distributions to account for the uncertainty in the measurements, which can obfuscate the radius valley, and for a sample with $P<25$ days to help mitigate the effects of incompleteness due to observational biases. Accounting for uncertainties and a smaller period cut both led to estimates of the slope that were larger and less precise than using the full sample with nominal values for the radii, but still consistent with each other and our other estimates. We also repeat our analysis on the data sets from Fulton et al. (2017), Van Eylen et al. (2018), and Fulton & Petigura (2018), measuring slopes that are consistent with our other analyses. We present all of our measurements of the radius valley for each method and data set in Table 2.

We conclude the following results:

• •

Powerful machine-learning techniques allow us to detect and characterize the radius valley using a large sample of exoplanets.

• •

The transitional radius between super-Earth and mini-Neptune decreases as a function of orbital period. We find a negative slope which is consistent with models of photoevaporation and core-powered mass loss, but is inconsistent with late formation in a gas-poor disk (characterized by a positive slope).

• •

We measure the location of the radius valley as a power law: $\textrm{log}~{}R_{p}=m\cdot\textrm{log}~{}P+b$, where $m=-0.319^{+0.088}_{-0.116}$ and $b=1.26^{+0.28}_{-0.17}$.

• •

Given the lack of a clear valley, our sample consists of planets with a wide range of core compositions or planets which have formed beyond the snow line. Our measured slope is consistent at $1\sigma$ with theoretical models of photoevaporation for cores consisting of a significant fraction of iron and for cores consisting of 1/3 ice, 2/3 silicates for planets with $P>2.5$ days. For smaller orbital periods, our measurements are more consistent with icy cores.

Given the exoplanetary missions that have recently launched or will soon launch and the missions that have been proposed for the future, such an analysis as presented herein should be repeated after more data have been gathered. The addition of data, especially data that are less biased or have different biases, could grandly expand the significance of this study, providing additional information about the formation and evolution of exoplanetary systems.

Acknowledgments

We thank the referee for helpful comments and suggestions that have improved this manuscript. We thank Kathryn Disher for her helpful discussion. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1255832. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.