Predicting Exoplanet Masses and Radii: A Nonparametric Approach

There have been several M–R relations proposed in the exoplanet literature; all of them assume that exoplanet masses and radii follow a power law:

• 1.  
  Lissauer et al. (2011) fit a relation to Earth and Saturn and obtained M = R^2.06, where M and R are in Earth units.
• 2.  
  Wu & Lithwick (2013) calculated the masses of 22 planet pairs based on the amplitudes of their sinusoidal transit timing variations (TTVs) and found that M = 3R.
• 3.  
  Weiss & Marcy (2014) used the 42 planets chosen by the Kepler team to be followed up with radial velocity measurements (Marcy et al. 2014) and found the power law is M = 2.69R^0.93 for planets with 1.5 < R < 4 R_⊕.
• 4.  
  Hadden & Lithwick (2014) revisited the M–R relation for 56 low-eccentricity TTV planets and fit M = 14.9 (R/3R_⊕)^0.65 to them.
• 5.  
  Wolfgang et al. (2016; hereafter denoted as WRF16) selected planets with radii <8 R_⊕ from the NASA Exoplanet Archive (up to the date 2015 January 30) and found two different power-law relations: M = 2.7R^1.3 for <4 R_⊕ planets and M = 1.6R^1.8 for <8 R_⊕ planets. Unlike the earlier results, which used basic least squares regression, WRF16 built a hierarchical Bayesian model that incorporated observational measurement uncertainty, the intrinsic, astrophysical scatter of planet masses, and non-detections and upper limits.
• 6.  
  Mills & Mazeh (2017) fit power laws to various subsets of exoplanet masses and radii, testing for differences between RV and TTV M–R relations in two different period bins.
• 7.  
  Chen & Kipping (2017) built on the hierarchical Bayesian model from WRF16 to fit a continuous broken power-law model across a much larger range of masses and radii. They found that a four-segment M–R relation best described their data, with Terran, Neptunian, Jovian, and Stellar regions. They fit a different power law and intrinsic scatter to each segment.

The power-law model is simple to fit to exoplanet masses and radii, and its parameters are easy to interpret. However, there are several reasons to move beyond it. First, we already know that a single power-law relation will not be able to describe the M–R relation for the entire population, as degeneracy pressure causes the radius–mass relation to flatten out around a Jupiter mass. Chen & Kipping (2017) addressed this by assuming a broken power law with multiple segments. As soon as one expands the space of models beyond those with a handful of parameters, it is useful to consider more flexible nonparametric models as well.⁶ Indeed, there is no particular reason to believe that exoplanet masses and radii should follow power laws as a population, because the processes that dominate planet formation for small rocky planets are different from those that dominate the formation of gas giants. Furthermore, we do not expect the distribution of masses at a given radius to be Gaussian, as was assumed by WRF16 and Chen & Kipping (2017), or even symmetric. There is already evidence from the hierarchical modeling checking performed by WRF16 that a Gaussian distribution does not completely reproduce the observed exoplanet masses and radii.

Moving forward, we adopt a nonparametric approach that allows us to relax these assumptions. There are many benefits of using nonparametric models, including

• 1.  
  They can take on variety of shapes to fit the data, which can be advantageous for making predictions that are less model-dependent. Therefore, there is no need to impose abrupt breaks for modeling M–R relations across a wide range of sizes.
• 2.  
  They can model the joint distribution of mass and radius, and thus provide a self-consistent way to predict both mass from radius and radius from mass over the entire exoplanet mass and radius range.
• 3.  
  Eventually we want to understand masses and radii as a function of many other star and planet properties like stellar mass, metallicity, orbital period, planet multiplicity, and even eccentricity. We do not have clear guidance from theory about the functional form of the dependence in these additional dimensions, and so we should not expect the mass–radius-X distribution to follow a power law.

Our nonparametric model for the joint probability distribution of mass and radius is a bivariate distribution that consists of a tensor product of sequence of location-scaled and transformed beta densities that can also be viewed as Bernstein polynomials. Each marginal distribution can be expressed as a linear combination of Bernstein polynomials (BPs), leading to a mixture of beta densities. When normalized, BPs have the same functional form as beta distributions (see Equation (3)). For a visual representation of Bernstein polynomials, see Figure 9, where we plot BPs as a function of their degree (denoted by k or l in Equation (7)).

The properties of BPs and a discussion of their advantages over other choices for basis functions are described in Appendix A. In addition, we note that we do not use Gaussian process regression, which is a popular nonparametric prediction method in exoplanet science, for a number of reasons. First, a joint density estimation approach is more appropriate for our purposes than a regression approach, as argued previously. Second, the marginal distribution of the M–R is not in general normally distributed, as we demonstrate in Section 5. Third, fitting a mixture of Gaussians to describe the joint density would involve three free parameters per component (mixture weight, mean, and variance), rather than the one free parameter per component that we use here (the mixture weight). Fourth, a mixture of Gaussians produces an identifiability problem, wherein different components can be interchangeable, that is not present in the Bernstein polynomial model due to the ordered nature of the different terms (see Equation (7) and the figures in Appendix A for how the Bernstein polynomial terms; i.e., the different β_k for a given d are distinguishable from each other).

Before introducing our model, we first define some mathematical notations that will be used throughout the paper. Let ${M}_{i}^{\mathrm{obs}}$ and ${R}_{i}^{\mathrm{obs}}$ be the reported mass and radius measurements for ith planet in our data set, and ${\sigma }_{{M}_{i}}^{\mathrm{obs}}$ and ${\sigma }_{{R}_{i}}^{\mathrm{obs}}$ be the reported standard deviations for their measurement errors. We denote M_i and R_i as the true mass and radius for ith planet that would have been observed if there were no measurement errors. We denote their joint distribution as f(m, r), and the M–R relation as the conditional distribution $f(m| r)$. We let ${{\boldsymbol{M}}}^{\mathrm{obs}}$ stand for a vector containing the observations $({M}_{1}^{\mathrm{obs}},\,...,\,{M}_{n}^{\mathrm{obs}});$ similarly, ${{\boldsymbol{R}}}^{\mathrm{obs}}$, ${{\boldsymbol{\sigma }}}_{M}^{\mathrm{obs}}$, and ${{\boldsymbol{\sigma }}}_{R}^{\mathrm{obs}}$ stand for the set of their respective observations. We let ${ \mathcal N }(\mu ,\sigma )$ stand for a normal distribution with mean μ and standard deviation σ, and denote the probability density of a Beta distribution with shape parameters i and $d-i+1$ by

$$\begin{eqnarray}&&{\beta }_{\mathrm{id}}(a)=d\displaystyle \left(\genfrac{}{}{0em}{}{d-1}{i-1}\right){a}^{i-1}{(1-a)}^{d-i},\end{eqnarray} \tag{ 3 }$$

where 0 ≤ a ≤ 1, and d > 0 is an integer that relates to the degree of a Bernstein polynomial.

Our hierarchical model can be written as follows:

$$\begin{eqnarray}&&{M}_{i}^{\mathrm{obs}}\mathop{\sim }\limits^{\mathrm{ind}}{ \mathcal N }({M}_{i},{\sigma }_{{M}_{i}}^{\mathrm{obs}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{R}_{i}^{\mathrm{obs}}\mathop{\sim }\limits^{\mathrm{ind}}{ \mathcal N }({R}_{i},{\sigma }_{{R}_{i}}^{\mathrm{obs}}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&({M}_{i},{R}_{i})\mathop{\sim }\limits^{\mathrm{iid}}f(m,r| {\boldsymbol{w}},d,d^{\prime} ),\end{eqnarray} \tag{ 6 }$$

where the symbol $\mathop{\sim }\limits^{\mathrm{ind}}$ stands for "independently distributed as," and the symbol $\mathop{\sim }\limits^{\mathrm{iid}}$ stands for "identically and independently distributed as." Let ${\boldsymbol{w}}=({w}_{11},\,...,\,{w}_{{dd}^{\prime} }$) be a set of weights that describes how much each corresponding term in the following series contributes to the overall joint distribution f(m, r):

$$\begin{eqnarray}&&f(m,r| {\boldsymbol{w}},d,d^{\prime} )=\displaystyle \sum _{k=1}^{d}\displaystyle \sum _{l=1}^{d^{\prime} }{w}_{\mathrm{kl}}\displaystyle \frac{{\beta }_{{kd}}\left(\tfrac{m-\underline{M}}{\overline{M}-\underline{M}}\right)}{\overline{M}-\underline{M}}\displaystyle \frac{{\beta }_{{ld}^{\prime} }\left(\tfrac{r-\underline{R}}{\overline{R}-\underline{R}}\right)}{\overline{R}-\underline{R}}.\end{eqnarray} \tag{ 7 }$$

To ensure that Equation (7) remains a probability distribution (i.e., that it integrates to 1), we impose the constraint ${\sum }_{k=1}^{d}{\sum }_{l=1}^{d^{\prime} }{w}_{{kl}}=1$ and that w_kl ≥ 0 for all values of k and l. The notations $\underline{M}$, $\overline{M}$, $\underline{R}$, and $\overline{R}$ are used to denote the lower and upper bounds, respectively, for mass and radius. The values of the upper and lower bounds can be determined by a set of observed values of masses and radii; for example, one could choose the lower mass bound to be the minimum (${M}_{i}^{\mathrm{obs}}-{\sigma }_{{M}_{i}}^{\mathrm{obs}}/n$) in the data set and the upper mass bound to be the maximum (${M}_{i}^{\mathrm{obs}}+{\sigma }_{{M}_{i}}^{\mathrm{obs}}/n$). The bounds could also be set as the minimum and maximum mass and radius theoretically expected for an exoplanet. For regions with no observations, the BP model reverts to the overall mean of the whole function; this happens when the values of $\underline{M}$, $\overline{M}$, $\underline{R}$, and $\overline{R}$ are chosen to be far away from the nearest observations (see the discussion in Appendix B.1).

Note that Equations (4)–(6) form a multi-level model, in that the measurement process is modeled as a separate random process from the underlying true distribution of exoplanet masses and radii. It therefore has a similar hierarchical structure as the hierarchical Bayesian models defined in WRF16 and Chen & Kipping (2017).

Although the model may look complex at first glance, the inference is relatively straightforward. There are only two types of unknown parameters in the model: the weights ${\boldsymbol{w}}$ and the degrees d and d'. Once d and d' are determined using the method that we will introduce in the later part of this section. The two parameters i and $d-i+1$ in each beta distribution (see Equation (3)) are also determined. This is advantageous compared to other mixture models; for example, a mixture of Gaussians would require that each mean and standard deviation of component distribution be estimated by the data. Moreover, it avoids the common identifiability issues created by label switching within the mixture of Gaussians. Another advantage of our mixture of beta distributions is that parameter space (which is a simplex) is a bounded convex set, which helps to guarantee the existence of the optimally estimated parameters. Here we take a maximum likelihood approach to estimate the parameters by explicitly deriving the likelihood using one-dimensional numerical integration.

The likelihood function of the model (4)–(6) can be written as

$$\begin{eqnarray}&&\begin{array}{l}L({\boldsymbol{w}},d,d^{\prime} | {{\boldsymbol{M}}}^{\mathrm{obs}},{{\boldsymbol{R}}}^{\mathrm{obs}},{{\boldsymbol{\sigma }}}_{M}^{\mathrm{obs}},{{\boldsymbol{\sigma }}}_{R}^{\mathrm{obs}})\\ =\ {\int }_{\underline{M}}^{\overline{M}}{\int }_{\underline{R}}^{\overline{R}}f({{\boldsymbol{M}}}^{\mathrm{obs}},{{\boldsymbol{R}}}^{\mathrm{obs}},m,r| {\boldsymbol{w}},d,d^{\prime} ,{{\boldsymbol{\sigma }}}_{M}^{\mathrm{obs}},{{\boldsymbol{\sigma }}}_{R}^{\mathrm{obs}})\ {drdm}\\ =\ \displaystyle \prod _{i=1}^{n}{\int }_{\underline{M}}^{\overline{M}}{\int }_{\underline{R}}^{\overline{R}}f({M}_{i}^{\mathrm{obs}}| m,{\sigma }_{{M}_{i}}^{\mathrm{obs}})f({R}_{i}^{\mathrm{obs}}| r,{\sigma }_{{R}_{i}}^{\mathrm{obs}})\\ \quad \times \ f(m,r| {\boldsymbol{w}},d,d^{\prime} )\ {drdm}\\ =\ \displaystyle \prod _{i=1}^{n}\displaystyle \sum _{k=1}^{d}\displaystyle \sum _{l=1}^{d^{\prime} }{w}_{{kl}}{\int }_{\underline{M}}^{\overline{M}}\displaystyle \frac{1}{{\sigma }_{{M}_{i}}^{\mathrm{obs}}}{\varphi }\left(\displaystyle \frac{{M}_{i}^{\mathrm{obs}}-m}{{\sigma }_{{M}_{i}}^{\mathrm{obs}}}\right)\displaystyle \frac{{\beta }_{{kd}}\left(\tfrac{m-\underline{M}}{\overline{M}-\underline{M}}\right)}{\overline{M}-\underline{M}}\ {dm}\\ \quad \times \ {\int }_{\underline{R}}^{\overline{R}}\displaystyle \frac{1}{{\sigma }_{{R}_{i}}^{\mathrm{obs}}}{\varphi }\left(\displaystyle \frac{{R}_{i}^{\mathrm{obs}}-r}{{\sigma }_{{R}_{i}}^{\mathrm{obs}}}\right)\displaystyle \frac{{\beta }_{{ld}^{\prime} }(\tfrac{r-\underline{R}}{\overline{R}-\underline{R}})}{\overline{R}-\underline{R}}\ {dr},\end{array}\end{eqnarray} \tag{ 8 }$$

where ${\varphi }$(·) is the standard normal distribution, the two integral terms are essentially two constants and can be calculated numerically by using integrate function in ${\mathsf{R}}$ or by any other Gaussian quadrature method available for one-dimensional numerical integration. Therefore, we denote

$$\begin{eqnarray*}\begin{array}{rcl}{c}_{{kl},i} & \equiv & {\int }_{\underline{M}}^{\overline{M}}\displaystyle \frac{1}{{\sigma }_{{M}_{i}}^{\mathrm{obs}}}{\varphi }\left(\displaystyle \frac{{M}_{i}^{\mathrm{obs}}-m}{{\sigma }_{{M}_{i}}^{\mathrm{obs}}}\right)\displaystyle \frac{{\beta }_{{kd}}\left(\tfrac{m-\underline{M}}{\overline{M}-\underline{M}}\right)}{\overline{M}-\underline{M}}\ {dm}\\ & & \times \ {\int }_{\underline{R}}^{\overline{R}}\displaystyle \frac{1}{{\sigma }_{{R}_{i}}^{\mathrm{obs}}}{\varphi }\left(\displaystyle \frac{{R}_{i}^{\mathrm{obs}}-r}{{\sigma }_{{R}_{i}}^{\mathrm{obs}}}\right)\displaystyle \frac{{\beta }_{{ld}^{\prime} }\left(\tfrac{r-\underline{R}}{\overline{R}-\underline{R}}\right)}{\overline{R}-\underline{R}}\ {dr},\end{array}\end{eqnarray*}$$

and we use the symbol ${{\boldsymbol{c}}}_{i}=({c}_{11,i},\,...,{c}_{{dd}^{\prime} ,i})$ to simplify the expression.

Then, for given values d and d', we find an estimator of ${\boldsymbol{w}}$ that maximizes the log-likelihood,

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{maximize}:\quad \mathrm{log}L=\displaystyle \sum _{i=1}^{n}\mathrm{log}({{\boldsymbol{c}}}_{i}^{T}{\boldsymbol{w}}),\\ \mathrm{subject}\ \mathrm{to}:\quad \displaystyle \sum _{k=1}^{d}\displaystyle \sum _{l=1}^{d^{\prime} }{w}_{{kl}}=1,{w}_{{kl}}\geqslant 0,\\ \qquad \qquad \mathrm{for}\ \mathrm{all}\ k=1,\,...,d,l=1,\,...,d^{\prime} .\end{array}\end{eqnarray} \tag{ 9 }$$

This is a convex optimization problem that can be readily solved using any standard numerical optimization methods. For our applications, we have used the ${\mathsf{R}}$ package Rsolnp to solve this optimization problem.

Now we discuss how we choose d and d', which act as tuning parameters relative to the smoothness of the underlying density function. We first choose a set of candidate values for d and d', such as 2, 3, ..., $n/\mathrm{log}(n)$, where n is the sample size (i.e., number of exoplanets). We use 10-fold cross validation to choose the optimal value of degrees d and d' from those candidate values (see more detailed discussions on choosing the number of folds in the cross-validation method in Chapter 7.10 of Hastie et al. 2001). To conduct the 10-fold cross validation, we separate the data set randomly into 10 disjoint subsets with equal sizes. Then we leave out the sth subset, denoted by D_s ($s=1,\,...,\,10$), and use the remaining nine subsets of data to estimate the parameters. Doing this for each D_s in turn results in 10 estimated sets of weights ${\hat{{\boldsymbol{w}}}}^{(s)}$ based on observations $i\notin {D}_{s}$. We plug in each ${\hat{{\boldsymbol{w}}}}^{(s)}$ along with the corresponding data that was used to estimate each set of weights to obtain an estimated value for the log-likelihood. Mathematically, we are calculating the quantity

$$\begin{eqnarray*}&&\mathrm{log}{L}^{\mathrm{pred}}=\displaystyle \sum _{s=1}^{10}\displaystyle \sum _{i\notin {D}_{s}}\mathrm{log}({{\boldsymbol{c}}}_{i}^{T}{\hat{{\boldsymbol{w}}}}^{(s)})\end{eqnarray*}$$

for different possible values of d and d'. The optimum degrees for the BP model are then the d and d' which give the largest value of $-\mathrm{log}{L}^{\mathrm{pred}}$.

Besides using the cross-validation method, other popular methods such as AIC (Akaike information criterion) and BIC (Bayesian information criterion) are also possible when data samples are abundant. When the sample size is relatively small, both methods will under-select the degrees, as in the example provided in Turnbull & Ghosh (2014). This is the reason why we choose to use the 10-fold cross-validation method.

The model may require d and d' to be large, and thus the number of weights that need to be estimated is also large. In such a situation, one may be concerned whether we have enough data to estimate the weights. Fortunately, most of the estimated weights turn out to be numerically zero, due to the fact that parameter space is simplex and therefore convex. This drastically reduces the number of effective parameters in our model. For example, mass predictions computed using the largest 25 weights of the fit performed in Section 5 are indistinguishable from predictions produced using the full set of weights, to the 3rd or 4th decimal in log(M).

Once we obtained $\hat{d}$, $\hat{d}^{\prime}$, and $\hat{{\boldsymbol{w}}}$, the estimated values of d, d', and ${\boldsymbol{w}}$, the M–R relation can be derived following Equation 1:

$$\begin{eqnarray}\begin{array}{rcl}f(m| r,\hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} ) & = & \displaystyle \frac{f(m,r| \hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} )}{{\int }_{\underline{M}}^{\overline{M}}f(m,r| \hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} ){dm}}\\ & = & \displaystyle \frac{f(m,r| \hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} )}{f(r| \hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} )},\end{array}\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray*}&&f(r| \hat{{\boldsymbol{w}}},\hat{d},\hat{d}^{\prime} )=\displaystyle \sum _{k=1}^{\hat{d}}\displaystyle \sum _{l=1}^{\hat{d}^{\prime} }{\hat{w}}_{{kl}}\displaystyle \frac{{\beta }_{l\hat{d}^{\prime} }\left(\tfrac{r-\underline{R}}{\overline{R}-\underline{R}}\right)}{(\overline{R}-\underline{R})}.\end{eqnarray*}$$

From Equation (10), mean, variance, and prediction intervals can be obtained analytically.

As noted in Section 2.1, we could also estimate our nonparametric model from a Bayesian framework by assigning priors to the parameters ${\boldsymbol{w}}$, d, and d'. For example, we can assign a Dirichlet prior to the weights as

$$\begin{eqnarray*}&&\pi ({w}_{11},{w}_{12}...,\,{w}_{{dd}^{\prime} })\sim \mathrm{Dir}({\alpha }_{11},{\alpha }_{12},\,\ldots ,\,{\alpha }_{{dd}^{\prime} }),\end{eqnarray*}$$

where ${\sum }_{k=1}^{d}{\sum }_{l=1}^{d^{\prime} }{\alpha }_{{kl}}=1$, and assign a Poisson or uniform prior for d and d', respectively. The inference can be carried out using a Markov chain Monte Carlo algorithm (MCMC; i.e., Petrone 1999a, 1999b; Petrone & Wasserman 2002). However, initial investigations into this approach indicated that the MCMC algorithm would take a much longer time to run. As a result, we employ a maximum likelihood method rather than a Bayesian method in this paper.

We choose to restrict our science data set (results presented in Section 5) to Kepler planets in an effort to minimize the influence of survey selection effects. Selection effects manifest in two ways for the mass–radius relation:

• 1.  
  The probability of detecting a planet is a non-constant function of either their mass (for RV detection) or radius (for transit detection), with smaller or less massive planets being more difficult to detect.
• 2.  
  The probability that a known planet has its radius (for RV detections) or its mass (for transit detections) measured by follow-up observations is a complicated, often unpublished function of the planet's discovery mass/radius, the predicted radius/mass, and the host star's brightness, spectral type, activity indicators, sky position, and so on.

These two selection effects impact the inferred mass–radius relation in different ways. The first affects the relative amount of data at large versus small radii (or at high versus low mass). Correcting for this effect becomes important when one's scientific goal is a joint mass, radius probability distribution g(m, r) that reflects the underlying, true distribution of exoplanet masses and radii. For example, if one tried to use the data from the Exoplanet Archive as is with no correction for detection completeness, they would incorrectly conclude that ∼1.1R_Jupiter and ∼1M_Jupiter is the most probable mass and radius for an exoplanet, as more Jupiter-sized planets have had their masses measured than smaller planets. We already know that this potential conclusion to be incorrect: the numerous occurrence rate studies which use Kepler data to correct for variable detection efficiency across the full range of planet radii (e.g., most recently Foreman-Mackey et al. 2014; Fulton et al. 2017; Hsu et al. 2018) have shown that smaller planets are much more prevalent than Jupiter-sized ones. Therefore, if one wishes to understand the probability of an exoplanet existing with a specific mass and radius (i.e., characterize g(m, r)), they must account for the different surveys' detection completeness.

Fortunately, this is not the goal of this work. We are scientifically interested in $g(m| r)$, the conditional distribution which gives us the M–R relation. In going from g(m, r) to $g(m| r)$, we marginalize over r, which divides out the completeness correction in that dimension. With this correction disappearing from the equation for $g(m| r)$, we can safely ignore the effect of the selection bias due to detectability of different radii planet. (see Equation (1)). This leaves only the second concern for our mass predictions: how transiting planets are selected for RV follow-up.

To address this concern, one needs to know how planets were chosen for follow-up. To date, no RV follow-up group has published their selection function in a quantitative way that would allow us to incorporate it into an analysis like this. Furthermore, given that the target list often evolves as RV data are acquired, this selection function is probably intractable to calculate for the current data set. Progress is being made toward this end by designing RV follow-up campaigns that have definable selection functions (i.e., Burt et al. 2018 and Montet 2018), but for the current data set, this remains practically out of reach. No other papers on the M–R relation have incorporated corrections for completeness or selection effects into their analysis. We follow this precedent, acknowledging that is a clear area for future work, and choose to focus instead on the novel development of a nonparametric analysis.

The M–R relation is the conditional distribution of log-scaled mass given log-scaled radius, derived from Equation (1):

$$\begin{eqnarray}&&f(\mathrm{log}m| \mathrm{log}r,{\boldsymbol{w}},d,d^{\prime} )=\displaystyle \frac{f(\mathrm{log}m,\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )}{f(\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )},\end{eqnarray} \tag{ 12 }$$

with $f(\mathrm{log}m,\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )$ as in (11) and

$$\begin{eqnarray*}&&f(\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )=\displaystyle \sum _{k=1}^{d}\displaystyle \sum _{l=1}^{d^{\prime} }{w}_{{kl}}\displaystyle \frac{{\beta }_{{ld}}\left(\tfrac{\mathrm{log}r-\mathrm{log}\underline{R}}{\mathrm{log}\overline{R}-\mathrm{log}\underline{R}}\right)}{\mathrm{log}\overline{R}-\mathrm{log}\underline{R}}.\end{eqnarray*}$$

We use 10-fold cross-validation (see Section 2.2) to select d (corresponding to mass) and d' (corresponding to radius), and find that $\hat{d}=55$ and $\hat{d}^{\prime} =50$. After we obtain $\hat{d}$ and $\hat{d}^{\prime}$, we then plug in their values to estimate ${\boldsymbol{w}}$ using the MLE method. Because the values of $\hat{d}$ and $\hat{d}^{\prime}$ may change under different realizations of the random number generator, we repeat the 10-fold cross-validation method five times to assess its stability, each time using a different random seed. We found that four out of five times, the cross-validation gives the same choices for d and d', with the dissenting realization suggesting a slightly larger number for $\hat{d}$. Thus we decide to use $\hat{d}=55$ and $\hat{d}^{\prime} =50$ to estimate the M–R relation. After we obtain $\hat{d}$, $\hat{d}^{\prime}$, and $\hat{{\boldsymbol{w}}}$, we derive the prediction distribution of the M–R relation and employ the bootstrap method described in Section 4 to obtain the bootstrap confidence intervals for the mean.

We plot the mean along with the prediction intervals and the bootstrap confidence intervals in Figure 3. In the plot, the gray region represents the 16% and 84% prediction intervals and the blue region represents the 16% and 84% bootstrap confidence intervals for the mean. The black line represents the estimated M–R relation. From the plot, we see that the M–R relation below 0.8 R_⊕ is unexpectedly similar to a step function. This is because there is only one isolated data point in that region, which the model tries to connect with the nearby points via smooth polynomials. Fortunately, we see that the bootstrap confidence intervals in this region are very wide, indicating that the M–R relation is not well understood in this region, and that more data are needed.

Zoom In Zoom Out Reset image size

The M–R relations can be roughly split into three parts: 0.8 R_⊕ ≤ r ≲ 5 R_⊕, 5 R_⊕ ≲ r ≲ 11 R_⊕, and ≳11 R_⊕. In the region between 0.8 R_⊕ and 5 R_⊕, the bootstrap confidence intervals narrow as data points become abundant for larger radii, which suggests the M–R relation is more accurate. The almost linear relation suggests that the power-law relation may not be a misleading assumption in this radius regime. Transitioning from the 5 R_⊕ region to the 11 R_⊕ region, the number of observations decreases, and we observe both the prediction intervals and bootstrap confidence intervals becoming slightly wider. In the region >11 R_⊕, the M–R relation is more flat or even decreases as the radius becomes larger. These findings are consistent with well-known features of the M–R relation (see the discussion in Section 6.1).

Figure 4 shows how the intrinsic scatter changes as a function of radius; because the intrinsic scatter is defined as the standard deviation of $f(\mathrm{log}m| \mathrm{log}r)$ (see Equation (12)), the intrinsic scatter is technically in units of $\mathrm{log}m$. We plot this on the left and the original scale on the right for a more intuitive comparison. Starting at the smallest radii, we see that the bootstrap confidence intervals are quite wide. This reflects the fact that there is only one data point below ∼0.8 R_⊕, and so the intrinsic scatter is quite uncertain in this regime. In the region between 0.8 R_⊕ and ≲5 R_⊕, the intrinsic scatter is an increasing function with radius. The increasing pattern becomes more obvious in the right plot. This finding is consistent with the result shown in Figure 2, even though the two results are based on different data sets. With more data points to constrain this region, we see that the size of the intrinsic scatter in Figure 4 is smaller than in Figure 2. In the >5 R_⊕ region, the intrinsic scatter behaves nonlinearly. It decreases first and then increases at radii dominated by inflated Jupiters. As discussed in Section 6.1, the M–R relation is difficult to interpret in the Jupiter regime because degeneracy pressure creates a vertical line in this space that is not well captured by the conditional. See the R–M relation in Figure 6 for a clearer picture of the model's behavior in the Jupiter regime.

Zoom In Zoom Out Reset image size

In Figure 5, we plot five conditional densities for mass given radius at different radius values: 1 R_⊕, 3 R_⊕, 5 R_⊕, 10 R_⊕, and 15 R_⊕. These curves represent the mass predictions for planets at those radii (i.e., the probability of a planet of that size having a certain mass). The shaded regions represent the uncertainty in that distribution and correspond to the 68% bootstrap confidence regions. While these distributions look roughly Gaussian in log scale, when transformed to a linear scale, all the densities are skewed to the right. The two densities within the region >5 R_⊕ are more skewed than the rest. Thus no matter the the size of the planets, the intrinsic scatter is definitely not Gaussian in linear scale, as was assumed by WRF16. The lognormal intrinsic scatter assumed by Chen and Kipping (2017) is more appropriate for mass given radius, although it fails to capture the low-mass tail corresponding to super-puffy planets at 7–10 R_⊕ and <30 M_⊕.

Zoom In Zoom Out Reset image size

As we are modeling the joint distribution of mass and radius (in log scale), the radius–mass relation (R–M relation) can also be derived easily. From Equation (2), the conditional distribution for $\mathrm{log}r$ given $\mathrm{log}m$ is

$$\begin{eqnarray}&&f(\mathrm{log}r| \mathrm{log}m,{\boldsymbol{w}},d,d^{\prime} )=\displaystyle \frac{f(\mathrm{log}m,\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )}{f(\mathrm{log}m| {\boldsymbol{w}},d,d^{\prime} )},\end{eqnarray} \tag{ 13 }$$

with $f(\mathrm{log}m,\mathrm{log}r| {\boldsymbol{w}},d,d^{\prime} )$ in (11) and

$$\begin{eqnarray*}&&f(\mathrm{log}m| {\boldsymbol{w}},d,d^{\prime} )=\displaystyle \sum _{k=1}^{d}\displaystyle \sum _{l=1}^{d^{\prime} }{w}_{{kl}}\displaystyle \frac{{\beta }_{{ld}}\left(\tfrac{\mathrm{log}m-\mathrm{log}\underline{M}}{\mathrm{log}\overline{M}-\mathrm{log}\underline{M}}\right)}{\mathrm{log}\overline{M}-\mathrm{log}\underline{M}}.\end{eqnarray*}$$

After we plug in $\hat{d}$, $\hat{d}^{\prime}$, and $\hat{{\boldsymbol{w}}}$, we plot this relation along with the prediction intervals and bootstrap confidence intervals in Figure 6. Note that the R–M relation is not a simple flip of the M–R relation around the 1:1 line; this is due to the fact that the M–R relation is actually a distribution that is asymmetric around a mean relation that is not a straight line.

Zoom In Zoom Out Reset image size

The R–M relation is highly uncertain for <0.8 M_⊕ due to the single isolated data point in that region; this is reflected in the wide bootstrap confidence intervals. The bootstrap confidence intervals narrow as mass increases, reflecting higher certainty in the mean relation, given the data. This mean R–M relation is an increasing function up to 100 M_⊕; above that, the R–M relation becomes flatter due to degeneracy pressure.

We also plot the 68% bootstrap confidence intervals for each of the conditional distributions for 1, 10, 50, 100, and 500 M_⊕. Almost all the densities are skewed to the right, especially when the displayed log scale is transformed to linear scale. In the 0.8 M_⊕ ≤ m ≤ 10² M_⊕ region, the three densities look quite different, indicating rapid change in the radius distribution from 1 to 50 M_⊕.

In this section, we apply the estimated M–R relation to predict masses of K2 planets, given their radius. The K2 planets data are obtained from the NASA Exoplanet Archive on 2017 June 5. In Figure 8, we first plot the prediction intervals from the M–R relation for all radii, then plot the K2 planets on top of the band. We see that for many of the planets with R < 10 R_⊕, the prediction intervals do not cover their masses. However, for planets with R > 10 R_⊕, the prediction intervals do span the masses for all the planets. This suggests that there is a significant bias for K2 planets. One explanation for this is that the most massive planets reach a high mass detection significance threshold first and thus are published first. Burt et al. (2018) discusses this bias and its implication for population mass–radius analyses in detail.

A crucial test of any method is how well it reproduces features that are well established by multiple studies using different techniques. To this end, we emphasize that our relation is not just the average line, but also the predictive distribution around the line. We also note that our goal was not to reproduce prior mass–radius relations, as one of our motivations for applying a nonparametric method was to assess how well the features identified in more parametric models persist when the strong assumptions in those models are relaxed.

To start with, our nonparametric M–R relation reproduces the well-established result that more massive planets ≲10² M_⊕ are on average larger (see Figure 3). This result is not a surprise: more massive planets are expected to form earlier in the protoplanetary disk lifetime and thus should able to accrete more gas-rich material, which quickly increases the planets' radii. Indeed, all M–R relations in the literature (see Section 1.1) show that mass (or radius) is an increasing function of radius (or mass) in the sub-Jupiter regime. However, above about half the mass of Jupiter, degeneracy pressure begins to noticeably affect the radii of hydrogen-dominated bodies (Zapolsky & Salpeter 1969), and the R–M relation should flatten out. Our nonparametric approach readily reproduces this result (see Figure 6).

The electron degeneracy pressure would manifest as a vertical line in the M–R relation of Figure 3; when taking the marginal of the joint mass–radius distribution to produce the $f(m| r)$ (mass given radius) conditional distribution, vertical features like this are collapsed, with their extent represented as the probability distribution at that radius. Figure 7 displays these probability distributions, and shows that there is some probability of more massive planets at 10 and 15 R_⊕. In fact, the difference between $f(m| r)$ and $f(r| m)$ in the Jupiter regime is an argument for why a two-way M–R relation is necessary: there are some features that are easier to capture one way as opposed to the other.

Zoom In Zoom Out Reset image size

Identifying transitions in the exoplanet mass–radius relation is important for our physical understanding of the types of planets that exist, based on their compositions.

No sharp transitions are visible in Figure 3. Nevertheless, there appears to be at least three segments to the M–R relation: ≲5 R_⊕, ∼5 to ∼11 R_⊕, and ≳11 R_⊕ radii, which roughly corresponds to sub-Neptunes, sub-Saturns, and Jupiters. A comparison of Figure 3, the M–R relation, with Figure 6, the R–M relation, shows that there is uncertainty in the sub-Neptune to sub-Saturn transition point. This is because we relax the assumption that the distribution around the average relation is symmetric and constant within a segment, as is assumed by WRF16 and Chen & Kipping (2017) via both studies' use of a Gaussian distribution to characterize the intrinsic scatter term. With this restriction lifted, the average mass given radius is no longer the inverse of the average radius given mass, as the conditionals shown in these two plots are produced by marginalizing over a joint mass and radius distribution that is asymmetric. When the joint distribution is much wider in one-dimension than it is in the other, as is the case around 4–7 R_⊕, then the conditionals can look quite different. A relatively large spread in radii for planets with 3–10 M_⊕ falls adjacent to a small spread in radii for planets around 30 M_⊕; this is what causes ambiguity in identifying a transition from Neptunes to Saturns.

Upon visual inspection, it is also clear that there are not currently enough mass measurements of extrasolar planets ≤2 R_⊕ to justify a super-Earth transition point. This result is at odds with some previously published results on the M–R relation—for example, Weiss & Marcy (2014), Rogers (2015), Chen & Kipping (2017), and Fulton et al. (2017). In the following discussion, we will address each of these works individually.

Weiss & Marcy (2014) fixed the transition point at 1.5 R_⊕ based on a visual inspection of the density versus radius plane. When their data set is plotted with the reported error bars, it becomes apparent that there is little empirical evidence for a transition at 1.5 R_⊕, given the size of the density or mass error bars for the smallest planets. Indeed, one could draw a straight line through the mass and radius points and have it fit the data just as well as the chosen relation with the transition.

To our best knowledge, Rogers (2015) provided the strongest evidence of a transition in the M, R plane among all four of these studies. This paper parameterized the transition as the radius below which 100% of planets must be rocky (i.e., their masses are greater than the minimum mass of a rocky planet at that size). This study tested both a step function for the fraction of planets that are rocky and a more gradual transition; given the size of the error bars, a step function was sufficient to describe this data set. While useful in its physical interpretation, this parametrization does not guarantee that there is a visible kink or transition in the M–R relation. Indeed, no such kink is visible in their data set, displayed in their Figure 1. This is because the Rogers parameterization of the transition tests for the absence of planets in one region of the M, R parameter space (that of very low-mass planets in the <1.5 R_⊕ radius regime); if their absence does not sufficiently alter the average M–R relation, then it will not produce a kink in the M–R plane. Considering the very large measurement uncertainties in this very small planet regime, it is indeed the case that the absence of these planets does not significantly shift the average M–R relation to produce a kink.

Fulton et al. (2017) looked only at the marginal radius distribution. There is no guarantee that the bimodality they see in the radius distribution will produce a kink, or a transition, in the M–R plane. Indeed, the transition is likely gradual with a period dependence. Additional follow-up to measure masses for a larger sample of these small planets is needed before a claim of a Fulton-like gap or transition in the M–R plane is warranted.

The tension between the existence of a super-Earth transition in Chen & Kipping (2017) and our lack of one arises from the fact that we do not include the minor Solar System bodies in our data set. It is not unreasonable to expect that the compositions of mainly refractory bodies in other planetary systems will be similar to those in ours. However, including this data without inflating the uncertainties to match those we obtain for extrasolar systems leads to misleadingly tight constraints for ∼1 M_⊕ planets. Indeed, there is only one planet between 1 and 3 M_⊕ in their data set, so most of the information about the "Terran" transition point comes from extrapolating up from Solar System minor bodies and extrapolating down from Neptunes, rather from measurements at the transition point itself. We choose a less Solar System–constrained approach to characterize the distribution of compositions in the 1–3 R_⊕ regime, choosing instead to develop a flexible method that can best let the extrasolar data speak for itself as more super-Earths are discovered and followed up by TESS.

The intrinsic dispersion in mass given radius slightly increases from 1 to 5 R_⊕ (see Figure 4), but varies much more in radius given mass (see the width of the gray region in Figure 6). This indicates that planet compositions in the 3–10 M_⊕ range exhibit more diversity than others, and that planets with 1–5 R_⊕ fall into a relatively narrow range of masses. The large spread in radius over a narrow mass regime could have a number of physical interpretations: perhaps these super-Earths/sub-Neptunes fall at a transition point where scaled-up terrestrial planet formation co-exists with scaled-down gaseous planet formation. Perhaps planet evolution is particularly dramatic in this mass regime; photoevaporation has been shown to be important for such low-density low-mass planets (e.g., Lopez & Rice 2016; Owen & Wu 2017).

In either case, it would be useful to search for multiple populations in this region of the M–R diagram, to test if the bimodal marginal radius distribution unearthed by Fulton et al. (2017) and verified by Hsu et al. (2018) extends to the joint radius and mass space. While there is not currently evidence for a bimodal mass and radius distribution, our approach of calculating the M–R and R–M relations by first estimating the joint distribution is capable of finding and quantifying such bimodalities. Before the joint distribution can be interpreted as an astrophysical result, however, both selection effects discussed in Section 3.1 would need to be incorporated into the density estimation.

Given the current data set, a lognormal distribution is a reasonable approximation for the conditional distribution of mass given radius in most radius regimes (see Figure 5). On the other hand, there is significant skewness in the conditional distributions for log-scaled radius, given log-scaled mass (see Figure 7). This result demonstrates that adopting the more flexible model was warranted, and that different mass regimes have different composition distributions. In particular, the super-puffy planets found by TTV analyses are visible as the skewed tail in the 10 M_⊕ conditional distribution.

As seen in Figure 5, mass predictions are the best constrained for planets at 3 R_⊕: the conditional distribution for mass at 3 R_⊕ is the most robust to the bootstrapping scheme we used to measure the uncertainty in these distributions. As a result, RV follow-up of TESS will yield the most scientific return for 1–2 R_⊕ planets and above ∼4 R_⊕ where there are intrinsically fewer planets.

While few Kepler target stars were sufficiently bright to enable follow-up, a larger proportion of K2 planet-hosting stars are not. One would hope that these brighter stars would enable more mass detections of smaller planets, yet we see in Figure 8 that there is significant bias in the K2 planets planetary masses, especially for the planets with ≤7 R_⊕ radii. This serves as a word of warning for those who aim to perform mass–radius analyses with TESS data: to build robust M–R relations, all mass measurements, even upper limits, must be published and incorporated into the analysis (see Burt et al. 2018 for an in-depth study of this effect). In addition, the second type of selection effect described in Section 3 must be quantified and published as part of major radial velocity follow-up programs.

Zoom In Zoom Out Reset image size

The Bernstein polynomials model requires specifying an upper bound and a lower bound for both mass and radius. Although in theory one could choose the upper bound and the lower bound to be any values, in practice setting the upper bounds and the lower bounds to different values will cause the results to be different from the results in Figures 1 and 3. In particular, when there are no data in a region, the Bernstein polynomials model will fit a smooth line toward the overall mean. This happens when the upper bound is chosen to be too large. To illustrate this, we offer Figure 10 as a comparison to Figure 1.

Zoom In Zoom Out Reset image size

In Figure 10, we set the upper bound of masses to be 8 R_⊕, which is the same value used in WFR16's paper, instead of 6.7 R_⊕. One can immediately tell that the region between 6.7 R_⊕ and 8 R_⊕ has a downward trend. This is because the overall mean of the data set's planet masses is smaller than the mass at 6.7 R_⊕, which is where the last observation lies. Because the model tends toward the overall mean when there is no data, a downward trend is produced. This trend might not provide a correct astrophysical interpretation of the M–R relation in this region, and so we strongly caution against extrapolating this nonparametric relation. That said, the prediction intervals are very wide in this region, which suggests that there is significant variability in the mean relation, which would be rectified with more observations.

Choosing the degrees d and d' is somewhat computationally expensive due to the cross-validation method's searching for the optimal d and d' from a large number of potential values. To decrease the computation time, one can force d = d'. We adopt this approach to choose d and then re-estimate the parameters in the model. Figures 10 and 11 are the results produced by choosing d = d'. The cross-validation method chooses d = d' = 35 for the first data set and d = d' = 55 for the second. Comparing these two figures to Figures 1 and 3, respectively, the results do not vary too much.

Zoom In Zoom Out Reset image size