© The Authors 2024

1 Introduction

Since the detection of the first exoplanet around a Sun-like star in 1995 (Mayor & Queloz 1995) over 5000 exoplanets have been discovered revealing a large diversity in their physical properties. The field of exoplanets is blossoming. We are at a stage where we move from exoplanet detection to exoplanet characterization, with work on both of these fronts being extremely active and reaching a new level. For exoplanet characterization, the two fundamental parameters are the planetary mass and radius. However, these two properties cannot be measured by the same method, and unfortunately in many cases only one of the two is available. To allow for a broader overview of exoplanets in a statistical sense, it is valuable to determine their mass-radius (M–R) relation and how this relation changes for different planetary types. Relating a measured radius to a planet’s mass can also help with radial-velocity follow-ups by estimating the expected radial velocity semi-amplitude of a transiting exoplanet.

The M–R relation depends on the planetary composition and therefore on the behavior of the different materials at planetary conditions (e.g., Seager et al. 2007; Chabrier et al. 2009; Grasset et al. 2009; Mordasini et al. 2012; Spiegel et al. 2014; Jontof-Hutter 2019). Theoretical M–R relations can be inferred from interior models that rely on equations of state (EOSs), which relate the density and pressure (and usually also the temperature) of a given composition. For simplicity, small terrestrial planets are often assumed to have constant densities and are thus expected to behave as R ∝ M^1/3 (e.g., Spiegel et al. 2014). On the other hand, in massive giant planets that are hydrogen-helium (H-He) dominated in composition, the gravitational pressure is high enough for the materials to be pressure ionized, and the electron degeneracy pressure becomes substantial (e.g., Helled et al. 2020). This causes the radius to decrease with increasing mass, with R ∝ M^−1/3. Of course, the planetary radius does not only depend on its mass, but also on other factors such as stellar irradiation or the planetary age.

Today, there are enough planets with mass and radius measurements to statistically infer an observed M–R relation. The observed M–R relations of exoplanets have been investigated in several studies and are often fit by a broken power law. The breakpoints are particularly interesting as they represent the transitions between different planetary types. For example, Weiss et al. (2013) inferred a mass-radius-incident-flux relation. Based solely on the visual inspection of the M–R and mass-density (M–ρ) diagrams, they found a transition point at 150 M_θ. They inferred relations R ∝ M^0.53F^−0.03 for M < 150 M_⊕ and R ∝ M^−0.04F^0.09 for M > 150 M_⊕, where F is the instellation flux. They suggested that for small planets the mass is the most important parameter for predicting the planetary radius, whereas for giant planets the incident flux is more important. Below 150 M_⊕ only 35 planets were available at the time, and the result for this region is less robust.

A different approach was taken by Hatzes & Rauer (2015). In this study, the transition for giant planets was explored. They used observations to infer the mass–density (M–ρ) relation. A minimum in density at ≃95 M_⊕ (0.3 M_J) and a maximum at ≃1.9 × 10⁴ M_⊕ (60 M_J) were inferred, suggesting that these values correspond to the transition into gas giant planets. This study also suggested there is no separation between brown dwarfs and giant planets as they display similar behaviors. The mass-density relation of giant planets was determined to be ρ = 0.78M^1.15.

Chen & Kipping (2017) presented an elaborate Markov chain Monte Carlo method to analyze the M–R relations of objects ranging from dwarf planets to stars. They inferred four distinct regions in their work: Terran worlds, Neptunian worlds, Jovian worlds, and stars. The corresponding transition points were fit and identified as 2.04, 132, and 2.66 ×10⁴ M_⊕, where the latter two values correspond to 0.41 M_J and 0.08 M_⊕. The M–R relations were identified as R ∝ M^0.28 for Terran worlds, R ∝ M^0.59 for Neptunian worlds, R ∝ M^−0.04 for Jovian worlds and R ∝ M^0.88 for stars. The results were obtained by a solely data-driven analysis rather than being derived from any physical assumptions. At the time of the study, only a few objects at ~ 1 M_⊕ had been observed, and therefore the transition point between the Terran and Neptunian worlds relied on a small number of planets.

Bashi et al. (2017) fit the M–R relations of two distinct regions using a total least-squares approach, where the transition point was assumed to be unknown. It was found that the transition occurs at ~124 M_⊕ and 12.1 R_⊕. The data were best fit by the relations R ∝ M^0.55 for small planets and R ∝ M^0.01 for large planets. In this study, the transition point was also fit instead of being imposed by some prior assumption.

Otegi et al. (2020) presented an updated catalog of exoplan-ets, for which robust measurements of both radius and mass are available, based on the NASA Exoplanet Archive catalog¹. They focused on finding the transition between rocky planets and those with a substantial gas envelope, and therefore only planets with masses up to 120 M_⊕ were considered. When displaying the planets in the M–R plane, two distinct regions were identified: the rocky and the volatile-rich populations. Because they overlap in mass and radius, the populations were separated by the pure-water composition line to distinguish volatile-rich planets from terrestrial ones. For the two distinct groups, the following (M–R) relations were inferred: R = 1.03 M^0.29 for the rocky population and R = 0.70 M^0.63 for the volatile-rich population.

Edmondson et al. (2023) showed that a discontinuous M–R relation, as well as a temperature dependence for giant planets, results in a good fit to the M–R measurements. Similarly to Otegi et al. (2020) they separate the rocky and the icy planets with a pure-ice EOS, while finding the transition of icy planets to gas giants at 115 M_⊕. For the rocky planets, they fit a relation of R ∝ M^0.34 and for the icy planets it is R ∝ M^0.55. When also considering the equilibrium temperature of the gas giants they find that R ∝ M^0.00T^0.35. This suggests that for the giant planets, the radius only depends on the temperature.

Recently, Mousavi-Sadr et al. (2023) used a machine-learning approach to analyze the exoplanet population. By applying various clustering algorithms, they identified the transition between small and giant planets at masses of 52.48 M_⊕ and sizes of 8.13 R_⊕. For the small planets, the M–R relation was found to be R ∝ M^0.50. They also showed that the radii of giant planets are positively correlated with the stellar mass.

In this work, we investigated the M–R relations of exoplanets and their transitions using the PlanetS catalog². We used a solely statistical approach to determine the breakpoints in the relation used to define the different planetary regimes and determine the distinct dependencies. We also analyzed the mass-density and radius-density (R-ρ) relations and investigated their validity in separating different planetary types.

2 Methods

In this work, we used the data from the PlanetS catalog. An earlier version was presented in Otegi et al. (2020); since then the catalog has been extended with additional discoveries and planets with masses up to 30 M_J. There was also an update on the planetary masses, and some planetary parameters were reanalyzed. The catalog only includes planets with relative measurement uncertainties on the mass and radius smaller than 25% and 8%. Since the updated catalog contains many more planets, it is more reliable and allows for additional analyses. The data we used were downloaded in July 2023 and contain the mass and radius measurements of 688 exoplanets. Figure 1 shows how the planets are distributed in the M–R plane.

Similarly to previous work our approach assumed a power-law relation between two planetary variables: for example, M and R. The first step was to transform the variables into the log-log plane to use a linear regression method. We further assumed that there is an unknown number of break points in the linear relation between the two log-transformed variables, that is, that a piece-wise linear function describes the data.

For a single breakpoint, the two-segmented piece-wise linear relation can be parametrized as (e.g., Muggeo 2003) $$y=c+\alpha \xi +\beta (\xi -\psi )H(\xi -\psi ),$$(1)

where c is the intercept, α the slope of the first segment, β the difference between the slopes of the two segments, ψ the breakpoint, and ξ the independent variable. H(x) is the Heaviside step function, defined as H(x) = 1 if x ≥ 0 and H(x) = 0 otherwise. For n breakpoints, Eq. (1) can be generalized as $$y=c+{\alpha }_1\xi +{\displaystyle \sum _i^n{\beta }_i}\left(\xi -{\psi }_i\right)H\left(\xi -{\psi }_i\right).$$(2)

Because of the Heaviside step function, Eqs. (1) and (2) are nonlinear, and ordinary linear regression methods cannot fit the parameters. Furthermore, the data used in this work has significant uncertainties on the dependent and independent variables, -which need to be considered. Therefore, in order to determine the piece-wise linear regression fit, we used orthogonal distance regression (ODR; as implemented in the Python package SciPy).

We treated the number of breakpoints n as an additional free parameter. To determine the number of breakpoints that yield the best fit to the data, we fit piece-wise linear functions with zero to four breakpoints. After fitting, we compared the models by calculating the Bayesian information criterion (BIC) in the form that is commonly used for linear regression: $$\text{BIC}=n\ln \left(\frac{1}{n}{\displaystyle \sum _i{\left(x_i-{\widehat{x}}_i\right)}^2}\right)+k\ln n,$$(3)

where x_i are the data points, ${\widehat{x}}_i$ the model predictions, and n, k are the number of data points and model parameters.

For the M–R relation, ξ ≡ log(M[M_⊕]) and y ≡ log(R[R_⊕]), for the M–ρ relation ξ ≡ log(M[M_⊕]) and y ≡ log ρ [g/cm³], and for the R-p relation ξ ≡ log(R[R_⊕]) and y ≡ log ρ [g/cm³]. The planetary bulk density was calculated using $\rho =\frac{3M}{4\pi R^3}$.

Since all fitting variables required the data to be transformed, the uncertainties had to be propagated. For transforming the mass-radius measurements into their logarithms, the error propagation is simply σ_ξ = ∂ (log x/∂x) σ_x = σ_x/x ln 10), with x the measured mass or radius. When the measurement uncertainties were asymmetric, we used their arithmetic mean as σ_x.

For vector-valued functions f of potentially covariant parameters, the more general form of the error propagation has to be used (e.g., Tellinghuisen 2001): $${\sigma }_f^2=g^{T}Vg,$$(4)

where g is the gradient of f whose ith element is ∂f/∂x_i, and V is the covariance matrix. In our case, this applied to the uncertainties of the density and of the power-law parameters that were calculated from the piece-wise linear fit.

Fig. 1 Mass-radius data from the PlanetS catalog and the derived relation. Data from PlanetS catalog is displayed in purple. The best-fit mass–radius relation (Eq. (5)) is represented by the solid black line and the shaded region shows the 1σ-uncertainty. The breakpoints are shown as the dashed lines.

3 Results

In this section, we first present our results for the M–R, M–ρ and R–ρ relations in Sects. 3.1–3.3. We then compare our results to previous studies in Sect. 3.4.

3.1 The mass-radius relation

By comparing the BIC of piece-wise linear models with zero to four breakpoints, we determined that two breakpoints provided the best fit to the M–R distribution of the exoplanets from the PlanetS catalog. This led to the following M–R relation: $$R=\{\begin{array}{ll}(1.02\pm 0.03)M^{(0.27\pm 0.04)}\hfill & M<(4.37\pm 0.72)\hfill \\ (0.56\pm 0.03)M^{(0.67\pm 0.05)}\hfill & (4.37\pm 0.72)<M<(127\pm 7)\hfill \\ (18.6\pm 6.7)M^{(-0.06\pm 0.07)}\hfill & M>(127\pm 7),\hfill \end{array}$$(5)

where R and M are in Earth units. The piece-wise linear fit with two breakpoints is shown together with the data in Fig. 1, and the fit parameters (see Eq. (2)) are listed in Table 1.

The M–R fit with two breakpoints is split into three different regimes (segments). These breakpoints correspond to small planets (M < 4.37 M_⊕), intermediate-mass planets (4.37 M_⊕ < M < 127 M_⊕), and giant planets (M > 127 M_⊕). The first breakpoint has a higher relative uncertainty than the second, with a transition mass of M₁ = (4.37 ± 0.72) M_⊕. Between the first and second segments, the gradient changes less than the transition from the second to the third segment, making it harder to identify the breakpoint. The fit is particularly suitable for determining both breakpoints and describing the data in the first and second segments. The description of the data in the third segment (the giant planets) is relatively uncertain. These planets show a large scatter, making it difficult to find a well-fitting gradient.

The first regime of exoplanets corresponds to planets with masses below 4.4 M_⊕ and roughly follows the relation of R ∝ M^0.27. These planets are most likely “rocky worlds” with compositions similar to the Earth’s. If terrestrial planets can be approximated as constant-density homogeneous spheres, they would follow the simple relation R ∝ M^1/3, which is very similar to what we find. The scatter around this relation in the actual data comes from the differentiated structures of planets and the diversity in their bulk densities due to, for instance, rocks-to-metals ratios and the possible existence of lighter elements such as water (e.g., Seager et al. 2007; Weiss et al. 2013; Zeng et al. 2016).

The change in the slope at ~4.4 M_⊕ defines the transition to the intermediate-mass planets, which could also have non-negligible H-He envelopes. Our fit to the data implies that the maximal mass of “rocky” exoplanets and possibly of naked planetary cores is somewhere between 4 and 5 M_⊕. This also implies that the minimum mass to accrete a substantial amount of volatile elements is about 4 to 5 M_⊕. The transition region is about half of the theoretical mass limit of about 10 M_⊕ for rocky exoplanets (e.g., Seager et al. 2007; Fortney et al. 2007; Charbonneau et al. 20091). The intermediate-mass planets between 4.4 M_⊕ and 127 M_⊕ correspond to planets with H-He envelopes that still have a large heavy-element mass fraction. Since the mass range is large, the diversity of the envelope mass fractions varies significantly and can range from very thin atmospheres to rather gaseous envelopes (e.g., Weiss et al. 2013; Hatzes & Rauer 2015; Ulmer-Moll et al. 2019). These planets have the steepest M–R relation following R ∝ M^0.67. An increase in mass results in a significantly larger radius, corresponding mainly to a larger envelope composed of volatile elements. The transition to the gas giants occurs around 127 M_⊕ and is where the planets start to be dominated by the H-He envelope. Interestingly, this transition mass is consistent with the suggested transition mass to giant planets based on recent giant planet formation models (Helled 2023).

As expected from their H-He dominated composition, for the giant planets we find that the radius is nearly independent of mass (R ∝ M^−0.06). For high-mass objects consisting of a degenerate electron gas, we expect a relation of R ∝ M^−1/3. In the giant planets, the gas is not completely degenerate, leading to a slightly compressible gas and a deviation from the expected relation. At the same time, we also observe a large scatter in radius due to stellar irradiation, different planetary ages, and metallicities that can strongly affect the radii of gas giants (e.g., Thorngren et al. 2016; Teske et al. 2019; Müller et al. 2020; Müller & Helled 2023).

3.2 The mass-density relation

In this subsection, we present a fit to the mass-density (M–ρ) relation using the planets from the PlanetS catalog. As in Sect. 3.1, we first determined the best-fit number of breakpoints. Similar to the M–R relation, we found that two breakpoints provided the best fit to the data. The fitting function was therefore given by Eq. (2)), with ξ = log(M[M_⊕]) and y = log(ρ[g cm⁻³]) and n = 2. The piece-wise linear function with two breakpoints yielded the following power-law M–ρ relation: $$\rho =\{\begin{array}{ll}(5.12\pm 0.36)M^{(-0.12\pm 0.11)}\hfill & M<(3.84\pm 0.76)\hfill \\ (16.6\pm 3.3)M^{(-0.76\pm 0.13)}\hfill & (3.84\pm 0.76)<M<(183\pm 7)\hfill \\ \frac{(4.39\pm 2.06)}{{10}^4}{M}^{(1.26\pm 0.19)}\hfill & M>(183\pm 7),\hfill \end{array}$$(6)

where ρ is in g cm⁻³ and M in M_⊕. The inferred M–ρ fit is shown in the top panel of Fig. 2 together with the data. Table 2 lists the values of the fitting parameters. We find that the first breakpoint in the M–ρ relation is similar, albeit slightly lower than that found in the M–R relation. The second breakpoint, the transition to the giant planets, is at a significantly higher mass (183 M_⊕) compared to what we found previously (127 M_⊕).

As discussed above, terrestrial planets can be approximated by a constant density, which is roughly consistent with our findings. The planets in the second and third segments have a large scatter in the M–ρ plane, implying that the planetary composition is rather diverse. Consequently, there is also a larger uncertainty in the power-law of the M–ρ density for the intermediate and giant planets.

To compare the M–ρ to the M–R relation (as derived in Sect. 3.1), we converted the density ρ(M) from Eq. (6) to a radius using $R={\left(\frac{4\pi \rho }{3M}\right)}^{1/3}$. Since the M–ρ and M–R relations were fit with the same data, they should yield similar results, although some deviations are expected. The results are shown in the bottom panel of Fig. 2. Indeed, we find that the M–ρ and M–R relations have very similar behaviors, and the relation derived from the M–ρ fit lies well inside the uncertainty of the fit to the M–R distribution. The fact that the power laws were similar while being found separately suggests that our approach yields consistent results. As noted earlier, the largest difference between the two fits is the transition mass between the intermediate and giant planets. Using the M–ρ relation yields a significantly higher transition mass.

Fig. 2 Mass–density and mass–radius data from the PlanetS catalog and the derived relations. Top: mass and density data of exoplanets from PlanetS catalog. The mass–density relation is represented by a solid black line. The dashed lines show the position of the breakpoints. The light shaded areas are the corresponding 1σ uncertainties. Bottom: mass versus radius distribution with best fit for relation and breakpoints in black, as in Fig. 1. The orange solid and dashed lines show the mass-radius relation and the corresponding breakpoints derived from the fit to the mass-density distribution.

3.3 The radius-density relation

When fitting the M–R and the M–ρ relations, the three different regimes were defined by a transition mass. Alternatively, is also possible to search for transitions in the R–ρ relation. Here, we calculated the mean density from the measured M and R (see Sect. 2) and attempted to find a piece-wise linear function that describes the R–ρ relation. The data are shown in the top panel of Fig. 3. Qualitatively, three regimes can be identified. For the smallest planets, the density seems nearly independent of radius. Then, there is a breakpoint where the density decreases steeply with increasing radius. The giant planets (around 13 R_⊕) show a strong dispersion in density. This is similar to what we already observed in the M–R relation. For the giant planets, the bulk density can vary greatly due to their age, instellation flux, and composition. The intermediate and the giant planets start to overlap around 7 R_⊕, and the density appears almost uncorre-lated with the radius. Therefore, when fitting the R–ρ relation, we excluded planets larger than 7 R_⊕.

We found that the best model (with the lowest BIC) to describe the R–ρ relation uses one breakpoint. This is unlike the two breakpoints for the M–R and M–ρ relations. However, it is somewhat expected since we excluded the giant planets. The resulting R–ρ relation is $$\rho =\{\begin{array}{ll}(5.11\pm 0.19)R^{(0.73\pm 0.15)}\hfill & R<(1.64\pm 0.05)\hfill \\ (17.9\pm 1.5)R^{(-1.80\pm 0.17)}\hfill & R>(1.64\pm 0.05),\hfill \end{array}$$(7)

where ρ and R are in g cm⁻³ and R_⊕, respectively. The R–ρ best fit and the data are shown in the bottom panel of Fig. 3.

The values of the parameters for the piece-wise linear function with one breakpoint (see Eq. (1)) are listed in Table 3.

A notable result is the breakpoint at around 1.6 R_⊕. Its relative uncertainty of 3% is significantly lower than for the mass threshold between the small and intermediate planets derived from the M–R relation (14%) or the M–ρ relation (20%). This shows that it is beneficial to consider the radius when distinguishing between different planetary types (Rogers 2015; Lozovsky et al. 2018). Similar to the results for the M–ρ relation (see Sect. 3.2), the relative uncertainty of the R–ρ in the first segment is rather high, but it is consistent with a constant density approximation.

Lozovsky et al. (2018) found threshold radii above which a certain composition is unlikely. For purely rocky planets, they found a threshold radius of ${1.66}_{-0.08}^{+0.01}{R}_{\oplus }$. Larger planets must consist at least partly of lighter elements, such as H and He. This is consistent with our result of a breakpoint at 1.64 R_⊕. However, they only distinguished between super-Earths and mini-Neptunes at ≃3R_⊕, because planets with a larger radius have a substantial H-He atmosphere (at least 2% mass fraction). In contrast, based on our data no distinction can be made there.

Our result of the radius breakpoint at 1.64 R_⊕ also coincides with the position of the radius valley around 1.5–2 R_⊕. The radius valley is a bimodal feature in the occurrence rate of planets as a function of their radii, which manifests itself as a scarcity of planets with R ≃ 1.5–2R_⊕. It has been observed for planets with short periods (e.g., Fulton et al. 2017) and is often used for the distinction between super-Earths (below the valley) and mini-Neptunes (above the valley).

Several previous studies have shown how photoevaporation or core-powered mass loss can lead to the depletion of the gaseous envelopes of planets at such radii (Chen & Rogers 2016; Owen & Wu 2017; Venturini et al. 2020) and therefore explain the radius valley. In particular, Kubyshkina & Fossati (2022) suggested that the M–R relation of intermediate planets is shaped by their thermal evolution and hydrodynamic escape. Additionally, it has been suggested that the planets at the upper edge of the radius valley are helium-rich (Malsky et al. 2023). This suggests that due to evaporation of the gaseous envelope for masses with R ≃ 1.3 R_⊕ planets are naked rocky cores, while around 2.6 R_⊕ they sustain at least part of their H-He envelope. As an alternative, it has also been suggested that the bimodal radius distribution of planets smaller than about 4 R_⊕ is due to different compositions of rocky super-Earths and ice- or water-rich mini-Neptunes (Zeng et al. 2019; Venturini et al. 2020; Izidoro et al. 2021, 2022).

Fig. 3 Radius–density data from the PlanetS catalog and the derived relation. Top: radius-density distribution of all exoplanets from the PlanetS Catalog. Bottom: the exoplanets from the PlanetS Catalog smaller than 7 R⊕ are displayed together with the mass–density relation (solid black line) and the corresponding breakpoints (dashed line). The light-shaded areas are the 1σ-uncertainties.

3.4 Comparison with previous studies

A comparison of our results with previous studies is presented in Table 4. To facilitate the comparison, the relations and breakpoints were converted to Earth units (M_⊕ and R_⊕). For the M–R relation, we use the results from Sect. 3.1. Also, the mass-density relation from Hatzes & Rauer (2015) was converted to a M–R relation. From Edmondson et al. (2023) we chose the M–R relation for the giant planets instead of their mass-radius-temperature relation.

Overall, it is clear that there is a rather good agreement between the results of the various studies despite the use of different methods and data. The M–R relations from the different studies are shown in Fig. 4 together with the data from the PlanetS catalog. It can be seen that the relations from Weiss et al. (2013) and Bashi et al. (2017) underestimate the radii of the smallest planets and are not a good fit. This is because they only use one breakpoint in the M–R relation, corresponding to the transition to giant planets. The relation by Chen & Kipping (2017) does not fit the planets around 10 M_⊕ very well because the location of the transition from small to intermediate-mass planets is underestimated. Hatzes & Rauer (2015) did not fit planets below 95 M_⊕ at all. The relation by Otegi et al. (2020) remains a good fit for the dataset. The main difference compared to our results is the transition mass from small to intermediate planets. Otegi et al. (2020) defined the transition with the water-composition line, while we used a statistical approach to find this transition. The benefit of our approach is that it does not rely on a priori assumptions or theoretical models to determine the transition (and the associated uncertainties in, e.g., the EOS). Perhaps surprisingly, despite the spread of the planets in the M–R diagram, our dataset yielded a rather small uncertainty in the transition masses. Similarly, Edmondson et al. (2023) used a pure-ice EOS to mark the transition between small and intermediate planets, which leads to a good description of the smallest planets and intermediate planets. However, from ~30– 90 M_⊕ their relation for icy planets significantly under-predicts the radii of the planets in the PlanetS catalog, leading to a poor fit. Compared to all the previous small-planet M–R fits listed, our uncertainty on the power-law index is slightly higher. This is likely because the small planets in our updated data have diverse radii. For the giant exoplanets, all the relations can qualitatively describe the M–R relation, although they are quite different and can even have a different sign of the gradient. The data in this regime show a strong dispersion, which leads to rather large uncertainties in the fit relation. Interestingly, using the M–ρ relation to determine the transition to the giant planets yielded a significantly larger mass (183 M_⊕) than both the M–R fit (127 M_⊕) and results from previous studies (95 to 150 M_⊕).

Fig. 4 Comparison of M–R relations from different studies with exo-planetary data from the PlanetS catalog.

4 Discussion and conclusions

In this work, we used the updated PlanetS catalog to infer M–R, M–ρ, and R–ρ relations and determined the transitions between different planetary types. While the presented analysis provides insight into the different planet regimes, it was simplified and did not consider all the subtleties related to exoplanetary data. First, we treated the data as one unit, although it is clear that the dataset is inhomogeneous and combines different observational methods with different biases. The effects of observational bias, for the most part, have not been considered. Other parameters affect the M–R relation that was not investigated in this work. For example, for giant planets, stellar age and irradiation are important. Giant planets are massive enough that their self-gravity causes them to contract over long timescales (~1 Gyr; e.g., Hubbard 1977; Burrows et al. 2001), and therefore their radius is expected to be correlated with their age. Additionally, high instellation fluxes inflate the radii of warm giant planets (e.g., Guillot et al. 2006; Fortney et al. 2007; Fortney & Nettelmann 2010; Thorngren et al. 2016; Müller & Helled 2023). This effect was included in Weiss et al. (2013) and Edmondson et al. (2023), where a third parameter (instellation flux or equilibrium temperature) was added to better fit the M–R of the giant planets. Recently, there have also been attempts to move beyond the two-parametric M–R relationship. For example, Kanodia et al. (2023) presented a framework to characterize exoplanets using up to four simultaneous parameters. In the future, such approaches may better constrain the transition from small to intermediate planets from observational data by considering additional parameters.

Out of the over 5000 detected exoplanets, only 688 have robust enough mass and radius measurements to be included in the PlanetS catalog. While this means that only a fraction of the currently detected exoplanets were used in this work, the results are also more robust, since only planets with low mass and radius uncertainties are included. More accurate data are needed to analyze the whole parameter space occupied by exoplanets.

The key results from our study can be summarized as follows:

1. Our analysis yielded a small-to-intermediate transition mass of (4.37 ± 0.72) M_⊕. Small planets below the transition mass follow R ∝ M^0.27. These are “rocky worlds” with different bulk compositions. The transition to the intermediate mass could imply a maximal mass of ~4.4 M_⊕ of “rocky” exoplanets and naked planetary cores;

2. The transition from rocky to volatile-rich planets can also be defined in terms of the radius. By fitting the radius-density relation, we found that the transition occurs around 1.64 ± 0.05 R_⊕. The transition in radius has a significantly lower relative uncertainty than the one in mass. Furthermore, the transition radius is consistent with the radius valley around 1.5 to 2 M_⊕;

3. Intermediate-mass planets ranging from about 4.4 M_⊕ to 127 M_⊕ behave as R ∝ M^0.67. They correspond to planets with H-He envelopes. The transition to giant planets occurs at (127 ± 7) M_⊕ and corresponds to planets that are H-He-rich;

4. Using the M–ρ relation to find the transition to the giant planets yielded a significantly higher transition mass of (183 ± 7) M_⊕;

5. The radii of giant planets are nearly independent of their masses, and the mass-radius relation in this regime follows R ∝ M^–0.06;

6. Overall, planets of different compositions and structures can have the same mass and radius. This leads to an intrinsic degeneracy of the mass-radius distribution of exoplanets.

Ongoing and future observations on the ground and in space will improve our understanding of exoplanets. The James Webb Space Telescope (Gardner et al. 2006) and the Ariel mission (Tinetti et al. 2018) will enable us to characterize the atmospheres of transiting planets, providing information about their chemical compositions. High-resolution spectroscopy from current (e.g., SPIROU; Artigau et al. 2014, CARMENES; Quirrenbach et al. 2016) and future (e.g., NIRPS; Bouchy et al. 2017; Wildi et al. 2017, CRIRES+; Kaeufl et al. 2004; Dorn et al. 2014, 2023) ground-based telescopes will provide further improvements with accurate radial-velocity measurements and atmospheric characterizations.

More exoplanets detected via direct imaging, for example by SPHERE at the Very Large Telescope (Beuzit et al. 2019), will facilitate studies of the properties of planets on wide orbits. Also, the upcoming PLATO mission (Rauer et al. 2014) will detect and characterize small terrestrial planets as well as intermediate-mass and giant planets. Theoretical studies to understand the key physical processes that shape the exoplanetary populations are also being developed, and we hope to be able to connect the properties of exoplanets with their origin and evolution. These ongoing and upcoming efforts are expected to reveal new insights into the population of planets beyond the Solar System.

Acknowledgements

We thank the anonymous reviewer for useful feedback. We acknowledge support from SNSF grant 200020_215634 and the National Centre for Competence in Research ‘PlanetS’ supported by SNSF. This research used data from 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. Extensive use was also made of the Python packages NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Jupyter (Kluyver et al. 2016), Matplotlib (Hunter 2007), and pandas (The pandas development team 2020).

References

All Tables

All Figures

	Fig. 1 Mass-radius data from the PlanetS catalog and the derived relation. Data from PlanetS catalog is displayed in purple. The best-fit mass–radius relation (Eq. (5)) is represented by the solid black line and the shaded region shows the 1σ-uncertainty. The breakpoints are shown as the dashed lines.
In the text

Fig. 2 Mass–density and mass–radius data from the PlanetS catalog and the derived relations. Top: mass and density data of exoplanets from PlanetS catalog. The mass–density relation is represented by a solid black line. The dashed lines show the position of the breakpoints. The light shaded areas are the corresponding 1σ uncertainties. Bottom: mass versus radius distribution with best fit for relation and breakpoints in black, as in Fig. 1. The orange solid and dashed lines show the mass-radius relation and the corresponding breakpoints derived from the fit to the mass-density distribution.

In the text

	Fig. 3 Radius–density data from the PlanetS catalog and the derived relation. Top: radius-density distribution of all exoplanets from the PlanetS Catalog. Bottom: the exoplanets from the PlanetS Catalog smaller than 7 R⊕ are displayed together with the mass–density relation (solid black line) and the corresponding breakpoints (dashed line). The light-shaded areas are the 1σ-uncertainties.
In the text

	Fig. 4 Comparison of M–R relations from different studies with exo-planetary data from the PlanetS catalog.
In the text