Spright: a probabilistic mass-density-radius relation for small planets

We present Spright, a Python package that implements a fast and lightweight mass-density-radius relation for small planets. The relation represents the joint planetary radius and bulk density probability distribution as a mean posterior predictive distribution of an analytical three-component mixture model. The analytical model, in turn, represents the probability for the planetary bulk density as three generalised Student's t-distributions with radius-dependent weights and means based on theoretical composition models. The approach is based on Bayesian inference and aims to overcome the rigidity of simple parametric mass-radius relations and the danger of overfitting of non-parametric mass-radius relations. The package includes a set of pre-trained and ready-to-use relations based on two M dwarf catalogues, one FGK star catalogue, and two theoretical composition models for water-rich planets. The inference of new models is easy and fast, and the package includes a command line tool that allows for coding-free use of the relation, including the creation of publication-quality plots. Additionally, we study whether the current mass and radius observations of small exoplanets support the presence of a population of water-rich planets positioned between rocky planets and sub-Neptunes. The study is based on Bayesian model comparison and shows somewhat strong support against the existence of a water-world population around M dwarfs. However, the results of the study depend on the chosen theoretical water-world density model. A more conclusive result requires a larger sample of precisely characterised planets and community consensus on a realistic water world interior structure and atmospheric composition model.


INTRODUCTION
The launch of the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) has enabled the bulk density measurement of hundreds of exoplanets thanks to its all-sky observing strategy of nearby, bright stars.Its contribution is already comparable to that of the Kepler/K2 mission for small planets ( < 4  ⊕ ) with mass determinations.Out of the thousands of small exoplanets discovered by Kepler/K2, only a few hundred have bulk density measurements (∼ 330, based on the NASA Exoplanet Archive 1 as of May 2023).On the other hand, approximately 130 small planets discovered by TESS have precisely determined masses and radii to date, with hundreds of candidates awaiting to be confirmed and characterised with ground-based facilities.For M dwarfs, the contribution is even larger -approximately 7% of all the Kepler/K2 small planets with measured bulk densities orbit M-dwarf hosts, while for TESS, the ratio is 40%.However, such in-depth characterisation is observationally expensive and it becomes harder as the planet-host mass and size ratios decrease.
Probabilistic mass-radius (M-R) relationships are useful not only for the purpose of predicting one quantity from the other but also as a means of understanding exoplanet compositions.On the one ★ E-mail: hannu@iac.es(HP) 1 https://exoplanetarchive.ipac.caltech.edu/index.htmlhand, they allow feasibility studies and efficient planning of radial velocity (RV) and transmission spectroscopy observations of transiting planets, which require an estimate of the mass given a radius measurement to predict the expected RV semi-amplitude and the planet's atmosphere scale height (which is inversely proportional to the planet's gravity and relates to its detectability), respectively.Reversely, upcoming microlensing discoveries with, e.g., the Roman Space Telescope (Spergel et al. 2015) will have mass estimates for which a direct radius measurement is impossible.
On the other hand, M-R relationships are a robust tool to identify demographic features of the exoplanet population, such as the transition from brown dwarfs to hydrogen-burning stars (e.g., Hatzes & Rauer 2015), Neptunian to Jovian planets (e.g., Wolfgang et al. 2016;Chen & Kipping 2017;Bashi et al. 2017), and rocky to volatile-rich planets (e.g., Weiss & Marcy 2014;Zeng et al. 2019;Otegi et al. 2020;Luque & Pallé 2022).Linking these trends with the physical and chemical processes at play during planet formation and evolution offers an avenue to constrain observationally such theories.This is especially relevant for sub-Neptune-sized exoplanetswith radii between 1.5 and 4  ⊕ -whose nature and origin are actively debated.Unlikely to be rocky in nature (Rogers 2015;Fulton et al. 2017), these planets reside in a degenerate part of the M-R parameter space where their bulk densities are equally well explained by solid rocky/iron cores with primordial gaseous hydrogen-rich atmospheres (sometimes referred to as "gas dwarfs", Lopez & Fortney 2014;Rogers et al. 2023) or a water-rich interior and atmosphere akin to the icy moons of the solar system (sometimes referred to as "water worlds", Léger et al. 2004;Dorn & Lichtenberg 2021;Aguichine et al. 2021).Both gas dwarfs and water worlds are naturally explained by current planet formation and evolution models (see e.g., Lee & Chiang 2016;Owen & Wu 2017;Ginzburg et al. 2018;Bitsch et al. 2019;Venturini et al. 2020;Burn et al. 2021), but with remarkably different implications about their location at birth.While the prevailing view is that sub-Neptunes are primarily gas dwarfs (see the review by Bean et al. 2021), the existence of water worlds appears strongly supported by recent individual planet discoveries (Bluhm et al. 2021;Diamond-Lowe et al. 2022;Piaulet et al. 2023;Acuña et al. 2022), observational demographic studies (Zeng et al. 2019;Neil et al. 2022;Luque & Pallé 2022), and advances in interior structure and global formation modeling (Venturini et al. 2020;Burn et al. 2021;Dorn & Lichtenberg 2021;Aguichine et al. 2021).
Most of the previous studies on M-R relations have assumed that exoplanet masses and radii follow one or multiple power-law segments of the form  ∝   (Lissauer et al. 2011;Weiss & Marcy 2014;Wolfgang et al. 2016;Mills & Mazeh 2017;Chen & Kipping 2017;Bashi et al. 2017;Otegi et al. 2020).Others have proposed a non-parametric approach instead (Ning et al. 2018;Kanodia et al. 2019).The power-law models are simple to fit, and their parameters are easy to interpret, but their rigidity also means that they may give an overly simplistic representation of the actual M-R relation.On the contrary, non-parametric models can take on a variety of shapes to fit the data and do not assume the distribution of masses at a given radius to be Gaussian or even symmetric, but their flexibility plays against their precision for small sample sizes.
In this paper, we present spright,2 a Python package that provides a lightweight probabilistic M-R relation for small planets.The relation models the joint planetary radius and bulk density distribution as a mean of the posterior predictive distribution of a simple analytic three-component mixture model.The approach is based on basic Bayesian inference and aims to overcome the shortcomings of the existing methods by 1) delivering robustness and flexibility not offered by parametric models while 2) avoiding the dangers of overfitting and the need for large sample sizes associated with non-parametric approaches.We detail how the model is constructed, compare our results with previous M-R relations, and study how the posterior radius-density model compares going from M-dwarf to solar-type hosts.Finally, we also explore the use of spright to study whether the current observations support the existence of water worlds as a separate population between rocky planets and sub-Neptunes, as suggested by Luque & Pallé (2022).

Overview
The spright package provides a numerical radius-density-mass relation (referred to as "numerical model" from now on) for small planets.The numerical model is constructed by averaging an analytical three-component mixture probability model ("analytical model", from now on) over its posterior parameter space given a catalogue of empirical planet mass and radius measurements.The analytical model is composed of three generalised Student's t-distributions with the distribution means and weights varying as functions of the planetary radius.The three components correspond to rocky planets, water-rich planets (water worlds), and hydrogen-rich sub-Neptunes, and the model parameterisation is designed so that the water world component is optional.That is, the analytical model can represent the observed mass-radius distribution either as a mixture of rocky planets and sub-Neptunes, or as a mixture of rocky planets, water-rich planets, and sub-Neptunes.Thus, the final numerical average model is agnostic to whether water worlds exist as a distinct population.

Probability distribution
The analytical radius-density relation implemented in spright models the probability distribution for a planet's bulk density, , given the planet's radius, , and model parameter vector,   , as a mixture of three generalised (scaled and transformed) Student's t-distributions with five degrees of freedom,3 as The distributions represent rocky planets ( r ), water-rich planets ( w ), and hydrogen-rich sub-Neptunes ( p );  r ,  w , and  p are the mixture weights with  r +  w +  p = 1 for all ; and    r ,    w , and    p are the distribution-specific parameter vectors.More precisely, the distributions are where  are the distribution means,  are the distribution scale parameters,  are the degrees of freedom, and  r and  w are the mean functions for rocky and water-rich planets, respectively.The package supports two theoretical radius-density models to represent the distribution means: the rocky-planet distribution mean follows the models by Zeng et al. (Z19, 2019) 4 and is parameterised by the iron-rock mixing ratio, ; the water-rich planet distribution mean follows the water-rich models either by Aguichine et al. (A21, 2021) 5 or Zeng et al. (2019) and is parameterised by the H 2 O-rock mixing ratio, ; and the sub-Neptune distribution mean is modelled as a power law with density at 2 ⊕ defined by  and exponent by .The choice to model the sub-Neptune population mean as a power law is motivated by the discussion in Lopez & Fortney (2014).

Mixture weights
The distribution weights vary as a function of the planetary radius, as shown in Fig. 2 The black lines show the models for rocky planets with iron-rock mixing ratio varying from 0% to 100%, and the blue lines show the water-rich planet models with H 2 O-rock mixing ratio varying from 5% to 100%.The analytical radius-density probability model creates the rocky and water-rich planet distribution mean functions by interpolating inside the theoretical models.
end radii,  1 and  2 , and water-sub-Neptune transition start and end radii,  3 and  4 .More precisely,  1 is the radius limit below which all planets are rocky,  2 is the maximum radius for a rocky planet,  3 is the minimum radius for a puffy sub-Neptune, and  4 is the radius limit above which all planets are sub-Neptunes.
The mixture weights are calculated by first mapping the planet's radius to a 2D triangle defined by the three composition classes, with rocky planets located at (0,0), water-rich planets at (1,0), and sub-Neptunes at (0,1).The (, ) coordinates for any  are where  denotes clipping the the value  between 0 and 1, that is,  = max(0, min(, 1)).Next, the (, ) coordinates are mapped to the mixture weights through linear interpolation inside the composi- where the weight calculation is simplified from the general case due to the choice of the vertex locations.

Parameterisation
The full parameterisation of the analytical model is shown in Table 1.
As mentioned, the rocky-planet transition start,  1 , stands for the radius below which all planets are rocky, while the sub-Neptune transition end,  4 , stands for the radius above which all the planets are puffy Neptune-like planets.The water world population strength and shape parameters,  and , define the water-rich planet population and are mapped to  2 and  3 as where  = 0.5 − | − 0.5|.As shown in Fig. 2, the mapping is chosen so that the model can represent scenarios where water-rich planets do not form a separate composition class of their own.For  = 0,  1 =  3 and  2 =  4 , and the rocky planet distribution transitions to the sub-Neptune distribution without a water world population in between; for  = 0.5, the water world population weight reaches unity at a single point between  1 and  4 ; and for  = 1,  1 =  2 and  3 =  4 , and the water world population weight is unity for all radii between  1 and  4 . 7he rocky planet iron ratio and water-rich planet water ratio,  and , respectively, map directly to the iron and water mass fractions in the Aguichine et al. (2021) and Zeng et al. (2019) models, and are used to interpolate the radius-density curves from the models.The parameters  and  define the location and exponent of the power-law mean function of the sub-Neptune distribution, respectively.Finally, the three pdf scale parameters define the log-scales of the Student's t-distributions.The model allows for three distinct small-planet populations: rocky planets (brown), water worlds (blue), and hydrogen-rich sub-Neptunes (orange).The density probability is a mixture of three generalised Student's t-distributions with five degrees of freedom, mean following theoretical radius-density models by Aguichine et al. (2021) or Zeng et al. (2019), and width being a free parameter in the fit.The mixture weights ( r ,  w , and  p ) are calculated based on two transition regimes defined by ( 1 and  2 ) and ( 3 and  4 ).In each of the three panels, the upper sub-panel shows the mean and one-sigma limits for the density distributions for each mixture component: the solid line shows the radius regime where the component explains the density fully (without the need for other components), and the dashed line shows an overlapping region between components.The lower sub-panels show the weights for the mixture components.The model can explain the observed radius-density distribution with or without the water world population.Panel a) shows a model realisation with  = 0 (that is,  1 =  3 and  2 =  4 using Eqs. 10 and 11), which leads to a direct transition from rocky planets to sub-Neptunes and excludes water worlds.Panel b) shows a realisation with  = 1 ( 1 =  2 and  3 =  4 ), where water worlds are explained as a well-defined separate population between rocky planets and sub-Neptunes.Finally, panel c) shows a realisation with  = 0.25 corresponding to a weak water world component.Further details about the component weights can be found from a Jupyter notebook at https://github.com/hpparvi/spright/blob/main/notebooks/A1_analytical_model_weights.ipynb.

Motivation
Figure 2 shows three individual realisations of the analytical mixture model that all agree with the observed radius-density distribution for small planets around M dwarfs within the observational uncertainties.Instead of choosing to use the best-fit model, we opt for a more robust approach and calculate a numerical radius-density probability model that is constructed by averaging the analytical model over its posterior parameter space.That is, the numerical model corresponds to the mean posterior predictive distribution (Gelman et al. 2013) of the analytical model given a set of planetary radius and mass observations with their uncertainties.This approach has several advantages.On the one hand, it avoids both the rigidity arising from representing the M-R relationship as a simple parametric model and the need for large sample sizes necessary to make non-parametric models reliable.On the other hand, it lets the underlying analytical model take advantage of the theoretical radius-density models for rocky and water-rich planets by Zeng et al. (2019) and Aguichine et al. (2021) but still ensures that the numerical model is not critically sensitive to the mean functions.The choice of these mean functions offers additional physical interpretability since the model can be parameterised by iron-rock and water-rock ratios, but they could be replaced by power laws -or a different set of planetary internal structure models (e.g., Dorn et al. 2015;Mousis et al. 2020) -with relatively minor effects on the posterior model.Further, the approach allows us to also average over different mean density models to increase the robustness of our prediction.8

Likelihood
Ignoring observational uncertainties, our analytical radius-density probability model would give the log likelihood directly as a sum of the log probabilities of the  observed (, ) points as However, the planet radius estimates from transit observations and mass estimates from RV observations have significant uncertainties that must be considered in the likelihood model.The code takes the uncertainties into account by drawing  sets of radius and mass samples for each planet from the probability distributions defined by the observational uncertainties and transforming the mass samples into densities, as shown in Fig. 3.After this, the code adopts the likelihood averaged over the  samples for each observation as the model log-likelihood, For now, the mass and radius estimate uncertainties are assumed to be normal and symmetric, but we are planning to allow the observations to be represented by freely chosen probability distributions.This will allow, for example, the use of observations with only an upper limit on the planetary mass.

Priors
We list the priors for the analytical model parameters in Table 1.The parameters controlling the transitions and the rocky planet iron-rock ratio have uninformative (uniform) priors.For Zeng et al. (2019), the rocky planets with low Fe content are degenerate with water-rich planets with low H 2 O content.To ensure that the water-world population actually represents water-rich planets, we set a normal prior centred at 0.5 with a width of 0.1 for the H 2 O-rock mixing ratio when using the Zeng et al. (2019)  The parameters defining the transition regions of the model are given additional constraints to ensure that  1 ≤  4 and that the sub-Neptune density at the beginning of the sub-Neptune population ( 3 ) is never larger than rocky-planet density at the same radius.

Creation of the numerical model
as illustrated in Fig. 4.After computing the numerical radius-density probability table, the code computes a discretised cumulative distribution function (CDF) for the planetary bulk density as a function of planet radius and, from this, a discretised inverse cumulative distribution function (ICDF) as a function of planet radius, as shown in Fig. 5.The posterior probability table, ICDF, parameter posterior samples, and the observational data are all then saved to a fits file used by the density and mass predicting part of the code.
The creation of the numerical model is relatively fast, even for large radius and mass data sets, and scales linearly with the number of samples.The evaluation of the log-likelihood function (Eq.13) is parallelised to take advantage of modern multi-core processors, so that, for example, calculating a new model for a data set with 158 planets (i.e. the TEPCat FGK catalogue discussed later) takes 3-7 minutes on a relatively modern eight-core desktop computer.

Evaluation of the numerical model
spright uses inverse transform sampling (Fig. 6) to draw a planet density sample given a planet radius with its uncertainties and a saved numerical radius-density model.The code draws a number of (, ) samples where  follows from the planet radius distribution 9 The number of total and warm-up iterations are defined by the user when calculating a new model.The models included with the package used 60000 iterations in total with a warm-up period of 50000 iterations, where the quality of the final samples was confirmed by inspecting the evolution of the chain population.
and  ∼  (0, 1) and obtains a density sample for each (, ) value by linearly interpolating the 2D ICDF table.After this, a mass sample is obtained from the density sample by multiplying the densities with the respective planet volumes.
The model evaluation is extremely fast since it consists only of the generation of  samples from the two distributions followed by bilinear interpolation inside the ICDF table.

Model creation
The creation of a new radius-density-mass relationship is carried out with the spright.RMEstimator class.At its simplest, the class can be initialised with the system names, a list of planetary radii with their uncertainties, and a list of planetary masses with their uncertainties.The initialisation is followed by model optimisation, parameter posterior estimation, and ICDF map computation: from spright import RMEstimator rme = RMEstimator ( names =names , radii =radii , masses = masses ) rme. optimize () rme.sample () rme.compute_maps () rme.save('map_name .fits') After the ICDF map is computed, it can be saved to be used in model evaluation.The class also allows the model creation to be tuned for specific interests by changing parameter priors or setting additional constraints.

Model evaluation
After a radius-density model has been computed, it can be evaluated using the spright.RMRelation class.The class offers methods to predict the planet's bulk density, mass, and RV semi-amplitude10 distributions given its radius with uncertainties, or the planet's radius given its mass and its uncertainty.For example, a mass distribution for an  = 1.8 ± 0.05 ⊕ planet can be obtained as from spright import RMRelation rmr = RMRelation ('map_name .fits') mds = rmr.sample ('mass ', (1.8 , 0.05)) mds.plot () where mds is a spright.Distribution object that offers utility methods for plotting the distribution, approximating it with analytical (mixture) distributions, calculating distribution percentiles, and so on.Figure 8 shows the plot created by the mds.plot() method visualising the actual numerical mass distribution, an analytical distribution fitted to the numerical distribution, and a set of distribution percentile limits.spright 7

Model evaluation from the command line
The package includes a command line tool spright that makes the model evaluation easy without any coding required.The example above can be evaluated from the command line as spright --predict mass --radius 1.8 0.05 where the script can also save a publication-quality plot of the predicted distribution and the radius-density map used to create the prediction.

Comparison between different catalogues
The flexibility of spright allows the user to quickly generate a new joint radius-density probability distribution based on different empirical data sets.In this section, we use spright to compare the distributions obtained using two catalogues of exoplanet properties and two sets of theoretical density models for the water-world population.The water-rich planet density models are the previouslymentioned models by Aguichine et al. (2021)  b) a sample of small planets orbiting FGK stars from the TEPCat catalogue (Southworth 2011) with better than 25% and 8% uncertainties in the mass and radius, respectively (containing 159 planets spanning a range of 380-2350 K in equilibrium temperature and 0.59 to 1.26  ⊙ in stellar host mass). 12  Figure 9 shows the radius-density and radius-mass distributions inferred from the two catalogues and water-world density models, and a comparison of the model parameter posterior distributions can be found in Figs.A1, A2, A3, A5, A4, A6, A7, A9, and A8.The values of most of the model parameters are consistent regardless of the catalogue used.However, the uncertainties of the model parameters decrease for the TEPCat FGK catalogue (particularly those related to the water world and sub-Neptune transition), which is the one with the largest number of planets.This result highlights the importance of increasing not only the precision and accuracy of exoplanet observed properties but also the number of planets with those properties measured in such a manner.
On the one hand, the iron-rock mixing ratios and the parameters defining the sub-Neptune population agree between the two waterrich planet density models, but differ between the STPM and TEPCat FGK catalogues.On the other hand, the water-rock mixing ratios differ between the density models but agree between the catalogues.The discrepancy in the water-rock mixing ratio between the A21 and Z19 models is however to be expected because the parameter has a different physical meaning for each model (see Zeng et al. 2019;   11 Including planets with measured radii and masses published between June 2021 and 2023.Available in the spright GitHub repository. 12The spright package also includes a model calculated for a sample of small planets ( < 4  ⊕ ) around M dwarfs ( eff < 4000 K) taken from the TEPCat catalogue (containing 54 planets spanning a range of 168-1323 K in equilibrium temperature and 0.089 to 0.65  ⊙ in stellar host mass), but the inferred parameter posteriors are so similar to the ones inferred from the STPM catalogue that we do not include a comparison here.Mousis et al. 2020;Aguichine et al. 2021;Dorn & Lichtenberg 2021, for details).
The water world population seems to shift to higher radii for the FGK hosts compared to the M dwarfs.Type I migration models, independent of the solid accretion mechanism (planetesimal-or pebblebased), predict a planet-to-disk-mass dependent mass-scale where planets migrate (Burn et al. 2021;Schlecker et al. 2022).Thus, in the M-dwarf case, lower-mass water worlds are able to migrate rapidly enough to reach the distances probed in the catalogues compared to the solar-type hosts.This dichotomy in migration timescales could be responsible for the shift to larger sizes of the water world population observed in the FGK versus M dwarf sample.
Compared to the STPM, the TEPCat FGK catalogue shows a stronger separate water-rich population from rocky worlds.The transition is also sharper, indicating a smaller overlap between rocky and water worlds.As discussed above, this effect can also be understood as a consequence of the higher minimum-mass water-rich planets able to migrate to the inner parts of the system for FGK hosts compared with M dwarfs.The lack of low-mass water-rich planets in the FGK sample, however, can be related to an observational bias.Water-rich planets around FGK stars with radii below 2  ⊕ generally have RV semi-amplitudes of 1 m/s or smaller,13 which limits their mass estimation with current RV instrumentation.Observational and detectability biases manifest also in the FGK sample at the Earth-and sub-Earth-size limit, where the lack of rocky planets with measured bulk densities is not because they are intrinsically rare (as clearly demonstrated by, e.g., Batalha et al. 2013), but due to their hardly detectable sub-meter-per-second RV signals.

Comparison with previous mass-radius relations
A significant amount of work has been invested during the last years to model the relationship between planetary radii and masses (Lissauer et al. 2011;Weiss & Marcy 2014;Wolfgang et al. 2016;Mills & Mazeh 2017;Chen & Kipping 2017;Bashi et al. 2017;Ning et al. 2018;Kanodia et al. 2019;Otegi et al. 2020).Most of these works have followed a parametric approach, modelling the mass-radius dependence as one or multiple power-law segments (e.g., Weiss & Marcy 2014;Wolfgang et al. 2016;Chen & Kipping 2017;Otegi et al. 2020); although non-parametric approaches have also been explored (e.g., Ning et al. 2018;Kanodia et al. 2019).The numerical M-R relation offered by spright aims to overcome the shortcomings of the existing parametric and non-parametric models.It offers the flexibility and robustness of non-parametric approaches (in particular, modelling the joint radius-density distribution rather than a single variable) while being conceptually and computationally simple to interpret and incorporate into a Bayesian framework.Figure 10 shows the radius-density and radius-mass relations obtained with spright compared to others that are widely used in the community (e.g., Weiss & Marcy 2014;Chen & Kipping 2017;Otegi et al. 2020).
Among these works, the results obtained with spright are particularly consistent with those by Otegi et al. (2020).Otegi et al. (2020) introduced two separate power law components based on the division in mass-radius space set by the equation of state of water (Dorn et al. 2015), allowing a better representation of the transition region between rocky and volatile-rich planets.The agreement between their volatile-rich power law segment and our hydrogen-rich model is remarkable.But, their rocky component seems to overestimate the size of the rocky planets with masses between 5 to 10  ⊕ .Furthermore, contrary to spright, the M-R relations from Otegi et al. (2020) require prior knowledge of the planet's bulk density to choose the adequate power-law, which is key to adequately predict the mass of the planets with radii between 1.5 and 3  ⊕ .
Regarding the other models, spright obtains consistent results for the rocky population compared to Weiss & Marcy (2014) and Chen & Kipping (2017) up to approximately 1.5  ⊕ .However, Chen & Kipping (2017) sets the transition between rocky and volatilerich planets at 2  ⊕ , thus failing to reproduce the high-mass tail of the rocky population with masses between 2 and 10  ⊕ that overlaps with the water worlds and is particularly prominent in the FGK catalogue.This difference is partially due to our use of new mass and radius data that was not available when the Forecaster model was fit, and partially due to Forecaster using a piecewise model to represent the mass-radius relation, while spright is using a mixture model.That is, Forecaster cannot model a situation with two overlapping populations, but spright assumes by default that the populations can overlap.Outside of the low-mass range, the agreement is very good.In particular, for the hydrogen-rich sub-Neptune population, spright infers a power law index of −0.8 ± 0.3 while Chen & Kipping (2017) obtained −0.77 ± 0.13.Conversely, Weiss & Marcy (2014) fails to predict the properties of the volatilerich population (overestimating the average size and, consequently, underestimating the density of this type of planet).

Evidence for water worlds
The analytical mixture model can explain the observed radius and mass distribution either as a mixture of rocky planets and sub-Neptunes, or as a mixture of rocky planets, water-rich planets, and sub-Neptunes (see Fig. 2).This flexibility allows the model to be used for studying whether observations support or contradict the presence of water-rich worlds as a distinct planetary population positioned between rocky planets and sub-Neptunes.
The analytical model uses the water world population strength, , to parameterise the significance of the water world population:  = 0 corresponds to a model without water worlds,  = 0.5 corresponds to a model where the water world population weight reaches unity for a single point in the planet radius space, and 0.5 <  ≤ 1 correspond to models where the water world population weight is unity for a fraction of the transition region between rocky planets and sub-Neptunes.This parametrisation allows us to present three competing hypotheses: H 0 ) water worlds do not exist as a distinct population, H 1 ) water worlds exist as a mixed population, or H 2 ) water worlds exist a significant distinct population, merely by setting different priors on .For H 0 , we can force the model to exclude the water world population by setting a delta function prior, (0) on .The mixed-population hypothesis, H 1 , represents a scenario where a population of water worlds exists mixed with rocky planets and sub-Neptunes, but there are no planetary radii for which all the planets would be purely water worlds.We encode H 1 by a uniform prior on  from 0.015 to 0.5, where the lower bound represents a weak water world population and the upper bound represents the case where the water world population weight reaches unity for a single point in the planet radius space.Finally, for H 2 , we choose to encode a "significant distinct population" by a uniform prior from 0.5 to 1.0.That is, the water world population weight must reach unity for at least one point in the radius space, and the upper Catalogue 2 ln  10 2 ln  20 2 ln  21 STPM Z19 -0.6 ± 0.7 -3.4 ± 0.7 -2.8 ± 0.7 STPM A21 -7.5 ± 0.6 -8.5 ± 0.7 -1.0 ± 0.6 TEPCat FGK Z19 2.5 ± 0.7 2.6 ± 0.7 0.2 ± 0.7 TEPCat FGK A21 -1.0 ± 0.7 -1.5 ± 0.7 -0.5 ± 0.7 Table 2. Bayes factors inferred from the STPM and TEPCat FGK catalogues for two theoretical water-rich-planet mean radius-density models and three hypotheses detailed in Sect.3.3.The Z19 scenarios use the Zeng et al. (2019) models to represent the mean radius-density function of the water-rich-planet population, while the A21 scenarios use the models by Aguichine et al. (2021).
limit marks the model with sharp transitions from rocky planets to water worlds and from water worlds to sub-Neptunes.After defining our hypotheses, H 0 , H 1 and H 2 , we can calculate their respective Bayesian evidences, Z 0 , Z 1 and Z 2 , by integrating over the model posterior spaces (Parviainen 2018;Gelman et al. 2013;Robert 2007;Kass & Raftery 1995). 14We carry out the integration using the Dynesty package (Koposov et al. 2023;Speagle 2020;Feroz et al. 2009;Skilling 2006Skilling , 2004) ) to estimate the Bayesian evidences via dynamical nested sampling.
We report the log Bayes factors -defined as 2 ln  10 = 2 ln  1 − 2 ln  0 and 2 ln  20 = 2 ln  2 − 2 ln  0 -in Table 2. Assuming the Zeng et al. (2019) water-rich planet density models, the evidence is insufficient to significantly support any of the scenarios over the others (Kass & Raftery 1995, p. 777). 15Assuming the Aguichine et al. ( 2021) models, H 0 (no water world population) is strongly favoured over H 1 and H 2 for the STPM catalogue, but only tentatively favoured for the TEPCat FGK catalogue.The reason for this discrepancy between Z19 and A21 models for the STPM catalogue comes from the impossibility of A21 models to match the slope of the waterrich planet population in mass-radius space for low-mass planets (< 5  ⊕ ), while it is easily reproduced by Z19 models.On the other hand, the higher average planet radius of the FGK sample is easier to reproduce with the A21 models compared to Z19.We note that spright only uses a set of A21 models that assume  irr = 500 K and an Earth-like composition for the core.A thorough exploration of our model comparison results as a function of these two parameters is beyond the scope of this paper.

Parameter posteriors
While the main use case for spright is in predicting planet masses given their radii and vice versa, the model parameter posteriors can also be used to study the physical properties of small-planet populations.For this use, however, we need to understand how the posteriors depend on factors such as the number of planets included in the model calculation.
To this end, we carried out a synthetic catalogue study for all combinations of four catalogue sizes (the number of planets included in the catalogue,  p ∈ {50, 100, 150, 200}) and three values for the water world population strength ( ∈ {0.0, 0.5, 1.0}).We created five realisations of synthetic (radius, mass) catalogues with realistic uncertainties on both quantities for each ( p , ) combination.The parameters other than  were fixed to  1 = 1.3,  4 = 2.4,  = 0.0,  = 0.2,  = 0.5,  = 3,  = −0.9, r = −0.4, w = −0.3, and  p = −0.4.The synthetic mass and radius data sets were created using the create_mock_sample helper function in spright, and the relative radius uncertainties were drawn from a uniform distribution from 1% and 8%, while the relative mass uncertainties were drawn from a uniform distribution from 3% and 24%.To simplify the analysis, we only use the Zeng et al. (2019) density models for the water-rich planets.
We created the numerical spright model for each synthetic catalogue, and show the inferred parameter posteriors in Fig. 11.The true parameter values are generally contained inside the inferred 68% central posterior limits, and only very rarely outside the 95% central posterior limits.The posterior uncertainties decrease with the increasing catalogue size as expected, except for the water-world population shape parameter, .This is not entirely surprising because  has the strongest effect on the shape of the water-world population for intermediate values of  and it will likely be well-constrained only for relatively strong water-world populations and large catalogue sizes.
We also studied how the posterior estimate for  changes for the three simulated  scenarios and five catalogue sizes, and show the results in Fig. 12.The water world population strength is relatively well constrained for both extreme cases even for  p = 50, and the distribution mode is in all cases close to the true  value.The intermediate  = 0.5 scenario is less well constrained, but the posteriors generally differ from the posteriors from the extreme cases.In practice, this means that the model can distinguish between the two extreme cases for , and a poorly-constrained  can be interpreted as support for the existence of a population of water-rich planets that is mixed with the rocky and sub-Neptune populations.

Water-world evidence
We repeated the Bayesian model comparison test in Sect.3.3 using synthetic catalogues.We calculated the Bayesian evidences for the H 0 , H 1 , and H 2 hypotheses for 20 catalogue realisations for each ( p , ) combination, where  p ∈ {50, 100, 150, 200} and  ∈ {0.0, 0.25, 0.5, 0.75, 1.0}, again restricting the simulations to the Zeng et al. (2019) density model for the water-rich planets. 16 We visualise the resulting evidence distributions in Fig. 13 and summarise them in Table 3.For  = 0, the log Bayes factors, 2 ln  10 and 2 ln  20 , generally support the no-water-world-population hypothesis H 0 over the two others.The evidence against H 1 is rather tentative for small  p , and can be weak even for large  p .The strongwater-world-population scenario, H 2 , can be ruled out in most cases with "positive" evidence already with small  p , and decisively with large  p .For  = 0.25, the evidence for the mixed-water-worldpopulation hypothesis, H 1 , reaches the level of "strong" evidence for some catalogue realisations, but in most cases, neither H 0 nor H 1 is strongly favoured over another.Both H 0 and H 1 are generally favoured over H 2 , but the level of support for H 1 does not increase significantly with increasing  p .For  = 0.5, 2 ln  10 is nearly always positive and reaches high levels of evidence for larger  p , while 16 We will repeat the simulation for the Aguichine et al. (2021) water-rich planet density models in the future and make the results publicly available from the spright GitHub repository, but this work is beyond the scope of this paper.the support for H 2 over H 0 varies from slightly negative to strongly positive.For  = 0.75, both H 1 and H 2 are significantly favoured over H 0 for all  p , and H 2 is generally favoured over H 1 .For  = 1, The strong-water-world-population hypothesis H 2 is favoured over H 0 and H 1 with decisive support already for  p = 50.All in all, the synthetic tests show that the spright model can be used to distinguish between the two extreme cases described by  values of 0 and 1 already with a catalogue consisting of ∼ 50 planets.Interestingly, the scenario with  = 1 can be identified much more securely than the scenario with  = 0.In most cases, 2 ln  20 < −2 for  = 0 and 2 ln  20 > 6 for  = 1, and the contrast between H 0 and H 2 increases quickly together with the number of planets included into the catalogue.The intermediate cases with  values of 0.25 and 0.5 are identified less securely, but for  ∼ 0.5, 2 ln  10 is still generally positive while being negative for  = 0.

CONCLUSIONS
We have presented spright,17 a Python package that implements a lightweight probabilistic radius-density-mass relationship for small planets based on basic Bayesian inference.The package represents the joint planetary radius and bulk density distribution as a mean of the posterior predictive distribution of a simple analytical threecomponent mixture model.The package offers tools to predict planetary masses, bulk densities, and radial velocity semi-amplitudes for planets orbiting M dwarfs (based on the revised STPM catalogue by Luque & Pallé 2022) and FGK stars (based on TEPCat catalogue).The package has been designed to be easy to install and use and also aims to make the computation of new M-R relations easy.Further, calculating a new M-R model takes only minutes with a modern desktop computer, even for large data sets containing hundreds of planets, and the computing time scales linearly with the number of planets.
We have also studied whether the current observational radius and mass estimates support the existence of water-rich worlds as a distinct planet population between rocky planets and sub-Neptunes.While the numerical M-R model is agnostic to what comes to the existence of water worlds, the analytical model can be used in a Bayesian model comparison setting to assess the level of evidence in favour of a distinct water world population.Our study finds that the inferred support for the existence of a water world population depends on the chosen theoretical water-rich planet density model.All in all, the TEPCat data set is insufficient to provide statistically significant evidence for or against the existence of a water world population around FGK stars.The STPM data set shows some evidence against the existence of a water world population around M dwarfs, but this evidence is not strong enough to be considered conclusive, in line with the recent results from Rogers et al. (2023).0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 Np = 200

Figure 1 .
Figure1.Theoretical models for the bulk planet density for rocky and waterrich planets byZeng et al. (Z19, 2019)  andAguichine et al. (A21, 2021).The black lines show the models for rocky planets with iron-rock mixing ratio varying from 0% to 100%, and the blue lines show the water-rich planet models with H 2 O-rock mixing ratio varying from 5% to 100%.The analytical radius-density probability model creates the rocky and water-rich planet distribution mean functions by interpolating inside the theoretical models.

Figure 2 .
Figure 2. Three realisations of the analytical density mixture probability model used to calculate the final numerical radius-density probability model.The empirical dataset used to fit the model is the STPM catalogue used in Section 3.1.The model allows for three distinct small-planet populations: rocky planets (brown), water worlds (blue), and hydrogen-rich sub-Neptunes (orange).The density probability is a mixture of three generalised Student's t-distributions with five degrees of freedom, mean following theoretical radius-density models byAguichine et al. (2021) orZeng et al. (2019), and width being a free parameter in the fit.The mixture weights ( r ,  w , and  p ) are calculated based on two transition regimes defined by ( 1 and  2 ) and ( 3 and  4 ).In each of the three panels, the upper sub-panel shows the mean and one-sigma limits for the density distributions for each mixture component: the solid line shows the radius regime where the component explains the density fully (without the need for other components), and the dashed line shows an overlapping region between components.The lower sub-panels show the weights for the mixture components.The model can explain the observed radius-density distribution with or without the water world population.Panel a) shows a model realisation with  = 0 (that is,  1 =  3 and  2 =  4 using Eqs. 10 and 11), which leads to a direct transition from rocky planets to sub-Neptunes and excludes water worlds.Panel b) shows a realisation with  = 1 ( 1 =  2 and  3 =  4 ), where water worlds are explained as a well-defined separate population between rocky planets and sub-Neptunes.Finally, panel c) shows a realisation with  = 0.25 corresponding to a weak water world component.Further details about the component weights can be found from a Jupyter notebook at https://github.com/hpparvi/spright/blob/main/notebooks/A1_analytical_model_weights.ipynb.

Figure 3 .
Figure 3. Generation of the mass and density samples for the likelihood evaluation.Panel a) depicts three planetary radius and mass measurements with their uncertainties, panel b) shows a set of radius and mass samples generated for each planet, and panel c) shows the radius and density samples derived from the radius and mass samples that are used in the model likelihood evaluation.
models.The meaning of the H 2 O-rock mixing ratio is different between the Zeng et al. (2019) and Aguichine et al. (2021) models (condensed versus supercritical water), and the latter do not have a problem with degeneracy with the Zeng et al. (2019) rocky-planet models.Because of this, we lift the normal prior constraint used with the Zeng et al. (2019) water-rich planet models and use a wide uniform prior on the H 2 O-rock mixing ratio when using the Aguichine et al. (2021) models.The sub-Neptune density normalisation factor and exponent have uninformative priors, and the logarithms of the Student's t-distribution widths' logarithms have loosely constraining normal priors.
When building the final numerical radius-density model, spright first finds the global mode of the posterior density given the observational radius and mass estimates, log (  |, ) = log  + log (  ), (14) where log  is the log-likelihood from Eq. 13 and log (  ) is the logprior.The optimisation is carried out using the Differential Evolution global optimisation method (DE, Price et al. 2005) implemented in PyTransit (Parviainen 2015), with the initial parameter vector population drawn from the model parameter prior distribution.After the optimisation, the code obtains a sample from the model parameter posterior using the emcee Markov chain Monte Carlo sampler (Foreman-Mackey et al. 2013).The emcee sampler is initialised using the parameter population clumped around the global posterior mode by the DE method, and the sampler is run until it has obtained a representative sample from the posterior distribution. 9Next, the code discretises the radius-density space into a two-dimensional array and computes the posterior probability for each discrete (, ) point by averaging the analytical probability model over the posterior samples, , |    ), and Zeng et al. (2019), and the two catalogues are: a) an updated 11 version of the Small Transiting Planets around M dwarfs (STPM) catalogue by Luque & Pallé (2022) (containing 48 planets spanning a range of 168-1089 K in equilibrium temperature and 0.089 to 0.63  ⊙ in stellar host mass),

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .
Figure 4. Construction process of the numerical radius-density probability model.The panels exhibit the mean values obtained from  s samples of the analytical mixture model, drawn from its posterior distribution.In the visualisation, the rocky planet population is represented by brown colour, the water world population by blue, and the sub-Neptune population by orange.It is important to note that the figure presents the averages of a single isocontour for each component for visual clarity, whereas the actual model considers averages over real-valued probabilities.The default models included with the spright have been averaged over 3000 posterior samples.

Figure 9 .
Figure 9. Numerical radius-density (left) and radius-mass (right) probability models fitted to the STPM M dwarf catalogue and the TEPCat FKG star catalogues using either Zeng et al. (2019, Z19, ) or Aguichine et al. (A21, 2021) water-rich planet density models to represent the density mean function for the water worlds, as explained in Sect.3.1.Grey data points show radius, density, and mass measurements with their uncertainties for all planets in each catalogue.The blue colour corresponds to the logarithm of the posterior probability, and the black lines show the posterior means for each of the three planet populations: the solid lines correspond to radius regimes where the component has a weight of unity (that is, all planets in this range belong to this component), while the dashed lines mark the transition regimes between the populations.

Figure 10 .
Figure 10.The spright radius-density and radius-mass relations for the STPM and TEPCat catalogues and Aguichine et al. (A21, 2021) and Zeng et al. (Z19, 2019) models for the water-rich planet densities plotted together with the ones by Otegi et al. (2020), Chen & Kipping (2017), and Weiss & Marcy (2014).The solid lines correspond to the spright posterior distribution means for planet radii where the component weight is larger than 0.1, the slashed lines to the Otegi et al. (2020) means, the dotted lines to the Chen & Kipping (2017) means, and the dash-dotted lines to the Weiss & Marcy (2014) relation.The blue shading shows the spright joint probability densities.We do not show the model uncertainties for visual clarity.

Figure 11 .
Figure 11.Posterior distributions for all the spright model parameters inferred from the 75 synthetic mass and radius catalogues described in Sect.3.4.1.The light and dark vertical bars show the central 95% and 68% posterior limits, respectively, and the horizontal dotted line shows the true value for all parameters except the water world population strength, .The simulations were carried out for five catalogue sizes ( p ) and three values of , and the posteriors are grouped first by  p (outer level grouping with separate colour for each  p ), then by  (a set of five posterior estimates), and finally by the catalogue realisation (a single vertical line).

Figure A1 .
Figure A1.Posterior distributions for the analytical mixture model parameters inferred from the updated STPM catalogue and TEPCat FGK star host catalogue, and quantities derived from the model parameters.RP refers to rocky planets, WP to water-rich planets, and SN to hydrogen-rich sub-Neptunes, while Z19 refers to theZeng et al. (2019) water-rich planet density models, and A21 to theAguichine et al. (2021) water-rich planet density models.The transition width is the width of the transition region between rocky planets and sub-Neptunes ( 4 −  1 ), WP population centre is the centre-of-mass radius for the water world population calculated from the water world population weights over the transition region, WP average weight is the average water world population weight over the transition region, and WP integrated weight is the weight of the water world population integrated over all planetary radii.

Figure A4 .
Figure A4.Joint posterior distributions for a set of analytical mixture model parameters and quantities derived from the model parameters inferred from the updated STPM catalogue using the Zeng et al. (2019) water-rich planet density models.Here,  4 −  1 is the width of the transition region between rocky planets and sub-Neptunes, ∫  w ( )d /( 4 −  1 ) is the mean water world population weight over the transition region, ∫   w d / ∫  w ( )d is the water world population centre calculated as a weighted mean of the planet radius, and ∫  w ( )d is the total integrated water world population weight.

Figure A8 .
Figure A8.Joint posterior distributions for a set of analytical mixture model parameters and quantities derived from the model parameters inferred from the updated STPM catalogue using the Aguichine et al. (2021) water-rich planet density models.Here,  4 −  1 is the width the transition region between rocky planets and sub-Neptunes, ∫  w ( )d /( 4 −  1 ) is the mean water world population weight over the transition region, ∫   w d / ∫  w ( )d is the water world population centre calculated as a weighted mean of the planet radius, and ∫  w ( )d is the total integrated water world population weight.

Figure A9 .
Figure A9.As in Fig. A6 but for the TEPCat FGK host star sample.

Table 1 .
Analytical mixture model parameters and priors.N( , ) stands for a normal prior with a mean  and standard deviation , and U(, ) stands for a uniform distribution from  to .

Table 3 .
Minimum and maximum Bayes factors (2 ln ) for hypotheses H1 and H2 in favour of H0 from the synthetic catalogue simulations visualised in Fig.13.
Figure 12.Posterior distributions for the water world population strength  inferred from the 75 synthetic mass and radius catalogues described in Sect.3.4.1.The columns show the posteriors for a given  t (where the subscript stands for "truth"), and the rows for a different catalogue size,  p .Figure 13.Log Bayes factors from the synthetic catalogue study described in Sect 3.4.2.The columns show the Bayes factors for hypotheses H 1 and H 2 in favour of H 0 , and H 2 in favour of H 1 , and the rows show the different  scenarios.The number of planets in the catalogue is shown on the x-axis, and the y-axis shows the 68% central intervals of the 2 ln  value distributions as blue boxes while black vertical lines mark the distribution minimum-tomaximum spans.The yellow shading follows the four levels of support by Kass & Raftery (1995): 2 ln  < 2) insignificant support, 2 < 2 ln  < 6) positive support, 6 < 2 ln  < 10) strong support, and 2 ln  > 10) very strong support.Negative ln  values mean the alternative hypothesis is favoured, so that 2 ln  10 < −10 would mean very strong support for H 0 over H 1 .
Zeng et al. (2019)stributions for a subset of the analytical mixture model parameters inferred from the updated STPM catalogue using theZeng et al. (2019)water-rich planet density models.Joint posterior distributions for a subset of the analytical mixture model parameters inferred from the TEPCat FGK catalogue using theZeng et al. (2019)water-rich planet density models.