A Gap in the Densities of Small Planets Orbiting M Dwarfs: Rigorous Statistical Confirmation Using the Open-source Code RhoPop

We use the open-source ExoPlex ⁹ mass–radius composition calculator (Unterborn et al. 2023) to build density–mass grids from pure-water worlds to pure-iron planets. ExoPlex finds the radius for a user-input planet mass and bulk composition by self-consistently solving the five coupled differential equations: the mass within a sphere, hydrostatic equilibrium, the adiabatic temperature profile, Gauss's law of gravity in one dimension, and the thermally dependent equation of state (EOS). We assume a pure-Fe core and adopt ExoPlex's default EOS for liquid iron (Anderson & Ahrens 1994). ExoPlex couples the thermodynamic database of Stixrude & Lithgow-Bertelloni (2011) with the Perple_X thermodynamic equilibrium software (Connolly 2009) to find the equilibrium mineralogy and corresponding material density throughout the mantle. For outer water layers, ExoPlex adopts the Seafreeze Fei et al. (1993) and with the empirical equations of Asahara et al. (2010) (Journaux et al. 2020) software for liquid water, Ice Ih, II, III, V, or VI. For Ice VI, ExoPlex couples the isothermal EOS of Journaux et al. (2020).

We build isocomposition rock and liquid/solid water grids over a mass range of 0.05–10 M_⊕ in increments of ${\rm{\Delta }}{\mathrm{log}}_{10}(M/{M}_{\oplus })\simeq 0.023$. For all planets, we assume the silicate mantle is Fe-free and fix the mantle molar ratios to Si/Mg = 0.79, Al/Mg = 0.07, and Ca/Mg = 0.09 per Unterborn et al. (2023), as the Fe/Mg ratio is the primary control on density for water-free rocky planets (Dorn et al. 2015; Unterborn et al. 2016). For pure-rock compositions, we vary ${\mathrm{log}}_{10}(\mathrm{Fe}/\mathrm{Mg})$ from −0.2 to 1.5 (Fe/Mg = 0.01–30) in increments of ${\rm{\Delta }}{\mathrm{log}}_{10}(\mathrm{Fe}/\mathrm{Mg})\simeq 0.035$. Our range of Fe/Mg corresponds to a core mass fraction (CMF) range of 0.006–0.95. For water-rich compositions, we build a grid of water mass fractions (WMF) for WMF = 0.001, 0.005, 0.01, and 0.025 and then in increments of ΔWMF = 0.025 from 0.025 to 1.0. We fix the interior rocky portion of these planets to have Fe/Mg = 0.71 corresponding to the average Hypatia values per Unterborn et al. (2023). For compositions between our computed grid lines, we linearly interpolate between the two nearest compositional neighbors to approximate the intermediate isocomposition line.

We note that there is an overlap between our rock grid and water grid between WMF = 0 and 0.1–0.15, depending on mass, and Fe/Mg = 0–0.71 (CMF = 0–0.29), meaning there are two grid solutions for a given mass and radius in this regime: one with water and one without. This is a result of a well-known degeneracy when inferring planet composition from density alone (e.g., Dorn et al. 2015; Rogers 2015; Unterborn et al. 2016), which persists despite our simplifications of a pure-Fe core and Fe-free silicate mantle. We circumvent this degeneracy by using a reduced rock grid from Fe/Mg = 0.71 to 30 (CMF = 0.29–0.95). For planet densities lower than a water-free planet with Fe/Mg = 0.71, RhoPop switches to the water grid. Where our models return a WMF = 0.0–0.15, however, we also report the corresponding best-fit Fe/Mg and CMF for a water-free planet.

Following the work of Luque & Pallé (2022), we normalize all planet densities to a scaled density parameter denoted as ρ_scaled. The scaled density parameter is a function of planet mass and a reference composition. It gives the density a planet of a given mass would have with the reference composition. Both our work and Luque & Pallé (2022) assume volatile-/water-free compositions for reference. Where Luque & Pallé (2022) uses a nominally "Earth-like" CMF = 0.325, however, we assume the average CMF = 0.29 of the Hypatia Catalog corresponding to Fe/Mg = 0.71, Si/Mg = 0.79, Ca/Mg = 0.09, and Al/Mg = 0.07. That is, ρ_scaled(M) = ρ(M, CMF = 0.29). We choose this value so any planet populations inferred with RhoPop may be directly compared with the NRPZ. We can then define a general density ratio parameter ρ_ratio ≡ ρ/ρ_scaled. For the ExoPlex-derived density–mass grids, this is a function of mass and composition,

$\begin{eqnarray}&&{\rho }_{\mathrm{ratio}}=\rho (M,\mu )/{\rho }_{\mathrm{scaled}}(M),\end{eqnarray} \tag{ 1 }$

where μ is the composition in terms of a CMF from 0.29 to 1.0 or a WMF from 0 to 1.0. Our water grid and reduced rock grid are shown relative to ρ_scaled in Figure 1.

Zoom In Zoom Out Reset image size

Inspired by the population insights of Chen & Kipping (2017) and Neil et al. (2022), RhoPop employs Gaussian mixture models in a hierarchical Bayesian framework. Hierarchical Bayesian models (HBMs) employ two sets of parameters, local and hyper. We treat each observed data point as a Gaussian random variable centered on its true value with a standard deviation equal to the observational uncertainty. Since the true values cannot be known a priori, these are treated as free model parameters called "local parameters." The parameters that govern the model that predicts the local parameters are the "hyperparameters." Here the hyperparameters are those that describe the underlying populations of planets.

The local parameters of our HBM are the true planet masses, M_t , and density ratios, ρ_ratio,t . The observed planet mass, M_ob, is modeled as a random Gaussian variable centered on the true mass, M_t , with a standard deviation equal to the observed mass uncertainty, ${\sigma }_{{M}_{\mathrm{ob}}}$. That is,

$\begin{eqnarray}&&{M}_{\mathrm{ob}}^{i}\sim { \mathcal N }\left({M}_{t}^{i},{\sigma }_{{M}_{\mathrm{ob}}^{i}}\right),\end{eqnarray} \tag{ 2 }$

where the superscript i denotes the ith planet in the sample. For planets with asymmetric observational mass uncertainties, we take the average of the upper and lower uncertainties to calculate ${\sigma }_{{M}_{\mathrm{ob}}^{i}}$. Evaluating the right-hand side of Equation (2) at ${M}_{\mathrm{ob}}^{i}$ gives the likelihood function for M_t ⁱ , which we denote as ${{ \mathcal L }}_{m}^{i}$, i.e.,

$\begin{eqnarray}&&{{ \mathcal L }}_{m}^{i}={ \mathcal N }\left({M}_{\mathrm{ob}}^{i}\left.\right|{M}_{t}^{i},{\sigma }_{{M}_{\mathrm{ob}}^{i}}\right).\end{eqnarray} \tag{ 3 }$

We link the local parameters to the hyperparameters that describe the planet populations using a Gaussian mixture model. Gaussian mixture models represent some global population as a combination of N_c normally distributed subpopulations. Each subpopulation has a mixing weight, w, which is its fractional size relative to the global population. The sum of the mixing weights over all N_c components must sum to unity, i.e.,

$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{N}_{c}}{w}^{j}=1,\end{eqnarray} \tag{ 4 }$

where we use the superscript j to denote the jth subpopulation. We assume that each j component of our Gaussian mixture models is centered on composition ${\mu }_{\mathrm{comp}}^{j}$ with an intrinsic scatter of density ratios ${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}$. Within a given population, ${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}$ represents the variation that cannot be accounted for by the observational uncertainties of the planets that belong to it. We then model the local parameter ${\rho }_{\mathrm{ratio},t}^{i}$ as a Gaussian mixture of all N_c components,

$\begin{eqnarray}&&{\rho }_{\mathrm{ratio},t}^{i}\sim \displaystyle \sum _{j=1}^{{N}_{c}}{w}^{j}\times { \mathcal N }\left(f({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j}),{\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}\right).\end{eqnarray} \tag{ 5 }$

The argument of the Gaussian in Equation (5), $f({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j})$, is calculated directly from the compositional grids (Section 2.1) and represents the predicted density ratio at M_t ⁱ given a central composition of population j, ${\mu }_{\mathrm{comp}}^{j}$. We then model the observed density ratio of the ith planet, ${\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}$, as being centered on ${\rho }_{\mathrm{ratio},t}^{i}$ with an error term corresponding to its observed uncertainty,

$\begin{eqnarray}&&\begin{array}{l}{\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}\sim \displaystyle \sum _{j=1}^{{N}_{c}}\left[{w}^{j}\times { \mathcal N }\left(f\left({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j}\right),{\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}\right)\right.\\ \qquad \left.+{ \mathcal N }\left(0,{\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}}\right)\right]\\ \qquad =\displaystyle \sum _{j=1}^{{N}_{c}}{w}^{j}\times { \mathcal N }\left(f\left({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j}\right),\right.\\ \qquad \left.\sqrt{{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}\right)}^{2}+{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}}\right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ 6 }$

As with Equation (2), evaluating the right-hand side of Equation (6) at ${\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}$ gives the likelihood function for ${\rho }_{\mathrm{ratio},t}^{i}$,

$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal L }}_{\rho }^{i}=\displaystyle \sum _{j=1}^{{N}_{c}}{w}^{j}\times { \mathcal N }\left({\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}\left.\right|f\left({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j}\right),\right.\\ \qquad \left.\times \,\sqrt{{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}\right)}^{2}+{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}}\right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ 7 }$

Finally, combining Equations (3) and (7), taking the natural logarithm, and summing over all N_p planets gives the log-likelihood function for the entire model as a function of the hyperparameters,

$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{i=1}^{{N}_{p}}\left[\mathrm{ln}{{ \mathcal L }}_{\rho }^{i}+\mathrm{ln}{{ \mathcal L }}_{m}^{i}\right]\\ \qquad =\displaystyle \sum _{i=1}^{{N}_{p}}\left[\mathrm{ln}\left(\displaystyle \sum _{j=1}^{{N}_{c}}{w}^{j}\times { \mathcal N }\left({\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}\left.\left.\right|f\left({M}_{t}^{i},{\mu }_{\mathrm{comp}}^{j}\right),\right.\right.\right.\\ \qquad \left.\left.\left.\times \,\sqrt{{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}\right)}^{2}+{\left({\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}}\right)}^{2}}\right)\right)\right.\\ \qquad \left.+\mathrm{ln}\left({ \mathcal N }\left({M}_{\mathrm{ob}}^{i}\right|{M}_{t}^{i},{\sigma }_{{M}_{\mathrm{ob}}^{i}}\right)\right)\right].\end{array}\end{eqnarray} \tag{ 8 }$

Our HBM framework is summarized schematically in Figure 2, and all parameters are summarized in Table 1. In this work, we compare models through Bayesian model evidence ${ \mathcal Z }$: the integral of the likelihood (Equation (8)) over the entire parameter prior space. Mathematically,

$\begin{eqnarray}&&P({\boldsymbol{D}}| M)\equiv { \mathcal Z }={\int }_{{{\rm{\Omega }}}_{{\boldsymbol{\Theta }}}}{ \mathcal L }({\boldsymbol{D}}| {\boldsymbol{\Theta }},{ \mathcal M }){\Pr }({\boldsymbol{\Theta }}| { \mathcal M })d{\boldsymbol{\Theta }},\end{eqnarray} \tag{ 9 }$

where Θ represents a generic set of local parameters and hyperparameters, Pr are the corresponding priors, and ${ \mathcal M }$ is the model. We choose wide uninformative uniform priors on all hyperparameters to remain agnostic when validating previous results in Sections 3.1 and 3.2. We summarize all hyperparameter priors in Table 2. RhoPop handles up to three populations, i.e., N_c = 1, 2, or 3. RhoPop makes no a priori assumptions about the central compositions of the populations ${\mu }_{\mathrm{comp}}^{j}$ other than ${\mu }_{\mathrm{comp}}^{1}\geqslant {\mu }_{\mathrm{comp}}^{2}$ for N_c = 2 and ${\mu }_{\mathrm{comp}}^{1}\geqslant {\mu }_{\mathrm{comp}}^{2}\geqslant {\mu }_{\mathrm{comp}}^{3}$ for N_c = 3. In the three-population scenario, we allow w¹ and w² to vary from 0 to 1. However, if the sum of these two is greater than 1, the likelihood function returns a $-\inf$, and that set of parameters is rejected since Equation (4) must be satisfied. We assume priors of ${M}_{t}^{i}\sim { \mathcal N }({M}_{\mathrm{ob}}^{i},{\sigma }_{{M}_{\mathrm{ob}}^{i}})$ for all i in N_p .

Zoom In Zoom Out Reset image size

Table 1. Summary of Parameter Definitions

Parameter	Parameter Type	Description
${\mu }_{\mathrm{comp}}^{j}$	Hyper	Central composition of population j in WMF or CMF
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{j}}$	Hyper	Intrinsic 1σ density ratio scatter of population j
wj	Hyper	Mixing weight of population j
${\rho }_{\mathrm{ratio},t}^{i}$	Local	True density ratio of planet i
Mt i	Local	True mass of planet i
${M}_{\mathrm{ob}}^{i}$	Data	Central observed mass of planet i
${\sigma }_{{M}_{\mathrm{ob}}^{i}}$	Data	1σ uncertainty in ${M}_{\mathrm{ob}}^{i}$
${\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}$	Data	Central observed density ratio of planet i
$\sigma {\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}$	Data	1σ uncertainty in ${\rho }_{\mathrm{ratio},\mathrm{ob}}^{i}$

Download table as:  ASCIITypeset image

Table 2. Hyperparameter Priors

Parameter	Nc = 1	Nc = 2	Nc = 3
${\mu }_{\mathrm{comp}}^{1}$	${ \mathcal U }$[1.0 WMF, 0.95 CMF]	${ \mathcal U }$[1.0 WMF, 0.95 CMF]	${ \mathcal U }$[1.0 WMF, 0.95 CMF]
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{1}}$	${ \mathcal U }$[0, 2]	${ \mathcal U }$[0, 2]	${ \mathcal U }$[0, 2]
w1	1	${ \mathcal U }$[0, 1]	${ \mathcal U }$[0, 1]
${\mu }_{\mathrm{comp}}^{2}$	⋯	${ \mathcal U }[0,1]\times {\mu }_{\mathrm{comp}}^{1}$	${ \mathcal U }[0,1]\times {\mu }_{\mathrm{comp}}^{1}$
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{2}}$	⋯	${ \mathcal U }$[0, 2]	${ \mathcal U }$[0, 2]
w2	⋯	1 − w1	${ \mathcal U }$[0, 1]
${\mu }_{\mathrm{comp}}^{3}$	⋯	⋯	${ \mathcal U }[0,1]\times {\mu }_{\mathrm{comp}}^{2}$
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{3}}$	⋯	⋯	${ \mathcal U }$[0, 2]
w3	⋯	⋯	1 − w1 − w2

Note. The default prior on ${\mu }_{\mathrm{comp}}^{1}$ spans the entirety of our water and reduced rock grids.

Download table as:  ASCIITypeset image

In practice, calculating ${ \mathcal Z }$ analytically is not tractable, and it must be computed numerically. As such, RhoPop uses the dynesty ¹⁰ dynamic nested sampling software package with the default weighting scheme, which is optimized to simultaneously estimate the parameter posteriors and model evidence. We point the reader to Speagle (2020) for more details.

2.2.1. Model Selection

The full sample of Luque & Pallé (2022) contains 34 planets with mass and radius precisions of better than 25% and 8%, respectively, and radii of less than 4 R_⊕. The authors identify seven planets as mini-Neptunes. We do not consider these planets further as they require an extended H/He envelope to explain and do not meet our definition of a small planet. We instead focus on the remaining 27 planets from which Luque & Pallé (2022) observe a bimodal density distribution: a higher-density population with 21 members and a lower-density population with six members. The authors interpret the higher-density population as being purely rocky in composition and the lower-density population as having a rocky interior surrounded by a 50% H₂O layer by mass.

We run this 27-planet sample through RhoPop with N_c = 1 and 2, i.e., one- and two-population models. We visualize our results in Figure 3 and summarize the hyperparameter posterior statistics in Table 3. We find that the two-population model is strongly preferred with a Bayes factor B₂₁ = 1042. In frequentist parlance, this corresponds to a p-value of 3.4 × 10⁻⁵ or a 4.1σ detection of two populations. For good measure, we run this sample using the three-population model. We find B₃₂ = 0.26, meaning the two-population model is indeed favored over the three-population model.

Zoom In Zoom Out Reset image size

Table 3. Summary of Parameter Posteriors, Model Evidence, and Model Selection Parameters for Our One- and Two-population Model Runs with the Luque & Pallé (2022) Sample

	Nc = 1		Nc = 2
Param.	Med.	68% CI	Med.	68% CI
${\mu }_{\mathrm{comp}}^{1}$	0.09 WMF	[0.14 WMF, 0.06 WMF]	0.02 WMF	[0.03 WMF, 0.017 WMF]
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{1}}$	0.21	[0.17, 0.25]	0.02	[0.005, 0.04]
w1	1	⋯	0.76	[0.67, 0.83]
${\mu }_{\mathrm{comp}}^{2}$	⋯	⋯	0.78 WMF	[0.87 WMF, 0.69 WMF]
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{2}}$	⋯	⋯	0.03	[0.008, 0.07]
w2	⋯	⋯	1-w1	⋯
	⋯	$\mathrm{ln}{{ \mathcal Z }}_{1}=7.20\pm 0.24$	⋯	$\mathrm{ln}{{ \mathcal Z }}_{2}=14.15\pm 0.32$

B21	p-value	Sigma	Interpretation
1042	3.4 × 10−5	4.1σ	Strong evidence for two populations

Note. CI = credible interval.

Download table as:  ASCIITypeset image

As further validation, our two-population model identifies the same 21 planets as Luque & Pallé (2022) belonging to a higher-density population. We find that this population is best fit by a central WMF of 0.02. This WMF falls within the degenerate zone described in Section 1. As such, we apply a χ² minimization with these 21 planets over our full rock grid to find the equivalent CMF assuming the planets are water-free. We find a CMF = 0.28 (Fe/Mg = 0.67), meaning this population is either slightly wet or has a ∼6% Fe depletion relative to the average volatile-/water-free NRPZ composition of CMF = 0.29 (Fe/Mg = 0.71).

We find the same six planets that Luque & Pallé (2022) identify as being 50% water worlds as indeed belonging to a separate lower-density population. We, however, find that this population is best fit by a central WMF of 0.78. While this WMF value is substantially more than that predicted by Luque & Pallé (2022), this is primarily due to the differences in the underlying material EOSs used in this work versus Luque & Pallé (2022). Both populations are found to have very little intrinsic scatter. Even at 10 M_⊕, where the difference in density ratios is minimized over the mass range considered, the two populations are separated by more than 12× their mutual intrinsic scatter parameters, meaning there is a robust density gap between these two small-planet populations.

The Adibekyan et al. (2021) sample includes 22 likely volatile-free rocky planets with M_p < 10 M_⊕ and mass and radius uncertainties of <30% around FGK stars. The authors identify a visually apparent but not statistically significant subpopulation of five iron-rich planets. We run our one- and two-population models with this sample and find B₂₁ = 0.39, meaning there is no support for the more complex two-population over the one-population model. Thus, the single-population model shown in Figure 4 is our preferred model for this data set. Given the lack of support for the two-population model for this data set, we do not consider the three-population scenario. The results for the two-population and preferred one-population model are summarized in Table 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Summary of Parameter Posteriors, Model Evidence, and Model Selection Parameters for Our One- and Two-population Model Runs with the Adibekyan et al. (2021) Sample

	Nc = 1		Nc = 2
Param.	Med.	68% CI	Med.	68% CI
${\mu }_{\mathrm{comp}}^{1}$	0.31 CMF	[0.02 WMF, 0.36 CMF]	0.36 CMF	[0.30 CMF, 0.53 CMF]
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{1}}$	0.16	[0.11, 0.22]	0.17	[0.06, 0.69]
w1	1	⋯	0.64	[0.07, 0.94]
${\mu }_{\mathrm{comp}}^{2}$	⋯	⋯	0.05 WMF	[0.53 WMF, 0.33 CMF]
${\sigma }_{{\rho }_{\mathrm{ratio},\mathrm{pop}}^{2}}$	⋯	⋯	0.19	[0.09, 0.77]
w2	⋯	⋯	1-w1	⋯
	⋯	$\mathrm{ln}{{ \mathcal Z }}_{1}=-25.991\pm 0.223$	⋯	$\mathrm{ln}{{ \mathcal Z }}_{2}=-26.93\pm 0.21$

B21	p-value	Sigma	Interpretation
0.39	1.0	0	No evidence for more than one population

Note. CI = credible interval.

Download table as:  ASCIITypeset image

We compare our results for the Luque & Pallé (2022; Section 3.1) and Adibekyan et al. (2021; Section 3.2) samples with the NRPZ (Unterborn et al. 2023). From Equations (24) and (27) of Unterborn et al. (2023), we estimate that the NRPZ covers a density ratio range of 0.85–1.1. We reiterate that the NRPZ represents the expected 3σ compositional range of volatile-free and water-free rocky planets. The preferred single-population model for the Adibekyan et al. (2021) sample yields a 1σ density ratio range that is similar in width to the NRPZ, meaning the 3σ range for this sample is ∼3× larger than the NRPZ. This result is not surprising given the number of exceptions to this simplistic picture of planet formation we outline in Section 1.

Our one-population model for the Luque & Pallé (2022) sample finds a central composition of 0.09 WMF with an intrinsic 1σ scatter of ${\sigma }_{{\rho }_{\mathrm{pop}}}=0.21$. This leads to a 1σ density range of ∼0.6–1.0, which is inconsistent with the NRPZ. This result is not surprising since we expected to find two populations a priori given the results of Luque & Pallé (2022), although we did not inject this expectation into our models. This tension with the NRPZ is reduced for the strongly preferred two-population model. We find a higher-density population with a central composition of WMF = 0.02 and ${\sigma }_{{\rho }_{\mathrm{pop}}}=0.02$. The 2% water mass fraction isocomposition line corresponds to a density ratio of ∼0.95, and the 3σ range in density ratios is ∼0.89–1.01. This density ratio range falls entirely within the NRPZ, meaning that this population of planets can be explained via volatile-free rocky compositions; i.e., we do not need to invoke water to explain the densities of these planets. More notable, however, is the fact that the 3σ range of this population is approximately half the width of the NRPZ, suggesting that the compositions of these planets may be less diverse than stellar abundances predict. We note the important caveat that the NRPZ is estimated using Fe, Mg, and Si abundances from the Hypatia Catalog database (Hinkel et al. 2014) and is thus strongly biased toward FGK stars, as there is a lack of M-type stars with abundance measurements, especially those with Fe, Mg, and Si measurements. This suggests that the Luque & Pallé (2022) sample has a smaller spread in the major rock-building elements than FGK stars. In fact, M dwarfs that host planets in general may have a smaller spread. Testing this possibility, however, requires more M dwarf abundances of Fe, Mg, and Si to build an M dwarf NRPZ, i.e., an M-NRPZ. While TOI-561 b in the Adibekyan et al. (2021) sample has a ρ_ratio ∼ 0.5, consistent with the lower-density population of Luque & Pallé (2022), the lack of evidence for a water-rich population around FGK stars further suggests that small-planet compositions and formation processes may be dependent on host-star type.

The primary goal of Section 3.1 is to validate the results of Luque & Pallé (2022) using their sample as published. We note, however, that the TRAPPIST-1 system accounts for roughly a quarter of the LP22 sample. Further, the seven TRAPPIST-1 planets have some of the smallest error bars with an average density uncertainty of 5.6% compared to the remainder of the sample, which has an average density uncertainty of 19%. Together, this suggests that our detection of two populations of M dwarf planets may be driven in large part by the TRAPPIST-1 system. To test this, we removed the TRAPPIST-1 planets from the LP22 sample and reran the one- and two-population models. Without these seven planets, the detection of two populations of planets is reduced from 4.1σ to 1.8σ, meaning the TRAPPIST-1 planets do indeed play a large role in the detection of two populations of small planets around M dwarfs. Further, for the N_c = 2 scenario, we find a median value of ${\mu }_{\mathrm{comp}}^{1}=0.327$ CMF, nearly identical to the CMF of Earth, in contrast to ${\mu }_{\mathrm{comp}}^{1}=0.02$ WMF, when the TRAPPIST-1 planets are included. This begs the question as to whether volatile-/water-free rocky planets around M dwarfs are generally Earth-like in composition and the TRAPPIST-1 planets are anomalously iron-depleted or wet. A larger sample of small M dwarf planets with density precisions comparable to the TRAPPIST-1 planets and/or an M-NRPZ are needed to fully investigate this question.

We note that our water grid assumes condensed water/ice phases. The equilibrium temperatures of the five LP22 planets in the lower-density population range from ∼1.5× to 4.2×Earth's. As such, it is likely that substantial fractions of their water exist in the vapor and supercritical phases. Our choice of condensed water/ice then overestimates the density of these water layers; therefore, our WMF estimates serve as an upper limit. We chose to work in terms of condensed water/ice for the most direct comparison with the interpretations made in Luque & Pallé (2022), but a future iteration of RhoPop will incorporate the supercritical and vapor phases of water and planet equilibrium temperature.

We also note that the population of planets interpreted as water worlds in Luque & Pallé (2022) and this work is not uniquely explained as such. Rogers et al. (2023) show that these planets are equally well explained as sub-Neptunes that have retained a small fraction of their primordial H/He envelopes and provide some observational evidence to support this interpretation. It is also plausible that such planets have outer gaseous envelopes comprised of species heavier than H/He (e.g., Piaulet et al. 2023). In short, the source M dwarf planet density gap remains to be fully explored. In a future update, we will expand RhoPop to account for these various scenarios and investigate whether the data favor one scenario over the others.

In Section 4, we estimated that approximately 154 planets with average ${\sigma }_{{M}_{p}}=10 \%$ and ${\sigma }_{{R}_{p}}=2.5 \%$ are needed to resolve a population of planets with Mercury-like compositions from a population with nominally Earth-like compositions. Interestingly, there are 154 small planets in the NASA Exoplanet Archive ¹¹ with density uncertainties of less than 100%. The average and median M–R uncertainties of these planets, however, are ${\sigma }_{{M}_{p}}\,=27 \%$ and ${\sigma }_{{R}_{p}}=6 \%$ and ${\sigma }_{{M}_{p}}=19 \%$ and ${\sigma }_{{R}_{p}}=5 \%$, respectively. These uncertainties are ≳2× the precisions required to resolve an Fe-rich population with high confidence for N_p = 154. As such, even if such a population exists, we do not expect it to be resolvable with the current sample of small planets. We note, however, that our analysis assumes there is one Fe-rich planet for every three nominally Earth-like planets (w^ml = 0.25) based on the work of Adibekyan et al. (2021). If Fe-rich planets are more abundant relative to Earth-like planets than 1:3, then the required sample size to resolve this population will be reduced, and vice versa. Further, we assume a constant bulk composition for the iron-rich population of ${\mu }_{\mathrm{comp}}^{\mathrm{ml}}=0.7$. While we use Mercury and the five potentially Fe-rich planets from Adibekyan et al. (2021) to justify this choice, the average CMF can be higher or lower, which will act to make this population more or less easily resolvable at a given mass and radius precision, respectively. We will consider a wider range of mixing weights and compositions for the Fe-rich population in future work.