EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A253 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A253
Number of page(s) 11
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202452904
Published online 25 March 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Gas dwarfs, planets smaller than Neptune with substantial H/He envelopes and radii between 1.5 and 4 Earth radii, represent a prevalent population of exoplanets in the stellar neighbourhood, standing in stark contrast to their absence within our solar system (Borucki et al. 2010; Howard et al. 2012; Rogers 2015; Fulton et al. 2017). Within this size range, another hypothesised population of planets is found, often called ‘water worlds’, which are theorised to harbour significant quantities of water in condensed or steam form (Zeng et al. 2019; Unterborn et al. 2018; Luque & Pallé 2022).

The potential presence of water worlds gains traction with the study by Luque & Pallé (2022), who analyses a population of 34 small planets orbiting M dwarfs. Using composition models from Zeng et al. (2019), they identify six candidates that could possess substantial condensed water layers. However, subsequent work by Rogers et al. (2023) introduces a layer of complexity by suggesting that the observed characteristics of these planets could instead be explained by atmospheric boil-off processes. Their models propose that planets initially endowed with H/He atmospheres may shed their outer layers, ultimately becoming stripped rocky cores. Rogers et al. argue that current mass and radius measurements alone are insufficient to establish the presence of a water world population around M dwarfs. Further insights are provided by Chakrabarty & Mulders (2024), who simulate planetary system evolution using formation and atmospheric escape models to reproduce the observed radius valley, a gap between rocky planets and gas dwarfs in the massradius diagram. Their findings indicate that 20–35% of M-dwarf planets without primordial H/He atmospheres could be waterice-rich, a fraction consistent with the upper limit proposed by Luque & Pallé (2022). Building upon this work, Parviainen et al. (2024) reanalyse an updated version of the dataset used in Luque & Pallé (2022) and find no statistical support for a distinct population of water worlds. Their results highlight the sensitivity of water world classification to the choice of theoretical density models, which emphasises that the current data sample is insufficient for definitive conclusions. They also stress the need for a standardised definition of water worlds to enable meaningful comparisons across studies. In addition to these efforts, Parc et al. (2024) conducts a comprehensive analysis using the updated PlanetS catalogue of transiting planets. Their findings align with those of Parviainen et al. (2024), and concludes that no distinct population of water worlds exists around M dwarfs. However, they observe a trend where the minimum mass of gas dwarfs increases with earlier spectral types across the M to F star range, which provides important insights into planetary formation processes across different stellar environments. These findings further reinforce the notion that the classification of water worlds remains uncertain and highly model-dependent.

In this work, we adopt a different approach by not relying on a specific compositional definition of water worlds. Instead, we aim to identify an isolated cluster of planets distinct from both the rocky planet and gas dwarf populations, seeking robust statistical evidence through data-driven classification techniques. In this study, we investigated the existence of a water world population within an updated sample of small planets by analysing their mass, radius, and density measurements using three different methods.

In Section 2, we described our sample and selection criteria. In Section 3, we applied unsupervised clustering techniques to the mass-density diagram to identify the number of distinct planetary populations and delineate their boundaries. In Section 4, we applied ExoMDN, a machine-learning-based interior inference model developed by Baumeister & Tosi (2023). This model allowed us to perform a classification based on the inferred internal structure of the planets, considering physical layers such as the core, mantle, water, and atmosphere rather than relying solely on mass-radius or mass-density measurements. By incorporating the posterior distributions of these physical layers, our classification leveraged the underlying physics encoded within the ExoMDN model. Finally, in Section 5, we discussed our findings and summarised the key conclusions of our analysis.

Table 1

Summary of the selection criteria of the 2022 and 2025 samples.

2 The 2025 sample

In their work, Luque & Pallé (2022) compiles a dataset of small planets around M-type host stars, which we here refer to as the 2022 sample1. We started with systems available on the NASA Exoplanet Archive as of January 8th, 2025 and applied similar selection criteria as Luque & Pallé (2022). See Table A.1 in the appendix, where the planets new to the 2025 sample are marked with an asterisk. The selection criteria for the two samples are listed in Table 1. When density values were not available in the listed references, we calculated them along with their uncertainties using the equations provided in Appendix B.

Figure 1 presents the samples in three panels. Panel A illustrates the masses and radii of planets in the 2022 sample, while panel B shows the masses and radii of planets in the 2025 sample. Panel C displays the masses and densities of the 2025 sample planets, comparing their bulk densities to an Earth-like composition model derived from Zeng et al. (2019)2. All planets below the grey line would have a lower density than a scaled Earth would have, given the model. The model by Zeng et al. (2019) assumes mass fractions of 32.5% iron and 67.5% silicates. In this work, we refer to this model as the Earth composition model, ECM for short, and the scaling of density by this model is referred to as the Earth composition model density, pECM.

thumbnail Fig. 1

2025 sample. Panel A shows the 2022 sample’s planets plotted by mass and radius, while Panel B displays the 2025 sample, with new planets added in 2025 marked in light grey and planets also present in the 2022 sample marked in dark grey. Panel C illustrates the masses and densities of the 2025 sample, comparing the planetary bulk densities to an Earth-like composition model based on Zeng et al. (2019).

3 Population clustering in mass-density diagram

In this section, we sorted the 2025 sample into populations based on mass and density. Our approach aims to determine whether the sample is best described by two or three planet populations. This work defines a planet population as a group of planets sharing similar characteristics in a specific parameter space, such as mass-density.

3.1 Gaussian mixture modelling

We started by applying Gaussian mixture models (GMM), a probabilistic machine learning approach that assumes the planets originate from a mixture of a finite number of Gaussian distributions, each characterised by distinct means and variance (Dempster et al. 1977; Xu & Jordan 1996). For this purpose, we used the GMM implementation available in the scikit-learn Python package (Pedregosa et al. 2011)3 This clustering algorithm is applied to group the planets based on their ρECM values by effectively fitting Gaussian distributions to the sample histogram. Using GMM enabled us to categorise the 2025 sample without requiring prior labelling. Instead, we specified the number of populations, allowing the algorithm to fit the corresponding number of Gaussian components to the dataset. To determine appropriate initialisation parameters, we first applied the K-Means algorithm to set suitable initial means and precisions. K-Means is an unsupervised clustering algorithm that assigns data points to one of K clusters based on their distance to the cluster centres (centroids) (Lloyd 1982; MacQueen 1967). While the GMM itself does not initialise randomly, its initialisation step contains inherent randomness due to the K-Means algorithm. To ensure we identify the best-performing model, we performed 100 initialisations during model training. GMM does not inherently take parameter uncertainties into account. To address this, we implemented a Monte Carlo sampling step during training. In this step, we created a new sample where each planet is drawn from a normal distribution based on its uncertainties (the mean of the upper and lower uncertainty) 100 times. These generated samples were used to train the algorithm, allowing it to account for uncertainties and improve the robustness of the clustering results.

We applied the GMM algorithm twice: once assuming the sample consists of two populations and again assuming it consists of three populations. The results for the two-population and three-population cases are shown in panels A and B of Figure 2, respectively4.

We used two metrics to evaluate the GMMs’ performance against each other and whether the 2025 sample is better described by two or three clusters. The first criterion we used is the Akaike information criterion (AIC). It measures the tradeoff between the model’s goodness of fit and its complexity (the number of clusters). The second metric is the Bayesian information criterion, BIC, which also includes the sample size. The lower the AIC and BIC, the more preferred the model is for this metric. The two criteria are described as follows: AIC=2p−2 ln(L^),$AIC = 2p - 2\ln (\hat L),$(1) BIC=p ln(n)−2 ln(L^),$BIC = p\ln (n) - 2\ln (\hat L),$(2)

where n is the number of planets in the sample, p is the number of parameters learned by the model, and L is the maximised value of the likelihood function of the model (e.g., Geron 2019). The clustering analysis indicates that the two-population model performs best according to both the AIC and BIC criteria. Specifically, for the GMM with two clusters, we obtain AIC = 9.37 and BIC = 19.76, whereas for the GMM with three clusters, AIC increases to 14.05 and BIC to 30.67. This means that the GMM model suggests that the planetary densities in this sample are best described by two distinct populations rather than three.

thumbnail Fig. 2

Gaussian Mixture Modelling. The 2025 sample is analysed using Gaussian Mixture Modelling, with planets divided into 21 evenly spaced bins between 0 and 1.5. Panel A assumed two populations, while panel B assumed three populations. Evaluation metrics for the fits are displayed in the top-right corner. GJ 367b is excluded from the analysis due to its high density.

3.2 Population density based clustering

With GMM, our analysis focused solely on ρECM, requiring us to predefine the number of populations and subsequently compare the results. To overcome this limitation and categorise the 2025 sample into different populations without assuming their number in advance, we employed density-based spatial clustering of applications with noise (DBSCAN) (Ester et al. 1996; Sander et al. 1998). The method is implemented using the scikit-learn Python package5. Unlike many other methods, DBSCAN does not assume any particular cluster shape and can handle outliers effectively, i.e., planets not categorised into a population (Khan et al. 2014; Schubert et al. 2017). We opted to apply DBSCAN to the dataset in the log10 mass versus ρECM plane, as this parameter space was used in the analysis by Luque & Pallé (2022), and it allows for clustering without the need for physically arbitrary normalisation. DBSCAN identifies clusters as continuous regions of high density, separated by lower-density regions. DBSCAN has two key parameters, epsilon, and MinPts. The former determines the radius around each data point to assess density, while the latter sets the minimum number of points required to form a dense region.

In order to qualify and record how well the DBSCAN model performs, we used two evaluation metrics: silhouette score (Rousseeuw 1987) and Calinski–Harabasz index (Caliński & Harabasz 1974).

The silhouette score is a measure of how well the planets averagely fit into their own cluster, compared to their nearest neighbour cluster. It is the mean of all planet’s silhouette coefficients, which is calculated as:  silhouette coefficient =b−amax(a,b),${\rm{silhouettecoefficient}} = {{b - a} \over {\max (a,b)}},$(3)

here, a represents the mean distance to a particular other planet in the same cluster (giving a measure of cohesion), and b is the mean distance to the centre of the nearest neighbouring cluster (giving a measure of separation). The max function simply picks the largest value between a and b. The silhouette coefficient ranges from −1 to +1. A value near +1 indicates that a planet is well within its assigned cluster and far from others, whereas a value near 0 suggests proximity to a cluster boundary. A negative value implies potential misclassification. The silhouette score is a mean of all silhouette coefficients in the sample and can thereby also vary between −1 and +1, with larger silhouette scores representing a more successful cluster representation of the sample (Geron 2019).

The second evaluation metric, the Calinski-Harabasz index (CH index), is also used to evaluate how well the planets are grouped within each cluster and how distinct the clusters are from one another. It measures the balance between cohesion, which refers to how closely packed the planets are within a cluster, and separation, which reflects how far apart the clusters are.  CH-Index = cohesion k−1 separation −1n−k,${\rm{CH - Index}} = {{{\rm{cohesion}}} \over {k - 1}}{{{\rm{separation}}{{\rm{}}^{ - 1}}} \over {n - k}},$(4)  cohesion =∑i=1kni‖ ci−c ‖2,${\rm{cohesion}} = \mathop \sum \limits_{i = 1}^k {n_i}{{\bf{c}}_i} - {\bf{c}}{^2},$(5)  separation =∑i=1k∑x∈Cink‖ x−ck ‖2.${\rm{separation}} = \mathop \sum \limits_{i = 1}^k \mathop \sum \limits_{x \in {C_i}}^{{n_k}} x - {{\bf{c}}_k}{^2}.$(6)

In these equations, k represents the number of clusters, and n is the total number of planets in the sample. The term c is the centroid of all planets in the sample. For each specific cluster i, the number of planets is ni, the set of planets in the cluster is Ci, and x represents an individual planet within the cluster. The centroid of cluster i is denoted as ci. A higher CH index value indicates that the clusters are compact and well separated. However, since this is a relative measure, there is no fixed threshold that defines a good clustering result (Liu et al. 2010).

The two parameters, epsilon and MinPts, are set before the learning process begins, and their values influence the learning algorithm’s behaviour. Therefore, it is important to find appropriate values for these two parameters before clustering the 2025 sample with DBSCAN. We, therefore, performed a grid search. Specifically, we tested 18 values for MinPts from 3 to 20 and 200 values for epsilon from 0.1 to 0.4, resulting in 3600 combinations, and only kept the results that split the sample into at least two populations6. For each combination of MinPts and epsilon, we recorded both the Silhouette Score and the CH index. To account for the uncertainties in the 2025 sample, we applied bootstrapping, drawing from the uncertainty distributions in mass and density. Specifically, we took 200 draws for each MinPts-epsilon combination, which resulted in a total of 720 000 iterations of DBSCAN. The combined results for the Silhouette Scores and the CH index are presented as heatmaps in Figure 3. Panel A shows the Silhouette Scores, while Panel B presents the CH index. The combination yielding the highest scores for both metrics is considered optimal. Across both metrics, the best-performing parameters are an epsilon of 0.18 and a MinPts of 9. We have marked this optimal combination with a dark red ‘x’ in both panels of Figure 3.

In Figure 4, the DBSCAN classification is plotted with three select combinations of the MinPts and epsilon parameters. Each panel corresponds to one combination of MinPts and its bestperforming epsilon value. In the figure, gas planets are depicted in beige, while rocky planets are shown in brown. Due to the nature of the DBSCAN algorithm, many planets may be identified as outliers, denoted in grey. These outliers are planets that do not fall within high-density regions and are not categorised into the two populations. This does not mean that the outlying planets do not potentially belong to a population; rather, it indicates that they are not situated within the dense regions where most other planets are found on the diagram. The fit with the best-performing parameters is seen in panel B of Figure 4. Here, it is seen that the 2025 sample is clustered into two populations, divided by a number of outliers.

In panel A of Figure 4, a fit with a MinPts parameter set to 7 is displayed, which is 2 below the optimal parameter found, which is 9. Here, we observe that with lower MinPts and epsilon values, fewer planets are clustered into populations, and more planets are considered outliers. This outcome is expected because when epsilon, the radius around each data point where density is assessed, is smaller, only planets in the highest density areas are included in the clusters. Consequently, fewer planets are classified into populations, resulting in more planets being designated as outliers. However, this clustering also finds two planet populations, not three. In panel C, a fit with the MinPts parameter of 11 is displayed. Continuing on the same trend, these panels show that choosing a higher MinPts and epsilon value increases the number of planets in the clusters, thereby decreasing the number of outliers.

We find that the 2025 sample consistently clustered into two distinct populations with some outliers in between. Using various combinations of epsilon and MinPts around the best-performing values does not change this result. The 2025 sample consistently contains two, not three, clusters with outliers in between.

A graphical summary of the results from this section using the different clustering methods is given in Figure 5. Here, the 60 planets of the 2025 sample are colour-coded based on the cluster they are assigned to using different methods. As shown in the figure, TOI-776 b and TOI-1695b are classified as rocky planets using the GMM assuming two populations but are identified as gas dwarfs by the population density-based clustering. With the exception of these two planets, all other planets classified as either rocky planets or gas dwarfs by the population density-based method receive the same classification when using the two-population GMM. Although the GMM assuming three populations also identifies 11 water worlds, the two-population model is preferred due to its superior performance based on both the AIC and BIC.

thumbnail Fig. 3

Grids showing the evaluation of the 2025 sample’s DBSCAN classification. The evaluation was performed in ρECM-log10(mass) coordinates with different MinPts and epsilon values. The evaluation is a collection of bootstrapping runs with 200 samples. For each bootstrap sample, two grids of evaluation scores are made by running the DBSCAN on the bootstrapped sample over a range of 200 epsilon values from 0.1 to 0.4 and MinPts from 3 to 20. One grid records the silhouette score, and one grid records the CHI. These two grids are used to find the optimal parameters for the 2025 dataset. In panel A, the total grid of all the recorded silhouette scores. In panel B, the total grid of all the recorded CHI. In both of the top panels, the best-performing epsilon value is marked with a red x.
thumbnail Fig. 4

DBSCAN classification of the 2025 sample in ρECM-log10(Mass) space, using the best-performing combinations of different MinPts and epsilon values. Each panel has a different combination of MinPts and epsilon parameters, indicated by the values in the panel. In panel A is seen the fit of the optimal epsilon for MinPts = 7. In panel B, the overall optimal hyperparameters were found from grid search, with MinPts = 9 and epsilon = 0.18, is found. In panel C, the optimal epsilon for MinPts = 11, is shown. The planets are coloured based on their clusters. Brown planets are rocky planets, and beige planets are gas-rich planets. Grey planets are classified as outliers. The densities are normalised by the Earth composition model (ECM, grey line).

4 Interior structure-based clustering

So far, our planet clustering has been based on mass, radius, and density. These parameters can be inferred directly from the observations. Here, we incorporate theoretical insights by applying a planetary structure model to infer the interior structure of each of the 60 planets. We then perform clustering on the internal structure found for each planet instead of on the mass radius or density measurements.

To model the interiors of the planets in the 2025 sample, we use ExoMDN developed by Baumeister & Tosi (2023). ExoMDN is an exoplanet interior inference model employing a mixture density network (MDN). This neural network outputs a posterior distribution of mixed Gaussians, allowing us to predict probability distributions for the individual planet’s interior structure based on its mass, radius, and equilibrium temperature. ExoMDN is trained using synthetic data generated by the TATOOINE code Baumeister et al. (2020); MacKenzie et al. (2023). TATOOINE divides a planet into four distinct layers: an iron core, a silicate mantle, a water-ice shell, and a hydrogen-helium gas layer. Each layer is characterised by its radius fraction, with the sum of the radii of all layers required to total the planet’s radius. In the future, one might want to incorporate the mixing of water into deeper layers of the planet. Luo et al. (2024) proposes the majority of the bulk water can be stored deep within the core and the mantle and not as only as a layer at the surface, specifically for more massive planets (> 6M⊕) and Earth-sized planets with smaller amounts of water.

Here, the core is assumed to consist entirely of molten, pure hcp-iron. The silicate mantle is segmented into an upper and lower part, distinguished by a phase transition occurring at 23 GPa. For the water layers, a composition of pure H2O is assumed, omitting additional compounds like methane or ammonia. The Equations of State (EoS) adopted for modelling span a vast range from 0.1 to 400 TPa and from 150 to 105 K, encompassing gas, liquid, and solid phases of water (Haldemann et al. 2020). Both liquid and solid layers are considered fully convective, with an adiabatic temperature profile, while water vapour is integrated into the atmosphere. The atmospheric layer is characterised as an isothermal gaseous H/He envelope, mirroring solar-like composition (71% hydrogen and 29% helium by mass), based on the EoS from Saumon et al. (1995), with a temperature equivalent to the equilibrium temperature.

ExoMDN outputs a single multi-dimensional posterior distribution, where each dimension corresponds to the thickness of a layer. Each layer is comprised of 20 mixture Gaussians, each characterised by a mean, spread, and mixing weight.

To categorise the 60 planets based on their four posterior planet layer distributions, we drew a discrete planet composition from each planet’s posterior distribution, ensuring that the thickness of each layer – core, mantle, water, and atmosphere – adds up to 1: 1=XCore +XMantle +XWater +XAtmosphere .$1 = {X_{Core{\rm{}}}} + {X_{Mantle{\rm{}}}} + {X_{Water{\rm{}}}} + {X_{Atmosphere{\rm{}}}}.$(7)

The posterior distributions of the planets in the 2025 sample are shown in Figure C.1. Using the predicted thicknesses of each planet, we applied the KMeans7 algorithm, implemented in scikit-learn, to group the planets into distinct populations. We ran the KMeans algorithm ten times to find the optimal clustering for the four-layer thickness predictions. From these, we selected the best-performing as the clustering with the lowest measure of the difference between points within the same group, Inertia.

We repeated this procedure 1000 times, i.e., randomly selecting a discrete makeup for each of the 60 planets for each iteration. This iterative process enabled us to identify which planets are most challenging to cluster correctly and provides a relative likelihood of a planet being in a particular cluster. We consider a planet part of a population if it is assigned to the same population in 900 out of 1000 iterations. We performed this overall routine twice, assuming two or three populations.

Figure 6 displays an overlay of 100 iterations. We assigned a colour to each planet group based on the average density of the planets within the group. Specifically, the group with the lowest ρECM is coloured yellow, the group with the middle mean ρECM is coloured blue, and the group with the highest mean ρECM is coloured red. For the three populations case (panel B) the routine exclusively classifies rocky planets, but cannot distinguish between the lower and middle clusters with sufficient accuracy to classify a single planet into one of these populations. The routine performs substantially better for the two population case (panel A), where a population of rocky planets and a population of gas planets are identified.

In Figure 7, the posterior thickness layer distributions of the 2025 sample are shown. The planets are sorted top-down by the mode of their predicted core fraction. The planets are colour-coded according to their population, as determined by the interior-based two-population K-Means clustering.

There is a clear distinction between the gas planets at the top and the rocky planets at the bottom. This paper’s work so far has been to try to draw these limits and find a division between populations. However, in this figure, no definitive horizontal line separates the two populations. Instead, we observe a gradual evolution in the core fraction of the planets. The top gas planets have small core radius fractions, and the bottom terrestrials are mostly made up of the core. Likewise, most of the large gas planets at the top have extensive atmospheres, while the bottom rocky planets have almost no envelope. The top gas planets also may have an extended water layer, while the bottom terrestrials do not.

thumbnail Fig. 5

Summary of clustering results. Each row represents a planet, and each column corresponds to a different clustering method. GMM2 is categorised using Gaussian mixture modelling, assuming two populations in the 2025 sample, and GMM3 is categorised using Gaussian mixture modelling, assuming three populations in the 2025 sample. DBSCAN is the population density based clustering using the DBSCAN algorithm on the 2025 sample.

5 Synthesis and conclusions

In this work, we have followed up on recent findings by Luque & Pallé (2022), who identify three distinct populations among a sample of 34 planets. We extended this work by applying three classification methods to an updated sample that includes 27 additional planets, selected using criteria similar to those in Luque & Pallé (2022).

We employed the following three methodologies to investigate the possible presence of a distinct water world population in the 2025 sample:

  1. Gaussian mixture modelling (GMM): this method, detailed in Section 3.1, assumed that planets are divided into two or three populations, each distributed around a mean density normalised by the Earth Composition Model (ρECM) and scattered in a Gaussian manner. GMM favours the existence of two populations evaluated by the AIC and BIC. However, it relies solely on the ρECM parameter, which has been questioned by Rogers et al. (2023) for its applicability in this context.

  2. Population density-based clustering (DBSCAN): introduced in Section 3.2, this approach allowed for flexible population identification in the mass-density plane without assuming a particular shape or requiring the number of populations as an input. Unlike GMM, DBSCAN incorporated the planets’ mass uncertainties using a bootstrapping grid search to optimise clustering parameters. The method robustly predicts the presence of two populations with a number of outliers, including L 98-59 d, Kepler-138 c and TOI-244b, which all are categorised as water worlds by the GMM assuming three populations. However, its reliance on normalisation choices and sensitivity to density variations within populations pose limitations. With an optimal MinPts value of 9 in our sample of 60 planets, very small populations, such as a potential water world population, may not be accurately identified.

  3. Interior structure-based clustering (ExoMDN): in contrast to the other methods, this approach (discussed in Section 4) did not rely on ρECM but instead uses the ExoMDN machinelearning model (Baumeister & Tosi 2023), which infers planetary interior compositions from mass, radius and equilibrium temperature. Clustering was performed using KMeans based on the inferred interior layers: core, mantle, water, and atmosphere. When using K = 2, the results consistently revealed two distinct populations with outliers that overlap with those identified by DBSCAN. When K = 3, the method can only distinguish rocky planets, with the remaining planets frequently switching between the two other clusters. This suggests that the sample lacks sufficient statistical support to confirm the presence of a distinct third population.

Among the three methods, the planets TOI-777 b and TOI-1695 b are classified differently. The population density-based clustering approach identifies them as outliers, whereas DBSCAN classifies them as gas dwarfs, and Gaussian Mixture Modelling categorises them as rocky planets. Overall, both the population density-based clustering and interior structure-based clustering methods yield consistent findings, identifying two dominant populations separated by a set of outliers. The agreement between these independent methods reinforces the conclusion that the planets in this sample do not form three distinct populations.

Our findings align with recent work by Parviainen et al. (2024), Parc et al. (2024), and Rogers et al. (2023), all of whom report no compelling evidence for a distinct water world population among M-dwarf planets. While Luque & Pallé (2022) initially suggest the presence of water worlds, our analysis using Gaussian Mixture Modelling, population density-based clustering, and interior structure-based clustering consistently identifies only two populations within the 2025 sample. If water worlds exist, their inferred properties do not form a statistically distinct cluster. This suggests that they either significantly overlap with other planet types or that the current sample lacks the resolution needed to detect them as a separate group. A larger, well-characterised sample, combined with improved interior modelling techniques, will be essential to further investigate the potential existence of water worlds.

thumbnail Fig. 6

KMeans classification applied to the planet layer thickness fractions derived from ExoMDN for the 2025 sample. (A) Classification performed with two assumed populations. (B) Classification with three assumed populations. The fitting algorithm was executed 100 times and stacked on top of each other. For each run, a new sample is made by randomly drawing a new composition for each planet from its four-layer posterior distribution, such that the total thickness of a planet adds up to 1. KMeans then clusters each drawn sample, and the planets are coloured according to their assigned cluster. The top cluster (cluster in which the mean density of all member planets is the highest) is coloured red, the bottom cluster is coloured yellow, and the middle (for three populations) is coloured blue. Planets that do not match the three primary colours have been assigned to different populations in various iterations.
thumbnail Fig. 7

Thickness layer predictions of the planets in the 2025 sample. The planets are sorted top-down by the mode of their predicted core fraction. The plots present the rankings for relative size, core, mantle, water, and atmosphere. The planets are coloured by the KMeans clustering method using two populations, with planets coloured based on their inclusion in the population drawn in at least 90% of clustering. To see the planet names of each set of posteriors, see the Figure C.1.

Acknowledgements

We thank Philipp Bauenmeister and Rafael Luque for their thoughtful comments and discussions, which improved this manuscript. We would also like to thank the anonymous referee for providing constructive and detailed feedback on the manuscript, which greatly enhanced its quality. We acknowledge the support from the Danish Council for Independent Research through a grant, No.2032-00230B. 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 work made use of the Python packages scikit-learn (Pedregosa et al. 2011), and ExoMDN (Baumeister & Tosi 2023).

Appendix A Sample data

Table A.1

Planets in the 2025 sample. Planets new to the 2025 sample and thus not in the 2022 sample are marked with an asterisk.

Appendix B Density calculations

We have calculated the densities and their uncertainties for planets that do not have listed densities by using standard uncertainty propagation. In the equations below, we describe the methods used to compute the densities, as well as the respective upper and lower uncertainties. ρ=3M4πR3,$\rho = {{3M} \over {4\pi {R^3}}},$(B.1) (Δρ)2=(δρδMΔM)2+(δρδRΔR)2=(34πR3ΔM)2+(9M4πR4ΔR)2.${({\rm{\Delta }}\rho )^2} = {\left( {{{\delta \rho } \over {\delta M}}{\rm{\Delta }}M} \right)^2} + {\left( {{{\delta \rho } \over {\delta R}}{\rm{\Delta }}R} \right)^2} = {\left( {{3 \over {4\pi {R^3}}}{\rm{\Delta }}M} \right)^2} + {\left( {{{9M} \over {4\pi {R^4}}}{\rm{\Delta }}R} \right)^2}.$(B.2)

For asymmetric (upper and lower) uncertainties, we calculate: ρ+=(34πR3M+)2+(9M4πR4R−)2,${\rho _ + } = \sqrt {{{\left( {{3 \over {4\pi {R^3}}}{M_ + }} \right)}^2} + {{\left( {{{9M} \over {4\pi {R^4}}}{R_ - }} \right)}^2}} ,$(B.3) rho−=(34πR3M−)2+(9M4πR4R+)2.$rh{o_ - } = \sqrt {{{\left( {{3 \over {4\pi {R^3}}}{M_ - }} \right)}^2} + {{\left( {{{9M} \over {4\pi {R^4}}}{R_ + }} \right)}^2}} .$(B.4)

Where:

  • M is the mass of the planet, and ΔM is its uncertainty. M+ and M− denote the upper and lower uncertainties (errors) in the planet’s mass, respectively.

  • R is the radius of the planet, and ΔR is its uncertainty. R+ and R− denote the upper and lower uncertainties (errors) in the planet’s radius, respectively.

  • ρ is the calculated bulk density of the planet, and Δρ represents the (symmetric) uncertainty in the calculated density.

  • ∂ρ∂M and ∂ρ∂R${{\partial \rho } \over {\partial M}}{\rm{and}}{{\partial \rho } \over {\partial R}}$ denote the partial derivatives of the density with respect to mass and radius, respectively.

In summary, equation B.2 gives the density, equation B.2 shows how to propagate symmetric uncertainties using partial derivatives, and equations B.3 and B.4 illustrate how to obtain the upper and lower bounds on the density, when the uncertainties in mass and radius are asymmetric.

Appendix C Posterior distributions of the 2025 sample

thumbnail Fig. C.1

Planet layer predictions for the planets in the 2025 sample. The planets are sorted top-down by the mode of their predicted core fraction.

References

  1. Agol, E., Dorn, C., Grimm, S. L., et al. 2021, PSJ, 2, 1 [NASA ADS] [Google Scholar]
  2. Almenara, J. M., Bonfils, X., Otegi, J. F., et al. 2022, A&A, 665, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Astudillo-Defru, N., Cloutier, R., Wang, S. X., et al. 2020, A&A, 636, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Baumeister, P., & Tosi, N. 2023, A&A, 676, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Baumeister, P., Padovan, S., Tosi, N., et al. 2020, AJ, 889, 42 [Google Scholar]
  6. Benneke, B., Wong, I., Piaulet, C., et al. 2019, ApJ, 887, L14 [Google Scholar]
  7. Bluhm, P., Luque, R., Espinoza, N., et al. 2020, A&A, 639, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bonfanti, A., Brady, M., Wilson, T. G., et al. 2024, A&A, 682, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [Google Scholar]
  10. Brady, M., Bean, J. L., Seifahrt, A., et al. 2024, AJ, 168, 67 [Google Scholar]
  11. Burt, J. A., Dragomir, D., Mollière, P., et al. 2021, AJ, 162, 87 [NASA ADS] [CrossRef] [Google Scholar]
  12. Burt, J. A., Hooton, M. J., Mamajek, E. E., et al. 2024, ApJ, 971, L12 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cadieux, C., Plotnykov, M., Doyon, R., et al. 2024, ApJ, 960, L3 [NASA ADS] [CrossRef] [Google Scholar]
  14. Caliński, T., & Harabasz, J. 1974, Commun. Stat. Theory Methods, 3, 1 [Google Scholar]
  15. Castro-González, A., Demangeon, O. D. S., Lillo-Box, J., et al. 2023, A&A, 675, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Chakrabarty, A., & Mulders, G. D. 2024, AJ, 966, 185 [Google Scholar]
  17. Chaturvedi, P., Bluhm, P., Nagel, E., et al. 2022, A&A, 666, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cherubim, C., Cloutier, R., Charbonneau, D., et al. 2023, AJ, 165, 167 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cloutier, R., Charbonneau, D., Stassun, K. G., et al. 2021, AJ, 162, 79 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cloutier, R., Greklek-McKeon, M., Wurmser, S., et al. 2024, MNRAS, 527, 5464 [Google Scholar]
  21. Cointepas, M., Almenara, J. M., Bonfils, X., et al. 2021, A&A, 650, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Cointepas, M., Bouchy, F., Almenara, J. M., et al. 2024, A&A, 685, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Dai, F., Howard, A. W., Halverson, S., et al. 2024, AJ, 168, 101 [NASA ADS] [CrossRef] [Google Scholar]
  24. Demangeon, O. D. S., Zapatero Osorio, M. R., Alibert, Y., et al. 2021, A&A, 653, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, J. R. Stat. Soc. Ser. B Methodol., 39, 1 [Google Scholar]
  26. Diamond-Lowe, H., Kreidberg, L., Harman, C. E., et al. 2022, AJ, 164, 172 [NASA ADS] [CrossRef] [Google Scholar]
  27. Ducrot, E., Gillon, M., Delrez, L., et al. 2020, A&A, 640, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Essack, Z., Shporer, A., Burt, J. A., et al. 2023, AJ, 165, 47 [NASA ADS] [CrossRef] [Google Scholar]
  29. Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. 1996, in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (KDD’96), eds. D. W. Pfitzner, & J. K. Salmon, 226 [Google Scholar]
  30. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
  31. Geron, A. 2019, Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems, 2nd edn. (California: O’Reilly Media, Inc.) [Google Scholar]
  32. Goffo, E., Gandolfi, D., Egger, J. A., et al. 2023, ApJ, 955, L3 [NASA ADS] [CrossRef] [Google Scholar]
  33. Goffo, E., Chaturvedi, P., Murgas, F., et al. 2024, A&A, 685, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. González-Álvarez, E., Zapatero Osorio, M. R., Caballero, J. A., et al. 2023, A&A, 675, A177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Hamann, A., Montet, B. T., Fabrycky, D. C., Agol, E., & Kruse, E. 2019, AJ, 158, 133 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hobson, M. J., Bouchy, F., Lavie, B., et al. 2024, A&A, 688, A216 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [Google Scholar]
  39. Kemmer, J., Stock, S., Kossakowski, D., et al. 2020, A&A, 642, A236 [EDP Sciences] [Google Scholar]
  40. Khan, K., Rehman, S. U., Aziz, K., Fong, S., & Sarasvady, S. 2014, in The Fifth International Conference on the Applications of Digital Information and Web Technologies (ICADIWT 2014), 232 [Google Scholar]
  41. Kossakowski, D., Kemmer, J., Bluhm, P., et al. 2021, A&A, 656, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Lacedelli, G., Pallè, E., Luque, R., et al. 2024, A&A, 692, A238 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Lam, K. W. F., Korth, J., Masuda, K., et al. 2020, AJ, 159, 120 [NASA ADS] [CrossRef] [Google Scholar]
  44. Leleu, A., Delisle, J. B., Udry, S., et al. 2023, A&A, 669, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Liu, Y., Li, Z., Xiong, H., Gao, X., & Wu, J. 2010, in 2010 IEEE International Conference on Data Mining, 911 [Google Scholar]
  46. Livingston, J. H., Crossfield, I. J. M., Petigura, E. A., et al. 2018, AJ, 156, 277 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lloyd, S. 1982, IEEE Trans. Inf. Theory, 28, 129 [Google Scholar]
  48. Luo, H., Dorn, C., & Deng, J. 2024, Nat. Astron., 8, 1399 [Google Scholar]
  49. Luque, R., & Pallé, E. 2022, Science, 377, 1211 [NASA ADS] [CrossRef] [Google Scholar]
  50. Luque, R., Pallé, E., Kossakowski, D., et al. 2019, A&A, 628, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Luque, R., Serrano, L. M., Molaverdikhani, K., et al. 2021, A&A, 645, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Luque, R., Fulton, B. J., Kunimoto, M., et al. 2022a, A&A, 664, A199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Luque, R., Nowak, G., Hirano, T., et al. 2022b, A&A, 666, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. MacKenzie, J., Grenfell, J. L., Baumeister, P., et al. 2023, A&A, 671, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. MacQueen, J. B. 1967, in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability (Berkeley, CA, USA: University of California Press), 1, 281 [Google Scholar]
  56. Mahajan, A. S., Eastman, J. D., & Kirk, J. 2024, ApJ, 963, L37 [NASA ADS] [CrossRef] [Google Scholar]
  57. Muirhead, P. S., Hamren, K., Schlawin, E., et al. 2012, ApJ, 750, L37 [NASA ADS] [CrossRef] [Google Scholar]
  58. Murgas, F., Pallé, E., Orell-Miquel, J., et al. 2024, A&A, 684, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Oddo, D., Dragomir, D., Brandeker, A., et al. 2023, AJ, 165, 134 [NASA ADS] [CrossRef] [Google Scholar]
  60. Palle, E., Orell-Miquel, J., Brady, M., et al. 2023, A&A, 678, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Parc, L., Bouchy, F., Venturini, J., Dorn, C., & Helled, R. 2024, A&A, 688, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Parviainen, H., Luque, R., & Palle, E. 2024, MNRAS, 527, 5693 [Google Scholar]
  63. Pass, E. K., Winters, J. G., Charbonneau, D., et al. 2023, AJ, 166, 171 [NASA ADS] [CrossRef] [Google Scholar]
  64. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825 [Google Scholar]
  65. Peterson, M. S., Benneke, B., Collins, K., et al. 2023, Nature, 617, 701 [NASA ADS] [CrossRef] [Google Scholar]
  66. Piaulet, C., Benneke, B., Almenara, J. M., et al. 2023, Nat. Astron, 7, 206 [Google Scholar]
  67. Rogers, L. A. 2015, ApJ, 801, 41 [Google Scholar]
  68. Rogers, J. G., Schlichting, H. E., & Owen, J. E. 2023, ApJ, 947, L19 [NASA ADS] [CrossRef] [Google Scholar]
  69. Rousseeuw, P. J. 1987, J. Comput. Appl. Math., 20, 53 [Google Scholar]
  70. Sander, J., Ester, M., Kriegel, H.-P., & Xu, X. 1998, Data Min. Knowl. Discov., 2, 169 [Google Scholar]
  71. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  72. Schubert, E., Sander, J., Ester, M., Kriegel, H. P., & Xu, X. 2017, ACM Trans. Database Syst., 42, 1 [CrossRef] [Google Scholar]
  73. Shporer, A., Collins, K. A., Astudillo-Defru, N., et al. 2020, ApJ, 890, L7 [Google Scholar]
  74. Soto, M. G., Anglada-Escudé, G., Dreizler, S., et al. 2021, A&A, 649, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Stefansson, G., Mahadevan, S., Maney, M., et al. 2020, AJ, 160, 192 [NASA ADS] [CrossRef] [Google Scholar]
  76. Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018, Nat. Astron, 2, 297 [Google Scholar]
  77. Van Eylen, V., Astudillo-Defru, N., Bonfils, X., et al. 2021, MNRAS, 507, 2154 [NASA ADS] [Google Scholar]
  78. Weiner Mansfield, M., Xue, Q., Zhang, M., et al. 2024, ApJ, 975, L22 [NASA ADS] [CrossRef] [Google Scholar]
  79. Xu, L., & Jordan, M. 1996, Neural Comput., 8, 129 [Google Scholar]
  80. Xue, Q., Bean, J. L., Zhang, M., et al. 2024, ApJ, 973, L8 [NASA ADS] [CrossRef] [Google Scholar]
  81. Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, PNAS, 116, 9723 [Google Scholar]
1

L 168-9 b is included in the 2022 sample but excluded from the 2025 sample in this work. This exclusion is due to the planetary parameters reported by Hobson et al. (2024), which differ from those in Astudillo-Defru et al. (2020) (used in Luque & Pallé 2022), having uncertainties in the planet radius that exceed the criteria set for inclusion in this study.

2

This composition model has been made available at https://lweb.cfa.harvard.edu/~lzeng/planetmodels.html#mrrelation by Zeng et al. (2019).

4

We exclude GJ 367 b in both GMM runs. It has an extremely high density, 10.2 ± 1.3 g cm°, (Goffo et al. 2023) for its size (0.699 ± 0.024 R⊕) and therefore clearly does not represent a water world nor a gas dwarf. It would be classified as its own population, taking up its own Gaussian, and it is not relevant to our discussion here.

6

Fractionally, most of the saved combinations result in two populations. Specifically, 85% of the saved combinations yield two planet populations, while 10% produce three populations, 3% result in four populations, and 2% lead to five or more populations.

All Tables

Table 1

Summary of the selection criteria of the 2022 and 2025 samples.

In the text
Table A.1

Planets in the 2025 sample. Planets new to the 2025 sample and thus not in the 2022 sample are marked with an asterisk.

In the text

All Figures

thumbnail Fig. 1

2025 sample. Panel A shows the 2022 sample’s planets plotted by mass and radius, while Panel B displays the 2025 sample, with new planets added in 2025 marked in light grey and planets also present in the 2022 sample marked in dark grey. Panel C illustrates the masses and densities of the 2025 sample, comparing the planetary bulk densities to an Earth-like composition model based on Zeng et al. (2019).
In the text
thumbnail Fig. 2

Gaussian Mixture Modelling. The 2025 sample is analysed using Gaussian Mixture Modelling, with planets divided into 21 evenly spaced bins between 0 and 1.5. Panel A assumed two populations, while panel B assumed three populations. Evaluation metrics for the fits are displayed in the top-right corner. GJ 367b is excluded from the analysis due to its high density.
In the text
thumbnail Fig. 3

Grids showing the evaluation of the 2025 sample’s DBSCAN classification. The evaluation was performed in ρECM-log10(mass) coordinates with different MinPts and epsilon values. The evaluation is a collection of bootstrapping runs with 200 samples. For each bootstrap sample, two grids of evaluation scores are made by running the DBSCAN on the bootstrapped sample over a range of 200 epsilon values from 0.1 to 0.4 and MinPts from 3 to 20. One grid records the silhouette score, and one grid records the CHI. These two grids are used to find the optimal parameters for the 2025 dataset. In panel A, the total grid of all the recorded silhouette scores. In panel B, the total grid of all the recorded CHI. In both of the top panels, the best-performing epsilon value is marked with a red x.
In the text
thumbnail Fig. 4

DBSCAN classification of the 2025 sample in ρECM-log10(Mass) space, using the best-performing combinations of different MinPts and epsilon values. Each panel has a different combination of MinPts and epsilon parameters, indicated by the values in the panel. In panel A is seen the fit of the optimal epsilon for MinPts = 7. In panel B, the overall optimal hyperparameters were found from grid search, with MinPts = 9 and epsilon = 0.18, is found. In panel C, the optimal epsilon for MinPts = 11, is shown. The planets are coloured based on their clusters. Brown planets are rocky planets, and beige planets are gas-rich planets. Grey planets are classified as outliers. The densities are normalised by the Earth composition model (ECM, grey line).
In the text
thumbnail Fig. 5

Summary of clustering results. Each row represents a planet, and each column corresponds to a different clustering method. GMM2 is categorised using Gaussian mixture modelling, assuming two populations in the 2025 sample, and GMM3 is categorised using Gaussian mixture modelling, assuming three populations in the 2025 sample. DBSCAN is the population density based clustering using the DBSCAN algorithm on the 2025 sample.
In the text
thumbnail Fig. 6

KMeans classification applied to the planet layer thickness fractions derived from ExoMDN for the 2025 sample. (A) Classification performed with two assumed populations. (B) Classification with three assumed populations. The fitting algorithm was executed 100 times and stacked on top of each other. For each run, a new sample is made by randomly drawing a new composition for each planet from its four-layer posterior distribution, such that the total thickness of a planet adds up to 1. KMeans then clusters each drawn sample, and the planets are coloured according to their assigned cluster. The top cluster (cluster in which the mean density of all member planets is the highest) is coloured red, the bottom cluster is coloured yellow, and the middle (for three populations) is coloured blue. Planets that do not match the three primary colours have been assigned to different populations in various iterations.
In the text
thumbnail Fig. 7

Thickness layer predictions of the planets in the 2025 sample. The planets are sorted top-down by the mode of their predicted core fraction. The plots present the rankings for relative size, core, mantle, water, and atmosphere. The planets are coloured by the KMeans clustering method using two populations, with planets coloured based on their inclusion in the population drawn in at least 90% of clustering. To see the planet names of each set of posteriors, see the Figure C.1.
In the text
thumbnail Fig. C.1

Planet layer predictions for the planets in the 2025 sample. The planets are sorted top-down by the mode of their predicted core fraction.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.