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

Using mass–radius composition models, small planets (R ≲ 2 R ⊕) are typically classified into three types: iron-rich, nominally Earth-like, and those with solid/liquid water and/or atmosphere. These classes are generally expected to be variations within a compositional continuum. Recently, however, Luque & Pallé observed that potentially Earth-like planets around M dwarfs are separated from a lower-density population by a density gap. Meanwhile, the results of Adibekyan et al. hint that iron-rich planets around FGK stars are also a distinct population. It therefore remains unclear whether small planets represent a continuum or multiple distinct populations. Differentiating the nature of these populations will help constrain potential formation mechanisms. We present the RhoPop software for identifying small-planet populations. RhoPop employs mixture models in a hierarchical framework and a nested sampler for parameter and evidence estimates. Using RhoPop, we confirm the two populations of Luque & Pallé with >4σ significance. The intrinsic scatter in the Earth-like subpopulation is roughly half that expected based on stellar abundance variations in local FGK stars, perhaps implying M dwarfs have a smaller spread in the major rock-building elements (Fe, Mg, Si) than FGK stars. We apply RhoPop to the Adibekyan et al. sample and find no evidence of more than one population. We estimate the sample size required to resolve a population of planets with Mercury-like compositions from those with Earth-like compositions for various mass–radius precisions. Only 16 planets are needed when σMp=5% and σRp=1% . At σMp=10% and σRp=2.5% , however, over 154 planets are needed, an order of magnitude increase.


INTRODUCTION
Small planets (R ≲ 2R ⊕ ) are among the most common in the Galaxy.Standard formation models predict that the relative amounts of the major rock-building elements (Fe/Mg and Si/Mg) in small planets should reflect that of their host stars (e.g., Bond et al. 2010a,b;Elser et al. 2012;Johnson et al. 2012;Thiabaud et al. 2015).The inferred compositions of observed small exoplanets (Schulze et al. 2021;Adibekyan et al. 2021;Unterborn et al. 2023) and rocky exoplanetary material (Bonsor et al. 2021) are consistent with this host-planet compositional connection at current observational precisions except in extreme cases, e.g., Kepler-107 c (Bonomo et al. 2019) and HD 137496 b (Azevedo Silva et al. 2022).In the Solar System, Earth and Mars have compositions that are consistent with the solar abundances of the major rock-building elements (e.g., Wanke & Dreibus 1994;Lodders 2003; Wang et al. 2019;Yoshizaki & McDonough 2020).The composition of Venus is unknown but expected to be consistent with Earth (e.g., Zharkov 1983;Aitta 2012;Dumoulin et al. 2017).Mercury, in contrast, has a 200-400% overabundance of Fe relative to the Sun, indicative of a fundamentally distinct formation regime and/or subsequent evolution from the other small planets in the Solar System (e.g., Cameron 1985;Benz et al. 1988;Wurm et al. 2013;Johansen & Dorn 2022).If we, therefore, assume that small planets are barren and rocky, then inconsistencies between inferred planet composition and host composition indicate deviations from standard planet formation models/theories.
Recently, Luque & Pallé (2022) observed a bimodal density distribution for 27 small planets in orbit around M dwarf hosts.The authors interpret this bimodality as two compositionally distinct subpopulations, one with Earthlike compositions and one with 50% water and 50% rock by mass, but do not provide robust statistical evidence to support this claim or a comparison with host abundances of the major rock-building elements.There is a paucity of known M dwarf abundances of the major rock-building elements as measuring M dwarf abundances is a non-trivial problem and an active area of research (e.g., Allard et al. 2000;Souto et al. 2022).As such, it is unclear whether there are indeed two subpopulations of planets around M dwarfs and to what degree these planets deviate from the host-planet compositional connection.Adibekyan et al. (2021) show that 22 likely rocky small planets in orbit around FGK stars have compositions correlated with their hosts' but not with a 1:1 relationship.The authors identify a class of 5 planets with a Mercury-like overabundance of iron separated from the general star-planet compositional link by a visually apparent, but not statistically conclusive, gap in density.Barros et al. (2022) identify two additional planets around HD 23472 that appear to belong to the same iron-enriched population identified in Adibekyan et al. (2021).Together these studies suggest that the small-planet density distribution may not be a continuum.Instead, planet formation mechanisms might create up to three distinct classes of small planets: (1) those that follow the star-planet compositional link, (2) those enriched in water/atmosphere and/or depleted in iron, and (3) those enriched in iron and/or depleted in silicates.
Statistical host-planet composition comparisons like that of Adibekyan et al. (2021) provide a powerful tool for confirming iron-enriched/silicate-depleted or water-enriched/iron-depleted small planets (e.g., Santerne et al. 2018;Schulze et al. 2021;Rodríguez Martínez et al. 2022;Azevedo Silva et al. 2022;Unterborn et al. 2023).There is, however, a paucity of planets with known densities and host abundances of Fe, Mg, and Si, and no single catalog lists all of these quantities simultaneously.Combining the NASA Exoplanet Archive1 with the Hypatia Catalog stellar abundance database2 (Hinkel et al. 2014), we estimate there are approximately 43 planets with measured densities and host Fe, Mg, and Si abundances.Of these, only 12 are small planets.Direct host-planet comparisons are therefore limited to single planets or small samples (e.g., 11, 22, and 7 in Schulze et al. 2021;Adibekyan et al. 2021;Unterborn et al. 2023, respectively).
To overcome this lack of host-star abundances, some studies instead compare inferred small-planet compositions with a large sample of stellar abundances of the major rock-building elements (Scora et al. 2020;Plotnykov & Valencia 2020;Unterborn et al. 2023).Unterborn et al. (2023) used the Hypatia Catalog stellar abundance database (Hinkel et al. 2014) to quantify a nominally rocky planet zone (NRPZ): the expected 99.7% range of rocky planet densities owing to the intrinsic spread in stellar compositions.Where the density and corresponding uncertainty of a small planet exclude it from the NRPZ, its bulk composition likely deviates from that predicted through standard formation models.Excluded planets whose densities exceed that of the NRPZ are then interpreted as having an iron enrichment and/or silicate depletion.Excluded planets with densities less than the NRPZ are interpreted as being iron-depleted, water-rich, or volatile-rich.Applying this technique to a large sample of small planets, Unterborn et al. (2023) find 21 planets that are excluded from the NRPZ with ≥ 90% confidence: 8 are likely iron-enriched and 13 are iron-depleted, water-rich, or volatile-rich.
It remains an open question whether the densities of small planets reflect a continuous variation in composition or multiple, distinct subpopulations.Robust identification of these subpopulations, or lack thereof, is an imperative first step to better constraining their formation and subsequent evolution.Direct host-planet compositional comparisons are a powerful method for identifying planets that deviate from standard formation models but are sample-size-limited, making it difficult to resolve distinct subpopulations of small planets.This limitation is often circumvented by instead comparing planet compositions with a large sample of stellar abundances.While the latter has been used to identify a number of individual planets that are inconsistent with the NRPZ, no study has investigated whether these planets form a continuum with the NRPZ or are members of fundamentally different subpopulations.
In this work, we present the open-source RhoPop3 code for robust identification of compositionally distinct populations of small planets.RhoPop models the density distribution of an input planet sample as a Gaussian mixture of up to three compositionally-distinct subpopulations in a hierarchical Bayesian framework.Similar models have been used for exoplanet population analyses, but small planets are treated as a single population (Chen & Kipping 2017) or assume a fixed composition for rocky materials (Neil & Rogers 2020;Neil et al. 2022).RhoPop uses compositional grids that span from pure-iron to pure-water planets that are normalized relative to the average composition of the NRPZ, so that model results are easily comparable with the expected compositional range of small water-/volatile-free planets.We describe the inner workings of our software in Sect. 2. We use RhoPop to validate the results of Luque & Pallé (2022) and Adibekyan et al. (2021) in Sect.3.1 and 3.2, respectively.In Sect.4, we use RhoPop to estimate the minimum number of planets needed to identify a population of Mercury-like planets from a population of Earth-like planets for various mass and radius uncertainties.

Density-mass-composition grids
We use the open-source ExoPlex4 mass-radius-composition calculator (Unterborn et al. 2023) to build densitymass 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, adiabatic temperature profile, Gauss's law of gravity in one dimension, and the thermally dependent equation of state.We assume a pure-Fe core and adopt ExoPlex's default equation of state (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 (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) with the empirical equations of Asahara et al. (2010) and Fei et al. (1993).
We build isocomposition rock and liquid/solid water grids over a mass range of 0.05 to 10 M ⊕ in increments of ∆ log 10 (M/M ⊕ ) ≃ 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 log 10 (Fe/Mg) from -0.2 to 1.5 (Fe/Mg = 0.01-30) in increments of ∆ log 10 (Fe/Mg) ≃ 0.035.Our range of Fe/Mg corresponds to a core mass fraction (CMF) range of 0.006 to 0.95.For water-rich compositions, we build a grid of WMF = 0.001, 0.005, 0.01, 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., Rogers 2015;Dorn et al. 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-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 nominally rocky planet zone.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 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 Fig. 1.

Hierarchical Bayesian Model Framework
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 (HBM) 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 can not 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 .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, σ M ob .That is, where the superscript i denotes the i-th planet in the sample.For planets with asymmetric observational mass uncertainties, we take the average of the upper and lower uncertainties to calculate σ M i ob .Evaluating the right-hand side of Eq. 2 at M i ob gives the likelihood function for M i t , which we denote as L i m , i.e., We link the local parameters to the hyperparameters that describe the planet populations using a Gaussian mixture model (GMM).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., Nc j=1 where we use the superscript j to denote the j-th subpopulation.We assume that each j component of our Gaussian mixture models is centered on composition µ j comp with an intrinsic scatter of density ratios σ ρ j ratio,pop . Within a given population, σ ρ j ratio,pop 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 ρ i ratio,t as a Gaussian mixture of all N c components . (5) The argument of the Gaussian in Eq. 5, f (M i t , µ j comp ), is calculated directly from the compositional grids (Sect.2.1) and represents the predicted density ratio at M i t given a central composition of population j, µ j comp .We then model the observed density ratio of the i-th planet, ρ i ratio,ob , as being centered on ρ i ratio,t with an error term corresponding to its observed uncertainty, As with Eq. 2, evaluating the right-hand side of Eq. 6 at ρ i ratio,ob gives the likelihood function for ρ i ratio,t , Finally, combing Equations 3 and 7, taking the natural logarithm, and summing over all N p planets gives the loglikelihood function for the entire model as a function of the hyperparameters, 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 Z: the integral of the likelihood (Eq.8) over the entire parameter prior space.Mathematically, where Θ represents a generic set of local parameters and hyperparameters, P r are the corresponding priors, and 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 µ j comp other than µ 1 comp ≥ µ 2 comp for N c = 2 and µ 1 comp ≥ µ 2 comp ≥ µ 3 comp for N c = 3.In the 3-population scenario, we allow w 1 and w 2 to vary from 0 to 1.However, if the sum of these 2 is greater than one, the likelihood function returns a − inf, and that set of parameters is rejected since Eq. 4 must be satisfied.We assume priors of In practice, calculating Z analytically is not tractable and it must be computed numerically.As such, RhoPop uses the dynesty5 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.

Model selection
We adopt the methods outlined in Benneke & Seager (2013) which use Bayes factors to compare models of increasing complexity.A general Bayes factor B ij is defined as the ratio of the Bayesian evidence of model i to the Bayesian evidence of model j, i.e., B ij ≡ Z i /Z j .In terms of hypothesis testing, the denominator represents the evidence for the null hypothesis and the numerator is the alternative hypothesis.A value of B ij < 1 suggests the data favor model j over model i i.e., the data are consistent with the null hypothesis.A B ij from 1-3 corresponds to inconclusive support for the more complex model i, and the null hypothesis cannot be rejected.A B ij from 3-12, 12-150, and > 150 correspond to weak, moderate, and strong support, respectively, for the alternative hypothesis.For a given sample in this work, we first run the 1-and 2-population models and compare via B 21 = Z 2 /Z 1 where Z 1 and Z 2 correspond to the evidence for the 1-and 2-population models, respectively.Where there is significant support for the 2-population model (B 21 > 150), we then run the 3-population model and compare it with the 2-population model via B 32 .Benneke & Seager (2013), and references therein, provide formulae to convert Bayes factors to the frequentist p-value (see their Eq.10) and corresponding sigma significance n σ (see their Eq.11).For all Bayes factors in this work, we use these formulae to calculate and report the equivalent p-values and n σ significance levels for the frequentist reader.
Table 2. Hyperparameter priors.The default prior on µ 1 comp spans the entirety of our water and reduced rock grids.

Luque & Pallé (2022) sample
The full sample of Luque & Pallé (2022) contains 34 planets with mass and radius precisions of better than 25% and 8%, respectively, and radii less than 4 R ⊕ .The authors identify 7 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 6 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 2 O layer by mass.
We run this 27-planet sample through RhoPop with N c = 1 and 2, i.e., 1-and 2-population models.We visualize our results in Figure 3 and summarize the hyperparameter posterior statistics in Table 3.We find the 2-population model is strongly preferred with a Bayes factor B 21 = 1042.In frequentist parlance, this corresponds to a p-value of 3.4 × 10 −5 or a 4.1σ detection of 2-populations.For good measure, we run this sample using the 3-population model.We find B 32 = 0.26 meaning the 2-population model is indeed favored over the 3-population model.
As further validation, our 2-population model identifies the same 21 planets as Luque & Pallé (2022) belonging to a higher density population.We find this population is best fit by a central WMF of 0.02.This WMF falls within the degenerate zone described in Sect. 1.As such, we apply a χ 2 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 6 planets that Luque & Pallé (2022) identify as being 50% water worlds as indeed belonging to a separate lower-density population.We, however, find 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 equations of state used in this work versus Luque & Pallé (2022).Both populations are found to have very little intrinsic scatter.Even at 10M ⊕ , where the difference in density ratios is minimized over the mass range considered, the two populations are separated by more than 12x their mutual intrinsic scatter parameters, meaning there is a robust density gap between these two small-planet populations.

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

MINIMUM NUMBER OF PLANETS REQUIRED TO RESOLVE A POPULATION OF IRON-RICH PLANETS
Here we assume the 5 potentially iron-rich planets identified by Adibekyan et al. (2021) do indeed represent a distinct higher-density population from a lower-density population comprised of the remaining 17 planets and that our null  2022) sample of small transiting planets in orbit around M-dwarfs.We color code the highest-density population as green and the lowest-density population as purple.The planets are assigned the same color as the population to which they most likely belong.The solid lines are the isocomposition lines for each µ j comp .The bands are the corresponding intrinsic density ratio scatter parameter σ j ρ ratio,pop .For all µ j comp and σ j ρ ratio,pop we use the median values given in Table 3.  .
finding in Sect.3.2 is a result of observational uncertainty alone.We then estimate the minimum number of planets needed to resolve a population of iron-enriched planets from a lower-density population of volatile-/water-free rocky planets.Evidencing the Solar System scaled-density dichotomy between Mercury and the other rocky planets, we assume the lower-density population has an Earth-like composition, i.e., µ el comp = 0.325 CMF, and the higher density population has a Mercury-like core mass fraction of µ ml comp = 0.7 CMF.The latter choice is further supported by the fact that χ 2 is minimized between the 5 Fe-rich planets of Adibekyan et al. (2021) and our full rock grid when CMF∼ 0.73.For a sense of scale, the Hypatia catalog has Fe/Mg = 0.79 ± 0.18 corresponding to CMF = 0.29 ± 0.05, meaning  2021) sample.We color code the highest-density population as green and the lowest-density population as purple.The planets are assigned the same color as the population to which they most likely belong.The solid lines are the isocomposition lines for each µ j comp .The bands are the corresponding intrinsic density ratio scatter parameter σ j ρ ratio,pop .For all µ j comp and σ j ρ ratio,pop we use the median values given in Table 4.
this population is ∼ 8σ above the central CMF of Hypatia.Per the Adibekyan et al. ( 2021) sample, Fe-rich planets represent 0.23 of the total population.For simplicity, we round this up to a quarter and set the mixing weight of our Mercury-like population to w ml = 0.25.The corresponding mixing weight for the Earth-like population is then w el = 0.75.
For each sample, we randomly draw N p planet masses from M p ∼ U[1, 10]M ⊕ .We then assume N ml = w ml × N p of these are Mercury-like and N el = w el × N p are Earth-like in composition.We then calculate the density corresponding to µ ml comp = 0.7 CMF and µ el comp = 0.325 CMF at mass for each Mercury-like and Earth-like planet, respectively, and then calculate the corresponding planet radius, R p .Finally, we resample each mass and radius from the Gaussian distributions ∼ N (M p , σ Mp ) and ∼ N (R p , σ Rp ), respectively, to simulate observational uncertainties.We consider three sets of mass-radius uncertainties: (1) σ Mp = 10% and σ Rp = 2.5%, (2) σ Mp = 6% and σ Rp = 2%, and (3) σ Mp = 5% and σ Rp = 1%.For each set of mass-radius uncertainties, we calculate B 21 and n σ as a function of N p .
We find that the number of planets required to achieve a strong detection of both populations is highly dependent on the average mass and radius uncertainties.For the near best-case scenario of σ Mp = 5% and σ Rp = 1%, we estimate only ∼ 12 planets are needed to detect both populations with moderate support (B 21 > 12 and n σ > 2.7) and ∼ 16 for detection with strong support (B 21 > 150 and n σ > 3.6).For the most conservative, but still precise, case of σ Mp = 10% and σ Rp = 2.5%, a moderate detection requires approximately 117-118 planets, and a strong detection requires at least ∼ 154 planets, i.e., an almost factor of 10 increase in N p from the near best-case scenario.For comparison, the average mass and radius uncertainties for the Adibekyan et al. (2021) sample are 14% and 5%, respectively.As such, it is not surprising that they do not find a conclusive and distinct population of Fe-rich planets as it would require significantly more than 154 planets to do so at their average precisions.We also show our model results for the synthetic 20-planet sample with σ Mp = 5% and σ Rp = 1% in detail in Appendix 8 to demonstrate RhoPop's ability to recover known hyperparameters.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 Sect. 1.
Our 1-population model for the Luque & Pallé (2022) sample finds a central composition of 0.09 WMF with an intrinsic 1σ scatter of σ ρ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 2 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 2-population model.We find a higher density population with a central composition of WMF = 0.02 and σ ρ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 don't 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 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 towards 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 TRAPPIST-1 planets
The primary goal of Sect.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 our detection of 2-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 1-and 2-population models.Without these seven planets, the detection of 2-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 2-populations of small planets around M dwarfs.Further, for the N c = 2 scenario, we find a median value of µ 1 comp = 0.327 CMF, nearly identical to the CMF of Earth, in contrast to µ 1 comp = 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 irondepleted 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.

Alternative explanations for water worlds
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 − 4.2× Earth's.As such, it is likely substantial fractions of their water exist in the vapor and super-critical phases.Our choice of condensed water/ice then overestimates the density of these water layers, and, 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 super-critical 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 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 favors one scenario over the others.

Caveats to estimates on the minimum number of planets needed to resolve a Fe-rich population
In Sect.4, we estimated that approximately 154 planets with average σ Mp = 10% and σ Rp = 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 on the NASA Exoplanet Archive 6 with density uncertainties less than 100%.The average and median M-R uncertainties of these planets, however, are σ Mp = 27% and σ Rp = 6% and σ Mp = 19% and σ Rp = 5%, respectively.These uncertainties are ≳ 2× the precisions required to resolve a 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 1 Fe-rich planet for every 3 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 µ ml comp = 0.7.While we use Mercury and the 5 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.

CONCLUSIONS
We present the open-source software RhoPop for identifying distinct populations of small planets and use it to validate the recent results of Luque & Pallé (2022) and Adibekyan et al. (2021).While we show that a distinct population of Fe-rich planets is likely unresolvable with the current sample of small planets, RhoPop is designed to provide the community with a ready-to-go tool for identifying such populations as the number and precisions of small planets continue to increase.The next release of RhoPop will include pre-loaded compositional grids with atmosphere layers comprised of various molecular species.We will use these additional grids to investigate whether the population of planets interpreted by Luque & Pallé (2022) as 50% water worlds are better explained as having an atmosphere.In future work, we will apply RhoPop to the entire available sample of small planets and expand our analysis on the minimum number of planets required to resolve a population of iron-rich planets.

Figure 1 .
Figure 1.Water grid (cyan) and reduced rock grid (black).The ρ scaled parameter corresponds to the density at mass of a water-and volatile-free rocky planet with a liquid pure-iron core and Fe-free silicate mantle with Fe/Mg = 0.71, Si/Mg = 0.79, Ca/Mg = 0.09, and Al/Mg = 0.07.

Figure 2 .
Figure 2. HBM framework.Hyperparameters are shown in blue, local parameters in white, and observed parameters in gray.The observed and local parameters always have Np members corresponding to the number of planets in the sample.Nc = 1, 2, or 3 when the 1-, 2-, and 3-population model is considered, respectively.Observed masses and density ratios are modeled as Gaussian random variables centered on their unknown true values (local parameters) with a standard deviation equal to their observational uncertainties (Eq. 2 and 6, respectively).The local and hyperparameters are linked through the mixture model defined in Eq. 5.

Figure 3 .
Figure 3. Density ratio as a function of mass for our Nc = 1 (left) and Nc = 2 (right) models applied to the Luque & Pallé (2022) sample of small transiting planets in orbit around M-dwarfs.We color code the highest-density population as green and the lowest-density population as purple.The planets are assigned the same color as the population to which they most likely belong.The solid lines are the isocomposition lines for each µ j comp .The bands are the corresponding intrinsic density ratio scatter parameter σ j ρ ratio,pop .For all µ j comp and σ j ρ ratio,pop we use the median values given in Table3.

Figure 4 .
Figure 4. Density ratio as a function of mass for the preferred Nc = 1 model (left) and the Nc = 2 model (right) applied to the Adibekyan et al. (2021) sample.We color code the highest-density population as green and the lowest-density population as purple.The planets are assigned the same color as the population to which they most likely belong.The solid lines are the isocomposition lines for each µ j comp .The bands are the corresponding intrinsic density ratio scatter parameter σ j ρ ratio,pop .For all µ j comp and σ j ρ ratio,pop we use the median values given in Table4.

Figure 6 .
Figure 6.Summary of results for a synthetic 20-planet sample of planets with two-populations: a Mercury-like one with µ ml comp = 0.7 CMF and w ml = 0.25 and an Earth-like population with µ el comp = 0.325 CMF and w el = 0.75.

Figure 7 .
Figure 7. hyperparameter posterior corner plot for the 1-population model applied to our synthetic sample.The +/-values correspond to the 95% credible intervals.This figure was generated with the dynesty dynamic nested sampling software package (Speagle 2020).

Figure 8 .
Figure 8. hyperparameter posterior corner plot for the 2-population model applied to our synthetic sample.The +/-values correspond to the 95% credible intervals.This figure was generated with the dynesty dynamic nested sampling software package (Speagle 2020).

Table 1 .
Summary of parameter definitions.

Table 3 .
Summary of parameter posteriors, model evidence, and model selection parameters for our 1-and 2-population model runs with the Luque & Pallé (2022) sample.CI = credible interval

Table 4 .
Summary of parameter posteriors, model evidence, and model selection parameters for our 1-and 2-population model runs with the Adibekyan et al. (2021) sample.CI = credible interval