Mass-Radius relationship for M dwarf exoplanets: Comparing nonparametric and parametric methods

TL;DR

This study compares parametric and nonparametric methods to analyze the Mass-Radius relationship of exoplanets around M dwarfs, revealing differences in predicted masses and providing a new Python tool for future research.

Contribution

It introduces a nonparametric approach to modeling the M-R relationship for M dwarf exoplanets and releases a Python code, MRExo, for broader application and bias reduction.

Findings

Nonparametric method predicts smaller masses for small and large planets compared to Kepler-based relations.

Mass distribution for 3 Earth Radii planets is similar across datasets.

The approach reduces bias by not assuming a specific functional form.

Abstract

M dwarfs, though the most abundant star in the galaxy, form only a small subset of stellar hosts with exoplanets with measured radii and masses. In this paper we analyze the Mass-Radius (M-R) relationship of planets around M dwarfs using M-R measurements for 24 exoplanets. In particular, we apply both parametric and nonparametric models and compare the two different fitting methods. We also use these methods to compare the results of the M dwarf M-R relationship with that from the Kepler sample. Using the nonparametric method, we find that the predicted masses for the smallest and largest planets around M dwarfs are smaller than that given by the relation fit to the Kepler data, but that the distribution of masses for 3 Earth Radii planets does not substantially differ between the two datasets. With future additions to the M dwarf M-R relation from TESS and instruments like the Habitable Zone Planet Finder, we will be able to characterize these differences in more detail, which will help illuminate the process of planetary formation and evolution around these stars. We release a publicly available Python code called MRExo which uses the nonparametric algorithm introduced by Ning et al. (2018) to fit the M-R relationship. Such a nonparametric fit does not assume an underlying power law fit to the measurements and hence can be used on samples spanning large mass and radii ranges. Therefore by not assuming a functional form, the fit is less biased. MRExo also offers a tool to predict mass from radius posteriors, and vice versa. This functionality can help inform observational strategies for radial velocity campaigns, such as TESS follow-up studies, as well as predict radii with microlensing planet masses.

Figures

Tables

	$f(m,r|\text{w},d,d^{\prime },M^{obs},R^{obs},{\sigma }_M^{obs},{\sigma }_R^{obs})$
	$=\prod _{i=1}^{n}f\left(M_i^{obs}|m,{\sigma }_{M_i}^{obs}\right)f\left(R_i^{obs}|r,{\sigma }_{R_i}^{obs}\right)\times f(m,r|w,d,d^{\prime })$		(3)
	$=\prod _{i=1}^n\sum _{k=1}^d\sum _{l=1}^{d^{\prime }}{w}_{kl}\frac{{\beta }_{kd}\left(\frac{m-\underset{¯}{M}}{\overline{M}-\underset{¯}{M}}\right)}{\overline{M}-\underset{¯}{M}}𝒩\left(\frac{M_i^{obs}-m}{{\sigma }_{M_i}^{obs}}\right)\times \frac{{\beta }_{l{d}^{\prime }}\left(\frac{r-\underset{¯}{R}}{\overline{R}-\underset{¯}{R}}\right)}{\overline{R}-\underset{¯}{R}}𝒩\left(\frac{R_i^{obs}-r}{{\sigma }_{R_i}^{obs}}\right),$		(4)

Table 2. Table 1: Two-sample Comparison Tests between the Parametric and Nonparametric fits to the M dwarf sample. Note that these are relative p 𝑝 p -values based on generated datasets of 100 planets, not absolute p 𝑝 p -values (see § 5.3 for discussion).

Radius	Kolmogorov-Smirnov		Anderson-Darling
$R_{\oplus }$	Statistic	$p$-value	Statistic	$p$-value
1	0.26	1.3 $\times$ 10-3	8.30	10-4
3	0.30	2 $\times$ 10-4	8.65	10-4
10	0.15	0.1511	1.87	0.1046

Table 3. Table 2: 2-sample comparison tests between the nonparametric fit to the M dwarf and Kepler samples. Note that these are relative p 𝑝 p -values based on generated datasets of 100 planets, not absolute p 𝑝 p -values (see § 5.3 for discussion).

Radius	Kolmogorov-Smirnov		Anderson-Darling
$R_{\oplus }$	Statistic	$p$-value	Statistic	$p$-value
1	0.47	10-4	24.2	10-4
3	0.16	0.1057	2.03	0.0824
10	0.42	10-4	14.5	10-4

Table 4. Table A1: Masses and Radii of M dwarf planets used in this work.

Planet	Mass	${\sigma }_M$	Radius	${\sigma }_R$	Mass	Radius
Name	($M_{\oplus }$)	($M_{\oplus }$)	($R_{\oplus }$)	($R_{\oplus }$)	Reference	Reference
HATS-6 b	101.3878	22.2481	11.187	0.213	Hartman et al. (2015)	Hartman et al. (2015)
GJ 1214 b	6.26125	0.85814	2.847	0.202	Harpsøe et al. (2013)	Harpsøe et al. (2013)
K2-22 b${}^{\text{<a><span>7</span></a>}}$	444.962	0	2.5	0.4	Sanchis-Ojeda et al. (2015)	Sanchis-Ojeda et al. (2015)
K2-18 b	7.96	1.91	2.38	0.22	Cloutier et al. (2017)	Cloutier et al. (2017)
LHS 1140 b	6.98	0.89	1.727	0.032	Ment et al. (2018)	Ment et al. (2018)
NGTS-1 b	258.078	22.40702	14.908	5.268	Bayliss et al. (2018)	Bayliss et al. (2018)
K2-137 b777 Upper limits	158.915	0	0.89	0.09	Smith et al. (2018)	Smith et al. (2018)
K2-33 b${}^{\text{<a><span>7</span></a>}}$	1175.971	0	5.04	0.355	Mann et al. (2016)	Mann et al. (2016)
LHS 1140 c	1.81	0.39	1.282	0.024	Ment et al. (2018)	Ment et al. (2018)
GJ 436 b	23.1	0.8	4.191	0.109	Turner et al. (2016)	Turner et al. (2016)
GJ 1132 b	1.66	0.23	1.43	0.16	Bonfils et al. (2018)	Southworth et al. (2017)
Kepler-32 b	9.4	3.35	2.231	0.072	Hadden & Lithwick (2014)	Berger et al. (2018)
Kepler-32 c	7.7	4.4	2.112	0.071	Hadden & Lithwick (2014)	Berger et al. (2018)
GJ 3470 b	13.9	1.5	4.57	0.18	Awiphan et al. (2016)	Awiphan et al. (2016)
TRAPPIST-1 b	1.017	0.1485	1.121	0.0315	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 c	1.156	0.1365	1.095	0.0305	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 d	0.297	0.037	0.784	0.023	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 e	0.772	0.077	0.91	0.0265	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 f	0.934	0.079	1.046	0.0295	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 g	1.148	0.0965	1.148	0.0325	Grimm et al. (2018)	Grimm et al. (2018)
TRAPPIST-1 h	0.331	0.0525	0.773	0.0265	Grimm et al. (2018)	Grimm et al. (2018)
Kepler-45 b	159	19	11.634	0.362	Southworth (2012)	Berger et al. (2018)
Kepler-138 b	0.066	0.048	0.629	0.0275	Jontof-Hutter et al. (2015)	Berger et al. (2018)
Kepler-138 c	1.97	1.516	1.519	0.0975	Jontof-Hutter et al. (2015)	Berger et al. (2018)
Kepler-138 d	0.64	0.5305	1.323	0.043	Jontof-Hutter et al. (2015)	Berger et al. (2018)
Kepler-26 b	5.02	0.66	3.191	0.095	Hadden & Lithwick (2017)	Berger et al. (2018)
Kepler-26 c	6.12	0.71	2.976	0.253	Hadden & Lithwick (2017)	Berger et al. (2018)

Equations

$$\log \frac{M}{M_{\oplus }}\sim \mathcal{N}\left(\mu =C+\gamma \log \frac{R}{R_{\oplus }},\sigma ={\sigma }_M\right),$$

$$f(m,r|w,d,d^{\prime })$$

$$=\sum _{k=1}^d\sum _{l=1}^{d^{\prime }}{w}_{kl}\frac{{\beta }_{kd}\left(\frac{m-\underline{M}}{\overline{M}-\underline{M}}\right)}{\overline{M}-\underline{M}}\frac{{\beta }_{l{d}^{\prime }}\left(\frac{r-\underline{R}}{\overline{R}-\underline{R}}\right)}{\overline{R}-\underline{R}},$$

Full text

planets and satellites: composition

††software: MRExo (Kanodia, 2019), Astropy (Robitaille et al., 2013), Numpy (Oliphant, 2006–), Matplotlib (Hunter, 2007), SciPy (Oliphant, 2007), iPython (Pérez & Granger, 2007), Enthought Canopy Python.

1 Introduction

In the Galaxy, M dwarfs are the most common type of star ($\sim$ 75%; Henry et al., 2006). With the launch of the Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2014), the hope is that we will soon discover hundreds of exoplanet candidates around them. While the discovery of these planets is interesting in itself, the comparison between them and the Kepler planets provides insight into the differing formation pathways of planets around M and FGK stars (Lissauer, 2007). For example, the pre-main-sequence lifetime of the star varies with its mass, which has an effect on the planetary migration process during its formation. Young M dwarfs are extremely active and exhibit high-intensity XUV radiation, which affects the inner planets and can potentially strip away the atmospheres of gaseous planets to leave rocky cores (Owen & Wu, 2017; Owen & Lai, 2018, show this for Kepler planets). There is also empirical evidence from Kepler (Dressing & Charbonneau, 2015; Gaidos et al., 2016) and radial velocity (RV) surveys (Bonfils et al., 2013) that suggests that the mass and radius distribution of planets is not identical for M and FGK dwarfs. These suggest that different physical processes may be at play and pose a number of questions: Do smaller planets around M dwarf have more rocky compositions (Mulders et al., 2015a, b)? Is planet formation more efficient around M dwarfs (Dressing & Charbonneau, 2015; Ballard & Johnson, 2016; Ballard, 2018)? If so, how does planet formation impact exoplanet chemical composition?

Probabilistic mass–radius (M-R) relationships provide us with an empirical window into these questions, as they are closely related to distributions of exoplanet compositions. They also have very practical uses. For example, efficient planning of TESS follow-up RV observations requires an estimate of the planetary mass given a planetary radius. Future microlensing space missions like WFIRST (Green et al., 2012) will produce hundreds of exoplanets with the inverse problem of having a mass but not a radius. The M-R relationships can be used to predict one quantity from the other. For this purpose, one needs a model for the M-R relationship that best balances the trade-off between the prediction’s variance (i.e. the width of the range of possible masses for a given planet) and its bias in the predicted values (i.e. the difference between the true mass and the mean predicted mass).

Using transit spectroscopy due to the lower stellar brightness and relatively large planet-to-star radius ratio, M dwarfs will also offer potential targets for atmospheric characterization of Earth-like habitable zone planets. It has also been shown that the spot-induced RV jitter is reduced in the near infrared (NIR) (Marchwinski et al., 2015), which is where the spectral energy distributions for these stars peak. For these reasons, many more M dwarf planets will have their masses measured in the near future by instruments such as the HPF (Mahadevan et al., 2012), CARMENES (Quirrenbach et al., 2016), NIRPS (Wildi et al., 2017), IRD (Kotani et al., 2018), SPIRou (Artigau et al., 2014), iSHELL (Cale et al., 2018), GIANO (Claudi et al., 2017), and NEID (Schwab et al., 2016). Here we set the stage for these future datasets by assessing the dependence of the M-R model choices on mass and radius predictions, especially as a function of stellar type.

Substantial past efforts have been put toward studying M-R relations in recent years. A summary of this is presented in Ning et al. (2018). Several of the widely used M-R relationships include those of Weiss & Marcy (2014); Wolfgang et al. (2016); Chen & Kipping (2017), the latter of which also introduced a publicly available Python package called Forecaster. The M-R model underlying Forecaster uses a broken power law to fit the M-R relationship across a vast range of masses and radii, in recognition of potential changes in the physical mechanisms responsible for the planetary formation at different mass regimes. However, as has been discussed in Ning et al. (2018) such a restrictive parametric model can portray an incomplete picture, since we do not know the true functional form of the underlying relationship, such as whether it is a power law to begin with. Conversely, a nonparametric model offers more flexibility in the fit, which can be advantageous when the goal is to obtain predictions that best reflect the existing dataset. Ning et al. (2018) introduced a nonparametric model for the M-R relationship which uses Bernstein polynomials, a series of unnormalized Beta probability distributions. We apply this model to M dwarf exoplanets in preparation for future, larger datasets of small planets around small stars, which are less likely to be fully described by highly parametric models.

While adapting the methodology of Ning et al. (2018) to a new dataset, we have developed and are now offering a publicly available Python package called MREXo, inspired by the useful community tool Forecaster. Not only can MREXo be used as a predicting and plotting tool for our M dwarf and Kepler M-R relationships, it can also be used to fit an M-R relationship to any other dataset (but see §6 for a discussion about the minimum dataset size at which a nonparametric fit becomes useful). This makes it a powerful tool for exoplanet population studies and a probe for potential differences in composition across samples. This package uses the open-source tools of Python deployed with parallel processing for efficient computation. It also offers a fast predictive tool for the M dwarf and Kepler sample M-R fits used in this paper, so that either mass or radius can be predicted from the other.

The rest of the paper is structured as follows. In §2 we discuss the input dataset, and in §3 we discuss the parametric and nonparametric fitting process and algorithms followed. In §4, we describe the Python package MRExo that we release along with this paper. Then, in §5 we discuss the results of the fits, where we compare the parametric fit with the nonparametric and the M dwarf M-R relationship with that from the Kepler exoplanet sample. In §6, we explain the simulation performed to test the efficacy of the nonparametric method. We end with a discussion in §7 and conclude in §8.

2 Data Set

Fitting an M-R relationship requires a sample set with confirmed mass and radius measurements222The nonparametric framework employed here cannot handle measurement upper limits, since the MLE method used in Ning et al. (2018) does not allow censored data; hence, planets with only upper limits are excluded from the results presented in §5.2. Adapting this methodology to include upper limits is an area for future work.. The mass and radius values for the exoplanets used in this work are obtained from the NASA Exoplanet Archive, which we last accessed on 14th December 2018 (Akeson et al., 2013). Figure 1 shows the 24 M-R data points that we have used, colour-coded by host star temperature. The mean values for mass and radius, along with their respective measurement uncertainties, are shown in Table A1 in the Appendix. We include mass value estimates from both RV and Transit Timing Variations (TTVs). This sample is hereafter referred to as the M dwarf dataset in this paper. We specifically chose to omit the three planets discovered by direct imaging with both radius and mass constraints, as these planets have substantially larger orbital separations than the other planets in our sample. Furthermore, the directly imaged planets have their masses and radii modelled and not directly measured.

To limit ourselves to the M dwarfs, we restrict our host star sample to $T_{\rm{eff}}$ $<$ 4000 K. To exclude brown dwarf companions, we restrict ourselves to planetary masses ($M_{\rm{p}}<$ 10 $M_{\rm{J}}$). We chose to use the most recent mass and radius values from the Exoplanet Archive rather than the default values. We manually verified the reported values from the Exoplanet Archive with the literature references. For the TRAPPIST-1 system, we used the Grimm et al. (2018) values for the mass and radius measurements. Note that of the 24 exoplanet measurements in our sample, seven belong to the TRAPPIST-1 system and thus our M-R joint distribution is heavily influenced by the TRAPPIST-1 system in the small planet ($R$ $<$ 2 $R_{\oplus}$) regime. For purposes of comparison, we also use the Kepler dataset from Ning et al. (2018) to compare our M dwarf results with those from a larger Kepler sample, which consists of 127 M-R measurements and does not have any $T_{\rm{eff}}$ restrictions on the host stars. This sample is hereafter referred to as the Kepler dataset in this paper.

3 Fitting the Data set

To fit these data sets, we use two recently proposed methods: a parametric method proposed by Wolfgang et al. (2016), and a nonparametric method by Ning et al. (2018). The two methods and their M-R fits are detailed and compared in subsequent sections. We perform the following analysis.

Comparing the parametric and nonparametric results for the M dwarf sample. 2. 2.

Comparing the M dwarf and Kepler sample M-R relations using the nonparametric model. 3. 3. 4. 4.

Performing a simulation study to show the efficacy of the nonparametric method as (a) the sample size increases, and (b) the intrinsic astrophysical dispersion in the sample increases (assuming a power law with an intrinsic dispersion).

We compare the parametric method with the nonparametric technique because the planetary masses and radii of the M dwarf sample currently appear to resemble a power law. This may not remain the case as we accumulate more data for M dwarf planets; indeed, with over 100 mass and radius measurements, the Kepler dataset is not currently well described by a power-law model (Ning et al., 2018). Additionally, nonparametric methods offer less biased predictions, as long as the dataset they were fit to is representative of the underlying distribution. In contrast, nonparametric methods often perform poorly for small data sets and can easily overfit the data, while parametric methods are easier to implement and are valuable for gaining an intuitive understanding of the dataset. Acknowledging both the simplicity and easy insight from a parametric model and the improved predictive capabilities of a nonparametric model for larger datasets, we compute both and compare the outputs for the current M dwarf sample (see §5). In addition, we run simulations to assess at what dataset size does the nonparametric technique work well (we also test the fitting as the intrinsic dispersion in the dataset increases).

We note that the nonparametric method cannot handle upper limits in its current implementation. In the TESS era of planet discoveries, we will soon have many exoplanet candidates with just upper limits, and so this is a clear area for future statistical development of this nonparametric approach. For the time being, we incorporate the upper limits in the parametric approach, which can accommodate censored data like upper limits.

3.1 Parametric method: Hierarchical Bayesian Modelling

The parametric method we use is heavily borrowed from the probabilistic hierarchical Bayesian model (HBM) introduced in Wolfgang et al. (2016). This model consists of a power law to describe the mean planet mass as a function of radius plus a normal distribution around that line to describe the intrinsic astrophysical scatter in planetary masses at a specific radius. The model we employ here is very similar to that of Wolfgang et al. (2016, Equation 2), except that we use a lognormal distribution to describe the intrinsic scatter. This change causes the intrinsic dispersion to be constant in log space, which Chen & Kipping (2017) showed is a better descriptor of the scatter over a larger mass range than was considered in Wolfgang et al. (2016). The parametric model used here is therefore:

[TABLE]

where $C$ is the normalization constant of the mean power law (once it is converted to the linear mass scale), the $\gamma$ is the power-law index, and the $\sigma_{M}$ is the intrinsic scatter in terms of log($M$). The $\sim$ symbol implies that the masses are drawn from a probabilistic distribution (as opposed to an $=$ sign for a deterministic model); this distribution is what models the astrophysical scatter in planetary masses. Here it is parameterized as a Gaussian centered on a line, which when converted to linear scale is a lognormal distribution centered on a power law. Evaluating this model within a hierarchical framework allows the upper limits in our dataset to be incorporated into the inference of $C$, $\gamma$, and $\sigma_{M}$, along with an arbitrary measurement uncertainty for individual data.

3.2 The Nonparametric Method – Bernstein polynomials

As explained in the Introduction, we also employ the nonparametric model introduced by Ning et al. (2018) to fit M-R relationships, with an eye toward future M dwarf dataset sizes that will likely benefit from a more flexible approach. The nonparametric method fits the joint M-R distribution, and hence can treat either mass or radius as the independent variable, and can be used to predict one from the other. This method uses a sequence of Bernstein polynomials as the basis functions to fit a nonparametric M-R relationship; when normalized, these Bernstein polynomials are identical to beta probability distributions. Hence, the model is equivalent to a mixture of unnormalized Beta probability distributions. We translate the code presented in Ning et al. (2018) from R to Python and release it in a publicly available package called MRExo; this code is further discussed in §4.

The nonparametric approach uses these Bernstein polynomials as the basis functions to fit the joint distribution of masses and radii f(m,r). This is detailed in §2 of Ning et al. (2018):

[TABLE]

where $w_{kl}$ is the $kl$-th element of the matrix w, which is a set of weights corresponding to individual Bernstein polynomials. Each weight is positive, and the sum of individual weights equals unity. The weights describe how much each term in the series contributes to the joint distribution. The degrees for the polynomials in the mass and radius dimensions are represented by d and $d^{\prime}$, respectively. Here m and r depict the continuous variables which represent mass and radii. $\underline{M}$, $\overline{M}$, $\underline{R}$ and $\overline{R}$ represent the lower and upper bounds in mass and radius, respectively. To fit for the weights ($w_{kl}$), we use maximum-likelihood estimation (MLE). The likelihood (Equation 4) includes measurement uncertainties as normal distributions; hence, Equation 4 modifies subsection 3.2 to introduce a convolution of the normal and Bernstein polynomials. The measured values for mass and radii are assumed to be drawn from a normal distribution centered on the true value, with a standard deviation equal to the measurement uncertainty. This produces the joint distribution shown in Equation 4, where $M_{i}^{obs}$, $\sigma_{M_{i}}^{obs}$ are the mass observations and their uncertainties, and likewise $R_{i}^{obs}$, $\sigma_{R_{i}}^{obs}$ for radius. Here $\mathcal{N}$ denotes a normal distribution. After optimization for the weights using the MLE method via the Python package SciPy (Oliphant, 2007), the joint distribution for the M dwarf dataset can be plotted (see Figure 3 b).

A key consideration in fitting this nonparametric model to a dataset is identifying the optimum degree for the Bernstein polynomial series. We use the cross-validation method as explained in Ning et al. (2018) to find that the optimum values for d and $d^{\prime}$ are both 17. In total, there are $\textit{d}\times$ $d^{\prime}$ weights that serve as our ‘parameters’ to fit for. While $17^{2}=289$ parameters may seem excessive for a dataset size of 24 planets, in practice, there are only five non-zero weights in the series (see Figure 2 for a pictorial representation of this). This highlights one of the advantages of using Bernstein polynomials as our basis functions: estimating the weights in this series is self-regularizing, meaning that the smallest number of nonzero coefficients is automatically found. Additionally, using Bernstein polynomials efficiently reduces the number of nonzero free parameters; if one instead used a mixture of Gaussians to fit the joint distribution, that fit would require at least three times as many free parameters (amplitude, mean, and standard deviation per component, rather than just the polynomial coefficients).

After finding the weights via MLE, we repeat the process using the bootstrap method, which helps to account for the variation of the parameters in the model. We do this by resampling the dataset with replacement and running the fitting routine again for each bootstrap. In regions without data points, the Bernstein polynomials revert to the overall mean of the distribution.

We found that the Bernstein polynomials may behave counter-intuitively at the boundaries of the joint distribution (see Appendix A in Ning et al. (2018)). To address this issue, we fix the first and last row and column of Bernstein polynomials to have zero weights ($w_{0j},w_{dj},w_{i0},w_{id^{\prime}}$), and as such, they do not contribute to the joint distribution. This was not done by Ning et al. (2018) since their Kepler sample had enough samples near the boundaries to constrain the fit.

4 MRExo

In this section, we shall discuss an important contribution made in this paper. We translate Ning et al. (2018)’s R code 333https://github.com/Bo-Ning/Predicting-exoplanet-mass-and-radius-relationship. into a publicly available Python package called MRExo444https://github.com/shbhuk/mrexo. This can be easily installed using PyPI and has extensive documentation and tutorials to make it easy to use.

The MRExo package offers tools for fitting the M-R relationship to a given data set. In this package, we use a cross-validation technique to optimize for the number of degrees. We then fit the joint distribution ($\S$ 3.2) to the sample set; this can then be marginalized to obtain the conditional distribution, which we can use to predict one variable from the other. We bootstrap our fitting procedure to estimate the uncertainties in the mean M-R relation. Further, MRExo is equipped with dedicated and easy-to-use functions to plot the best-fit conditional M-R relationships, as well as the joint M-R distribution555We caution the user not to overinterpret the joint distribution; in the first version of this software we have made no attempt to correct for detection and selection bias, which is needed before the joint distribution can be interpreted as an M-R occurrence rate. We show the joint distribution in Figure 3 solely to illustrate how the behavior of the conditionals relates to the joint.. Crucially, MRExo also predicts mass from radius, and radius from mass. For example, in the case of planets discovered using the transit method, the feasibility of an RV follow-up campaign can be evaluated by predicting the estimated mass and its confidence intervals given the measured radius and its uncertainty. Another feature of this package is that it can accommodate radius (or mass) posterior samples from separate analyses, which are then used to compute the posterior predictive distribution for mass . Along with the MRExo installation, the results from the M dwarf sample dataset from this paper and the Kepler exoplanet sample from Ning et al. (2018) are included.

The degree of the Bernstein polynomials ($d$) approximately scales with the sample size ($N$). Since the number of weights goes as $d^{2}$, the computation time involved in the fitting a new M-R can soon start to become prohibitive. Therefore, we also parallelize the fitting procedure and the bootstrapping algorithm. As an example, the M dwarf sample ($N=24$; $d=17$) took about 2 minutes to perform cross-validation, fit a relationship, and do 100 bootstraps on a cluster node with 24 cores and 2.2 GHz processors. The Kepler sample ($N=127$; $d=55$) took about 36 hr for the cross-validation, fitting, and 48 bootstraps. We realize that the fitting computation time will start to become prohibitive as the sample size increases $\gtrsim 200$; therefore, we plan to optimize the code further by benchmarking, optimizing floating point operations, and correcting the precision requirements in the integration step. However, this time-intensive step of cross-validation and fitting is only necessary if the user needs to fit their own relationships. To run the prediction routine on the preexisting M dwarf or Kepler samples is fairly quick and takes less than a second for a prediction. In order to do a large number of predictions as part of a larger pipeline or simulation, the user can also generate a lookup table that makes the calculations even faster (the function to generate and use the lookup table is provided with the package).

5 Results

5.1 Parametric fit results

As discussed in §3, there is a trade-off between the flexibility and lower bias of nonparametric models and the interpretability and lower predictive variance of parametric models. At present, the M dwarf dataset only consists of 24 planets, and their masses and radii seem to be well approximated by a power law. As such, we fit the parametric model described by Equation 1 to the same M dwarf data set, to serve as a basis for comparison to the nonparametric results. Currently, the fit looks to be fairly reasonable by eye (see Figure 3e), but this is not guaranteed to be the case as more measurements are obtained (indeed, it is clear from the fit to the Kepler data, which consist of $>100$ planets, that a single power law is not a sufficient model for that data; see Figure 3f).

The best-fit parameter values for the model described in Equation 1 are $C=-0.130^{+0.081}_{-0.055},\gamma=2.14^{+0.11}_{-0.16}$, and $\sigma_{M}=0.184^{+0.077}_{-0.021}$; these values were found by identifying the mode of the joint three-dimensional hyperparameter posterior, and the error bars represent the marginal central 68% credible intervals. This power-law slope is steeper than that found by Wolfgang et al. (2016) for nearly all of the datasets they consider (most lie within $1.3<\gamma<1.8$). There are two possible explanations for this. First, the dataset used by Wolfgang et al. (2016) — and, in fact, nearly all previous M-R results except Neil & Rogers (2018) — is dominated by FGK dwarf planet hosts. This result could therefore be driven by intrinsic differences between the planets that form around Sun-like stars and those that form around M dwarfs. However, Wolfgang et al. (2016) also use only a limited radius range for their dataset ($R<4$ $R_{\oplus}$ or $R<8$ $R_{\oplus}$), as their focus was on super-Earths. As shown by Ning et al. (2018), the slope for $1$ $R_{\oplus}<R<4$ $R_{\oplus}$ is shallower than that for $4$ $R_{\oplus}<R<10$ $R_{\oplus}$, and so one would expect that a single power-law fit to the entire radius range would result in a steeper slope than that fit just to $R<4$ $R_{\oplus}$. This is corroborated by the fact that the slopes fit to the $R<8$ $R_{\oplus}$ radius range are steeper than those fit to the $R<4$ $R_{\oplus}$ range. Additionally, Chen & Kipping (2017) found that the segment of their broken power-law relation that most closely matches radius range of our M dwarf dataset ($1$ $R_{\oplus}\lesssim R\lesssim 11$ $R_{\oplus}$) has a power-law slope of $\frac{1}{0.59}=1.7$. We note that the fact that there is a noticeable difference between the slopes fit to different radius ranges is evidence that a more flexible relation, such as the nonparametric one developed in Ning et al. (2018) is needed.

Another notable difference between the parametric fit for the M dwarf planets and the FGK planet sample from Wolfgang et al. (2016) is the predicted mass for a 1 $R_{\oplus}$ planet. This information is contained in $C$, the mean power-law normalization constant. In linear units, the best-fit $C$ for the M dwarf dataset is $10^{-0.13}=0.74$ $M_{\oplus}$. This is significantly lower than the power-law constants found in Wolfgang et al. (2016): it is inconsistent with the $R<8$ $R_{\oplus}$ fit ($C=1.5$ $M_{\oplus}$) at 4$\sigma$, and with the $R<4$ $R_{\oplus}$ fit ($C=2.7$ $M_{\oplus}$) at 7$\sigma$. At least two effects are contributing to this difference. First, very few planets around FGK dwarfs with $R<1$ $R_{\oplus}$ have had their masses measured and constrained to be $M<1.5$ $M_{\oplus}$. Because of the publication bias quantified in Burt et al. (2018), wherein planetary masses are required to reach a certain significance threshold to be published, small, low-mass planets are more likely to be left out of the FGK dwarf datasets. This causes the relation to be fit to the more massive small planets that do appear in the literature, which in turn causes $C$ to be high. Conversely, the M dwarf dataset is strongly affected by the presence of the TRAPPIST-1 planets (Gillon et al., 2017), whose masses are reported by Grimm et al. (2018) to be both small and quite precise. This one analysis of a single system dominates the M dwarf dataset at $0.7$ $R_{\oplus}$ $<R<1.2$ $R_{\oplus}$ and therefore the best-fit value of $C$ — a caveat that all potential users of these M dwarf relations should keep in mind.

To determine whether the differences highlighted above are due to the differences in the host star population, models, or considered radius ranges, we also fit this simple parametric model to the Kepler dataset from Ning et al. (2018). While we do not expect this model to be a good fit to the data, it offers an apples-to-apples comparison to the above M dwarf parametric results, and its shortcomings highlight the advantages of more flexible nonparametric models. The best-fit parameter values for the Kepler dataset are $C=-0.0250^{+0.093}_{-0.111},\gamma=2.13^{+0.11}_{-0.14}$, and $\sigma_{M}=0.457^{+0.026}_{-0.044}$. Therefore, there is no statistically significant difference between the M dwarf planetary mean power law and the one that would be fit to the planets around FGK dwarfs. More specifically, the two slopes are nearly identical; while the Kepler normalization constant is larger than the M dwarf constant, the difference is at the $1\sigma$ level. This comparison is consistent with the result of Neil & Rogers (2018), which used a sample size of six M dwarf planets. On the other hand, there is a statistically significant difference in the intrinsic scatter for the M dwarf and Kepler datasets, with the Kepler dataset having more variation in planet mass at a given radius. However, this may be at least partially affected by the differences in the sample sizes ($N=24$ vs. $N=127$) and the fact that it takes more time to get significant mass measurements for planets on the low mass side of the intrinsic scatter distribution. As we obtain more masses for transiting M dwarf planets with instruments like HPF, we will be able to test whether the difference in the intrinsic scatter remains statistically significant.

All told, the comparison between the M dwarf and Kepler parametric model fits does not yield very much insight into the differences between the planetary populations. Because the parametric model is not able to capture detailed features of the M-R relation, it may hide some important differences between the two datasets once the M dwarf dataset becomes large enough to warrant detailed analyses of these features. Looking ahead to the M-R dataset we will have by the end of the TESS mission, we apply the nonparametric model to the current M dwarf dataset and perform the comparison between it and the nonparametric fits to the Kepler dataset to serve as a basis for comparison to this future work (see §5.4).

5.2 Nonparametric fit

We fit the M dwarf exoplanet dataset shown in Figure 1 using the nonparametric approach described in §3.2; the results of this fit are displayed in Figure 3 (a) and (b). Using the cross-validation method the optimum number of degrees selected are both equal to 17, giving 17${}^{2}=289$ total weights. We note that our algorithm automatically forces most of the weights to zero to prevent overfitting, especially for small datasets. For this dataset, only five weights are nonzero (see Figure 2).

We also plot the conditional distributions $f(m|r)$ and $f(r|m)$, where the conditional distribution $f(m|r)$ is the ratio of the joint distribution (subsection 3.2) to the marginal distribution of the radius $f(r)$; likewise for $f(r|m)$(see Equations 1 and 2 of Ning et al. (2018)). The conditional distribution spans regions of M-R space where there is no data. This is due to the symmetrical nature of the Bernstein polynomials and indicates that the current dataset is not yet large enough to be effectively described by a nonparametric method; this is further discussed in §6.

From Figure 3(b), we note that the red conditional distribution–$f(m|r)$ is not the same as the blue conditional distribution–$f(r|m)$. This is because when the joint distribution (see Figure 3 (b) and (d)) is marginalized to obtain the conditional, the distribution behaves differently for each axis. The two conditional distributions would be the same for a dataset with very little error and symmetric, equal scatter along the entire M-R relation (see §6 for an illustration of this). The Kepler dataset effectively illustrates how localized areas of large scatter in one dimension can drive differences between $f(m|r)$ and $f(r|m)$. In particular, at log($m$) $\sim 0.5$, there is large scatter in radius, which is fit by the mean M-R relation in $f(m|r)$ but is represented by the distribution (the width of the light shaded region) around that mean relation in $f(r|m)$. While both conditionals were computed from the same joint distribution (Figure 3d), asymmetries such as this in the joint distribution can produce conditionals that qualitatively look very different. We also note that truncating the probability distributions for integration (that is, using finite bounds to compute $f(r)$ and $f(m)$ for the conditionals) also contributes to some differences between the two conditional distributions. This can be further seen in the simulation in Figure 7, where the disparity between the two conditional distributions increases as the dispersion from the power law increases. Therefore, to predict one quantity from the other it is imperative to use the right conditional distribution.

A common concern about nonparametric methods is their ready ability to overfit data. The cross-validation method we adopt to choose $d$ and $d^{\prime}$ is designed to minimize overfitting while maximizing the predictive accuracy of the M-R relation. Both underfitting and overfitting produce high predictive error (that is, the predicted value is far from the true value); because the cross-validation method identifies the degree that minimizes the predictive error, it finds the optimum $d$ that balances the trade-off between the two. That said, we acknowledge that by eye, there are wiggles in the M dwarf M-R relation that do not appear to be supported by the dataset. This is likely due to the sparsity of the dataset in certain radius and mass ranges; where there is no data, the relation tends toward the mean of the closest nonzero term. This is why we perform the simulation study in §6, to find at what size dataset does our nonparametric model effectively describe the conditionals.

5.3 Comparing Nonparametric vs. Parametric fitting

We plot the conditional probability density functions (CPDFs) of planetary masses for planets with $R=1,$ $3$, and $10$ $R_{\oplus}$ in Figure 4. These CPDFs show the distribution of planetary masses that would be predicted for planets at those radii. Assessing the differences visually, the parametric and nonparametric M dwarf fits have the same median but different spreads. This difference is expected due to the bias vs. variance trade-off in parametric vs. nonparametric modeling: if the parametric model is a poor description of the data, the nonparametric fit will be less biased but have a higher variance than the parametric estimator. Because the medians of the CPDFs are similar between the nonparametric and parametric models, it appears that, with this current dataset, a power law is a decent fit to the data. That said, we emphasize that this will likely not remain the case as the M dwarf dataset grows, as we have seen for the Kepler dataset.

To quantify the differences between the mass predictions produced by these two model choices, we perform the two-sample Kolmogorov-Smirnov (K-S; Kolmogorov, 1933; Smirnov, 1948) and Anderson-Darling (AD; Anderson & Darling, 1952; Scholz & Stephens, 1987) tests on samples drawn from these CPDFs (see Table 1). These statistical tests assess whether we can reject the null hypothesis that two datasets were drawn from the same distribution, with the K-S test being more sensitive to the location of the median and the AD tests more sensitive to differences of the tails.

Importantly, the reported significance of these tests depends on the size of the datasets, with small differences between the samples becoming more significant as the datasets grow. Since we have the functional form of these CPDFs from our fitting procedure, we have a choice in how many samples we can draw from them, and therefore how significant we report them to be666For every comparison in Tables 1 and 2, we use 100 samples from each CPDF. We also note that the assignment of a $p$-value to a given statistic depends on a theoretical statistical distribution for the variation in that statistic given randomly drawn datasets from the same distribution. To test this theoretically assigned $p$-value, we generated reference distributions of K-S and AD statistics based on 20,000 randomly drawn datasets (each of which also contains $N=100$ planets) from the CPDFs, compared against themselves for a total of 10,000 comparisons (and 10,000 values for the K-S and AD statistics). We then identified at what quantile is the median of the nonparametric vs. parametric K-S statistic distribution with respect to the distribution of KS statistics produced by comparing the M dwarf nonparametric CPDFs against themselves. These values are what is reported in Table 1 (and an analogous comparison but between the Kepler and M dwarf nonparametric CPDFs with respect to the self-compared nonparametric M dwarf CPDFs is reported in Table 2). It turns out that this more careful approach is in good agreement with the theoretical p-values provided by the R package kSamples, built from the work of Scholz & Stephens (1987). If the median of the nonparametric vs. parametric K-S statistic distribution fell completely outside that of the nonparametric self-compared KS distribution, the p-value is reported as $<10^{-4}$, as we only generate 10,000 pairs of datasets for these K-S and AD statistic distributions.. As such, the results in Tables 1 and 2 should be interpreted on a comparative basis, not on an absolute basis: rather than concluding that the CPDFs for 1 $R_{\oplus}$ planets are different on a statistically significant level while the CPDFs for 10 $R_{\oplus}$ are not, the correct interpretation is that it will take a smaller number of mass measurements of planets at 1 $R_{\oplus}$ to distinguish between the parametric and nonparametric models than it will for planets at 10 $R_{\oplus}$ (see Table 1). Along the same vein, an even smaller number of mass measurements would be needed to distinguish between the parametric and nonparametric models for 3 $R_{\oplus}$ planets. Because we only have 24 planets in our M dwarf dataset in total, let alone at any single radius, these model comparisons will need to be reassessed in the future with more data.

5.4 Comparing M dwarf M-R vs. Kepler M-R

We also seek to compare the M-R relationship from a Kepler exoplanet sample to an M dwarf exoplanet sample. As discussed in §1, there is empirical evidence that the planetary radius and mass distributions differ between M and FGK dwarfs. This is further illustrated empirically by our fit to the M dwarf and Kepler sample (Figure 3): there are visual differences between the joint and conditional distributions for the M dwarfs and Kepler. Additionally, Figure 4 shows that the median value of the CPDFs of Kepler and M dwarf masses do not coincide for the smallest and largest planets; this is true whether the nonparametric M dwarf relation or the parametric M dwarf relation is used. Therefore, if the Kepler M-R relationship was used to predict the mass of a transiting M dwarf exoplanet like TRAPPIST-1 b, it would produce mass predictions that were too large, on average. Conversely, if used for a nontransiting planet like Proxima-b to predict its radius from its mass, the prediction would be erroneous. This illustrates that the conditional density functions are different for the two samples. We also perform the K-S and AD tests to quantify differences between the M dwarf and Kepler CPDFs (see §5.3 for discussion about these tests). We find that it will take $<<100$ planets to distinguish between the M dwarf and Kepler mass predictions at 1 and 10 $R_{\oplus}$, while even a dataset of 100 3 $R_{\oplus}$ planets does not illuminate any statistically significant difference between the mass predictions. Therefore, the M dwarf M-R relation could be most different from the FGK dwarf M-R relation at the smallest and the largest planet sizes.

6 Impact of sample size on nonparametric methods

As discussed in §3, we run a simulation to test the effectiveness of the nonparametric method as a fitting and predictive tool. In particular, we visually assess the ability of the Bernstein polynomial model to qualitatively reproduce a known underlying distribution as a function of sample size and astrophysical scatter. We set the known distribution to be a power law in M-R space with a slope of $\sim$ 2.3 (mass as a function of radius). We simulated synthetic datasets from this known distribution and added mass and radius uncertainties of 10$\%$. We then tested this dataset for three different values of intrinsic dispersion, i.e. the scatter of the data points around the original power law. The tested values of intrinsic dispersion are 0, 0.1, and 0.5 in units of log($m$).

We note that in the simulation, the nonparametric algorithm fits the simulated dataset well in the case of low intrinsic dispersion, even for small datasets ($\sim$ 20 points; see Figure 8). However, as we increase the intrinsic dispersion, the fit compensates by increased uncertainties and visual departures from the known underlying distribution (see Figure 7).

7 Discussion

7.1 Nonparametric vs Parametric methods

In §5.3, we compare the results of the nonparametric and parametric fit to the M dwarf sample. A parametric model is useful for small data sets, where it offers an easy means to develop an intuitive understanding for the data. It can also be used when there is a known physical process that is being fit or tested for. However, in cases where we need to explore the relationship and develop a forecasting or predicting routine, nonparametric methods offer a more flexible and unbiased option. Furthermore, in cases where the underlying distributions of the variables are unknown, or the data could have small-scale structure that can not be captured by a power law, the nonparametric method is a better predictor. For this M dwarf sample, the parametric and nonparametric methods offer similar median predictions in the realm of limited sample size (see Figure 4). However, as has been shown by Ning et al. (2018) with the Kepler sample, a simple parametric power law (or broken power law) does not necessarily capture the features in the data as the sample size grows. A parametric predicting method will not necessarily reflect the data with the same accuracy and has the potential for higher bias, meaning that the predicted masses or radii will be farther from the true masses or radii. . Therefore, we propose this nonparametric method for the M-R relationship and offer the Python package MRExo as a practical tool for making these predictions.

7.2 M dwarf vs Kepler

Neil & Rogers (2018) fit an M-R relationship on a sample of M dwarf exoplanets and compared it to a sample of exoplanets orbiting FGK dwarfs to try to find evidence for any dependence on the host star mass. They did so by expanding upon the HBM approach introduced in Wolfgang et al. (2016), to include another dimension that accounts for a possible host star mass dependence. As their M dwarf sample consisted of only six planets with radii and mass measurements, they found that their M-R power law was slightly shallower than that for their FGK sample, yet still consistent with no host star mass dependence. Our sample size of 24 planets agrees with this when fit with a similar parametric relationship.

On the other hand, we find in §5.4 that the Kepler and M dwarf predictive mass distributions are different for the smallest and largest planets, in the sense that the smallest and largest planets are, on average, less massive around M dwarfs than around Sun-like stars. This result holds for both the nonparametric and the parametric M dwarf fit. The difference appears in these comparisons but not in the parametric comparison. This is because, as shown by Ning et al. (2018) the parametric model is not flexible enough to fit the larger Kepler data set, and hence does not provide a reasonable baseline for comparison.

Astrophysically, this could indicate that the low-mass protoplanetary disks around M dwarfs result in lower-mass gas giants, which would explain the lower predictive masses for the 10 $R_{\oplus}$.

For the smaller planets that are likely terrestrial, a lower mass at a given radius means a higher silicon-to-iron ratio, which may suggest that M dwarf disks process their refractory elements differently than Sun-like protoplanetary disks.

Interestingly, the predictive mass distributions for 3 $R_{\oplus}$ planets are similar between M and FGK dwarf hosts. This could indicate that the process of sub-Neptune formation is independent of the mass of the protoplanetary disk. That said, the predictive mass distributions between the M dwarf and FGK planetary samples do overlap substantially, and the M dwarf dataset is still quite small and is dominated by the TRAPPIST-1 planets. Therefore, these preliminary results should be revisited once more transiting M dwarf planets are discovered and followed up with ground-based observations.

7.3 Future Prospects

As was seen with the Kepler sample, as the sample size increases, we can start to unearth interesting phenomena that would be hard to discern from small data sets. In addition, more precise measurements with smaller error bars would help in more accurate M-R model fits. This would further unveil features and regions that were indistinguishable earlier. This will particularly hold true as the M dwarf M-R space starts to fill up with radii from TESS (Ricker et al., 2014), and the advent of high-precision RV follow-up instruments like HPF (Mahadevan et al., 2012), NEID (Schwab et al., 2016), MAROON-X (Seifahrt et al., 2016), MINERVA-RED (Blake et al., 2015), CARMENES (Quirrenbach et al., 2016), SPIROU (Thibault et al., 2012), IRD (Kotani et al., 2018), ESPRESSO (Hernández et al., 2017), NIRPS (Wildi et al., 2017), SPIRou (Artigau et al., 2014), iSHELL (Cale et al., 2018), GIANO (Claudi et al., 2017) and EXPRES (Jurgenson et al., 2016).

8 Conclusion

In this paper, we fit the M-R relationship for a sample of 24 exoplanets around M dwarfs using nonparametric and parametric methods. Considering the small sample size, the measurements are currently well described by a power law, which we fit with a parametric hierarchical Bayesian model. We then discuss the deficiencies of parametric methods and the utility of nonparametric models, especially as the sample size increases. To further illustrate this point, we also run a simulation study showing how the nonparametric fit changes with a change in sample size and dispersion in the sample. We then compare the nonparametric and parametric results, finding them to be similar, on average, but currently with a larger variance for the nonparametric mass predictions. We also discuss differences in the M-R relationship for M dwarf versus Kepler exoplanets and note that the predicted conditional probability density functions differ for the smallest and largest planets. Using this comparison, we illustrate the utility of an M dwarf M-R relationship in an era of the exciting new discoveries with TESS and ground-based precision RV instrumentation.

We also introduce a new Python package called MRExo, which can be used as a predictive tool, as well as to fit nonparametric M-R relationships to new datasets. This code is available to the community to use in its own applications.

9 Acknowledgements

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. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. Computations for this research were performed on the Pennsylvania State University’s Institute for CyberScience Advanced CyberInfrastructure (ICS-ACI). This work was partially supported by funding from the Center for Exoplanets and Habitable Worlds. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. S.K. acknowledges helpful discussions with Joe Ninan regarding the code, and Michael G. Scott and Tippy Toes for help with the project. G.K.S wishes to acknowledge support from NASA Headquarters under the NASA Earth and Space Science Fellowship Program grant NNX16AO28H. A.W. acknowledges support from the National Science Foundation Astronomy and Astrophysics Postdoctoral Fellowship program under award No. 1501440. We acknowledge support from NSF grants AST-1006676, AST-1126413, AST-1310885, and AST-1517592 in our pursuit of precision radial velocities in the near-infrared to discover and characterize planets around M dwarfs.