© The Authors 2024

1. Introduction

Exoplanets with sizes between Earth and Neptune are the most abundant type known until now (e.g. Petigura et al. 2022). The size distribution of these objects is bimodal, with a “radius valley” present at ∼1.5 − 2 R_⊕ separating the smaller super-Earths from the larger mini- (or sub-)Neptunes (Fulton et al. 2017; Fulton & Petigura 2018; Van Eylen et al. 2018; Martinez et al. 2019). The exact location of the radius valley depends on the stellar mass and planets’ orbital period, with different data analysis yielding different slopes in the log(R_P)−log(M_⋆) and log(R_P)−log(P_orb) planes (with R_P, P_orb, and M_⋆ as the planet radius, orbital period, and stellar mass, respectively; Fulton et al. 2017; Fulton & Petigura 2018; Berger et al. 2020; Petigura et al. 2022; Luque & Pallé 2022; Ho & Van Eylen 2023; Bonfanti et al. 2024).

Understanding the origin of the radius valley and its dependence on orbital and physical parameters has become a crucial endeavour in modern exoplanetology (e.g. Owen & Wu 2017; Ginzburg et al. 2018; Gupta & Schlichting 2019; Venturini et al. 2020a; Lee & Connors 2021). The first mechanisms to account for the existence of the radius valley were the pure evolution models known as photoevaporation (e.g Lopez & Fortney 2013; Chen & Rogers 2016; Owen & Wu 2017; Jin & Mordasini 2018) and core-powered mass-loss (Ginzburg et al. 2018; Gupta & Schlichting 2019). Based on these models, super-Earths and mini-Neptunes would originate from the same single-composition population of rocky cores with H/He atmospheres. The distinction between super-Earths and mini-Neptunes arises because some planets lose their atmospheres completely during the evolution, while others retain ∼1% by mass in terms of H/He. Both models find slopes in the radius-period and radius-stellar mass plane in agreement with observations (Gupta & Schlichting 2019; Rogers et al. 2021). However, these pure evolution models are in contradiction with planet formation theory, which predicts most mini-Neptunes to have formed beyond the water ic line and to therefore be water-rich (Alibert et al. 2013; Raymond et al. 2018; Bitsch et al. 2018; Izidoro et al. 2021; Brügger et al. 2020; Venturini et al. 2020a). In Venturini et al. (2020a, hereafter V20), we performed pebble-based planet formation simulations for Sun-like stars, including a self-consistent treatment of pebble growth. We found that because of the different physical properties of icy versus rocky pebbles, icy cores are typically born bigger than rocky ones. In terms of total planet radius, this bi-modality is hidden at the time of disc dissipation because most planets accrete gas that puffs the planetary radii up. The separation between rocky and water-worlds becomes clear once atmospheric mass-loss sets in, stripping the atmospheres of the small rocky and icy planets and unveiling the primordial radius valley separating the two planetary types, at ∼1.5–2 R_⊕ (Venturini et al. 2020a). In this work, we adapt and apply our planet formation and evolution model to a wide range of stellar masses, from 0.1 to 1.5 M_⊙ (i.e. stellar spectral types spanning M dwarfs to A types). The upper limit on stellar mass is motivated by the earliest-stellar-type exoplanets discovered by Kepler (e.g. Berger et al. 2020). On the theoretical side, we note that no population study has addressed the origin of planets orbiting super-solar stellar masses until now. This is important to address given the new TESS data on exoplanets around F and A stars (Psaridi et al. 2022).

Regarding M planets, the motivation to study them is based on the increasing number of facilities targeting such objects, such as TESS (Ricker et al. 2015), CHEOPS (Benz et al. 2021), and NIRPS (Bouchy et al. 2022). In particular, regarding the radius valley, several studies are reporting masses and radii of super-Earths and mini-Neptunes (e.g. Demory et al. 2020; Van Eylen et al. 2021; Luque & Pallé 2022; Bonfanti et al. 2024), and the mere existence of the radius valley for M dwarfs remains controversial (Cloutier & Menou 2020; Luque & Pallé 2022; Bonfanti et al. 2024). From a theoretical perspective, Alibert (2017), Miguel et al. (2020), and Burn et al. (2021) showed that low-mass water-rich planets are expected to be common at close-in orbits around M dwarfs.

In this work, we conduct a large planet formation and evolution parameter study, focusing on the population of low- and intermediate-mass planets. The aim is to devise a theoretical description of the behaviour of the radius valley for different stellar types. Our main results are highlighted in the brief main text, supplemented by several in-depth analysis in the Appendices.

2. Methods

We computed the formation and evolution of planets around stars of mass 0.1, 0.4, 0.7, 1.0, 1.3, and 1.5 M_⊙. In our model (Venturini et al. 2020a,b), planets grow from a lunar-mass embryo by pebble and gas accretion, while embedded in a gaseous disc that evolves by viscous accretion and X-ray photoevaporation from the central star. The discs are adapted to different stellar masses following Burn et al. (2021) (details in Appendix A.1). The evolution of the pebbles was calculated considering pebble coagulation, fragmentation, drift, diffusion, and ice sublimation at the water ice line (model of Birnstiel et al. 2011; Drążkowska et al. 2016). The planets can accrete either rocky or icy pebbles, depending on its position with respect to the water ice line (which evolves in time, moving inwards as the disc cools down). Core growth halts when pebble isolation mass is reached (Lambrechts et al. 2014). As in V20, we adopted a fragmentation threshold velocity of icy and silicate pebbles as v_th = 10 m s⁻¹ and v_th = 1 m s⁻¹, respectively (e.g. Gundlach & Blum 2015). We note, based on experimental works, that these values are still quite uncertain (see detailed discussion in Appendix C.1.2). Our results hold as long as the fragmentation threshold velocity of icy pebbles is in the range of 5–10 m s⁻¹. More experimental work and sensitivity studies are required to better constrain these fragmentation velocities. For the disc viscosity, we consider α to either be 10⁻³ or 10⁻⁴. The details about the choice of the initial conditions are given in Appendix A.1. The planets also migrate along the disc, either via the type I-or type-II regime. The former includes the Lindblad, corotation, and thermal torques (Guilera et al. 2019).

As explained in V20, since our code can handle at the moment the formation of only one embryo per disc, and given the fast timescales of pebble accretion and type-I migration, most of our formed planets get stranded near the disc inner edge, which is set at an orbital period of 11.55 days for all the considered discs (stemming from the choice of disc inner edge at 0.1 au for solar-type stars, see Appendix A.1). This means that most of the results reported throughout this work refer to planets with orbital period of ∼11–18 days.

Because of the limitation of considering only one embryo per disc, the possibility of atmospheric mass loss due to giant impacts is in principle hindered in our model. Nevertheless, as in V20, we estimated this rate in a simplified way, following Inamdar & Schlichting (2015) and Ronco et al. (2017), by considering a potential impact between the simulated planet and another, less massive one formed in a different simulation, but under the same disc initial conditions (see Appendix A.1). We refer to this batch of simulations as the “collision case” (while the standard simulations without the estimation of atmospheric mass-loss by collisions is referred as “nominal case”). The details of this implementation and its limitations can be found in Appendix A.3.

Once the disc dissipates, the cooling of the planets (including the cooling of the core) is calculated during 2 Gyr with photoevaporation following Mordasini (2020). We refer to this phase as “evolution”, while “formation” refers to the disc-embedded phase when planetary accretion takes place. During evolution, we consider the planets to consist of a rocky core of an Earth-like composition, surrounded by an envelope of volatiles which can contain mixtures of H/He with water (with all the components and amount stemming from the formation process). The atmospheric mass-loss occurs for all the components of the volatile layer. Further details of the evolution model are described in Appendix A.4. Finally, we computed the planetary radii after 2 Gyr of evolution with state-of-the-art interior structure models (details in Appendix A.5). The evolution and structure calculations were performed for all the planets that ended with orbital periods below 100 days after formation and these are the cases analysed and presented throughout this study.

3. Results

3.1. Mass-radius diagrams: General trends

We analysed the output of our formation-evolution simulations at time 2 Gyr. Figure 1 shows the mass-radius (MR) diagram at that time, for the different stellar masses, colour-coded with the planet’s water mass fraction (total water mass divided by the total planet mass). The dotted lines represent the MR relations for an Earth-like composition (brown), and for planets composed of 50% water and 50% rocks by mass, with the water either in condensed form (blue) or in the form of steam (light-blue). Hereafter, we refer to these compositional curves as the “condensed-water line” and the “steam-line”, respectively. The derivation of these M–R relations is detailed in Appendix B.5.

Fig. 1. MR diagram for different stellar masses. The colour bar represents the total water mass fraction in the planets. Dots correspond to the nominal output, diamonds to collisions. The brown- and blue-dashed lines correspond to a theoretical composition of Earth-like and 50% water + 50% Earth-like, as defined in Appendix B.5, with dark-blue indicating condensed-water worlds and light-blue indicating steam worlds. The green-dotted curves correspond to 50% steam + 50% Earth-like from Aguichine et al. (2021) (Teq = 400 K for dark-green and Teq = 600 K for light-green). Grey dots are real exoplanets with mass and radius measurements better than 75 and 20%, respectively (data taken from the NASA Exoplanet Archive on 08.09.23). The real exoplanets in each panel correspond to the stellar masses of 0–0.25 M⊙ for M⋆ = 0.1 M⊙, 0.25–0.55 M⊙ for M⋆ = 0.4 M⊙, 0.55–0.85 M⊙ for M⋆ = 0.7 M⊙, 0.85–1.15 M⊙ for M⋆ = 1.0 M⊙, 1.15–1.4 M⊙ for M⋆ = 1.3 M⊙, and 1.4–1.7 M⊙ for M⋆ = 1.5 M⊙.

A general trend that we find for all panels, is a division given by the rocky and condensed-water lines. Bare rocky cores cluster around the Earth-like composition (as they should based on our assumptions) and wet planets tend to be either on the condensed-water line or above. This diagonal band between the rocky and condensed-water lines remains planet-depleted for all stellar masses (although some planets with intermediate compositions exist in between those lines, particularly for M dwarfs; see Fig. 1). Planets with radii larger than the one given by the condensed-water line tend to have substantial H/He, particularly when moving towards larger planet mass. This can be appreciated as a fading of the blue colour towards large planetary radii (the mass fraction of water decreases as the H/He mass fraction increases). The H/He mass fraction of each synthetic planet is also depicted in Fig. B.4. Another clear trend of Fig. 1 is that for M_⋆ ≥ 0.7 M_⊙, a large number of water-worlds cluster around the steam line (particularly for M_⋆ = 0.7 M_⊙). Indeed, for those cases the temperatures are high enough for water to be present in the form of steam throughout the atmosphere. This is demonstrated in Appendix B.4. Other aspect analysed in detail in Appendix B.7 is the presence of “red” outliers among “blue” planets. These are planets that actually formed beyond the ice line but lost most of their water during evolution.

We also note from Fig. 1 that the lower the stellar mass, the larger the overlap in mass and radius of dry and wet planets. As the stellar mass increases, the minimum mass of wet planets increases as well, taking values of M_min≃ 0.6, 1.4, 2.8, 4.5, 5.4 and 7.3 M_⊕ for M_⋆ = 0.1, 0.4, 0.7, 1.0, 1.3, 1.5 M_⊙, respectively¹. This is the imprint of orbital migration. Indeed, the lower the stellar mass, the lower the threshold mass to undergo type-I migration when embedded in a gaseous disc (Paardekooper et al. 2011). The same effect was reported in Burn et al. (2021). This produces a large overlap in mass and radius between dry and wet planets for M_⋆ = 0.1 M_⊙, filling the radius valley, as noted in the histograms of Fig. B.2 for M_⋆ ≲ 0.4 M_⊙. The radius valley starts to clearly emerge for M_⋆ ≳ 0.7 M_⊙, where water worlds have typical masses above 2–3 M_⊕.

3.2. Valley in the radius-stellar mass plane

Different analysis of exoplanet data have shown that the radius valley follows a power-law dependence with stellar mass, namely,

$$\begin{array}{c}\hfill R_{\mathrm{valley}}=R_{\mathrm{v},0}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^m.\end{array}$$(1)

For KGF stars from the Kepler data, Berger et al. (2020) found $m=0.{26}_{-0.16}^{+0.21}$, and Petigura et al. (2022), $R_{\mathrm{v},0}=1.{86}_{-0.03}^{+0.03}{R}_{\oplus }$, $m=0.{18}_{-0.07}^{+0.08}$. For M dwarfs, Luque & Pallé (2022) found m = 0.08 ± 0.12.

In Fig. 2, we plot the planet radius as a function of the stellar mass for our synthetic planets, colour-coded by f_H2O, as in Fig. 1. The blue-dotted line represents our best fit of the valley adopting the power-law dependence of Eq. (1). The location of the valley was found using the open tool https://github.com/parkus/gapfit (Loyd et al. 2020).

Fig. 2. Radius valley fit as a function of the stellar mass (blue line with errors in light-blue) for all our synthetic planets (colour-coded as a function of the water mass fraction). Dots correspond to the nominal output, diamonds to collisions. The green dashed line with the light-green shaded area show the best fit found by Petigura et al. (2022), the violet dashed line and the lilac shaded area represent the results found by Berger et al. (2020), and the mustard dashed line and the light-mustard shaded area represents the best fit by Luque & Pallé (2022). Our best fit was found using https://github.com/parkus/gapfit (Loyd et al. 2020).

By adopting R_V0, guess = 1.82 and m_guess = 0.09 in gapfit, our best fit yields R_v, 0 = 1.8197 and $m=0.{14}_{-0.01}^{+0.02}$. We note that we are not plotting the 1σ error computed with gapfit, but the best fit with the corresponding errors of the constants R_v, 0 and m. We also over-plot the fits from previous studies with their corresponding errors in Fig. 2 (see figure caption). Our theoretical results are in excellent agreement with the observational results, for all the stellar masses ranging from 0.1 to 1.5 M_⊙. Figure 2 also shows that for M_⋆ ≲ 0.4 M_⊙, the radius valley fades (or, in other words, it gets filled), because of the increasing overlap between rocky and water worlds (see Fig. 1). A vanishing radius valley for M dwarfs is observationally supported by the findings of Luque & Pallé (2022) and Ho et al. (2024). Another remarkable aspect of Fig. 2 is that for M_⋆ ≳ 0.7 M_⊙, the radius valley is actually empty, particularly for M_⋆ ≥ 1.3 M_⊙. This is in agreement with the findings of Van Eylen et al. (2021) and Ho & Van Eylen (2023).

3.3. Valley in the radius-orbital period plane

We additionally analysed the location of the valley as a function of orbital period using gapfit for KGF stars and M dwarfs. The result of the best fit for KGF-stars is $R_{\mathrm{valley}}(R_{\oplus })=2.5704P_{\mathrm{orb}}^{-0.16\pm 0.01}$ and is shown in Fig. 3, together with the fit derived by Petigura et al. (2022), $R_{\mathrm{valley}}(R_{\oplus })=2.{37}^{\pm 0.04}{P}_{\mathrm{orb}}^{-0.11\pm 0.02}$. We note that including A type stars (M_⋆ = 1.5 M_⊙) does not alter the fit. Despite the limitations of our model (too many planets concentrated at P_orb ≈ 11 − 15 days, see Sect. 2 and Appendix A.6), we find a remarkable good agreement with observations, with a negative slope carved by photoevaporation. For M dwarfs (M_⋆ = 0.1 and 0.4 M_⊙), we find a slope of $-0.{04}_{-0.03}^{+0.08}$, i.e. a flat slope, in line with Luque & Pallé (2022) – who found a slope of −0.02 ± 0.05 –, and also in agreement (within errors) with Bonfanti et al. (2024), who derived a slope of $-0.{065}_{-0.013}^{+0.024}$. However, we caution that in this case the package gapfit cannot easily find a gap, with different initial guesses yielding different slopes. This is not a problem related to a lack of points (with other pairs of stellar bins this problem does not arise), which reinforces the conclusion that the valley gets indeed filled when moving towards M dwarfs.

Fig. 3. Radius valley fit as a function of orbital period (blue line with errors in light-blue). The colour bar indicates total water mass fraction. Top panel: Planets around KGF-type stars. The green line and shade show the best fit found by Petigura et al. (2022) for these stellar types. Bottom panel: Planets around M dwarfs. The green line and shade show the best fit found by Luque & Pallé (2022). The particular clustering in the orbital period in both panels is due to the model limitations (Appendix A.6).

4. Discussion: Compositional valley

We computed the growth and post-formation evolution of planets orbiting starts from 0.1 to 1.5 M_⊙. Several model limitations exist, which we discuss in detail in Appendix A.6. Overall, our findings reinforce our conclusion of V20 that the radius valley is a separator in planetary composition between bare rocky cores whose atmospheres were stripped off by photoevaporation and water worlds that formed beyond the ice line and retained some or none of their primordial H/He. However, the radius valley fades with decreasing stellar mass (Figs.1–3). This is a consequence of orbital migration, which happens for less massive planets at lower stellar masses (Paardekooper et al. 2011; Burn et al. 2021). Hence, for low mass-stars, originally small icy planets that formed beyond the ice line can reach the inner regions of the system, rendering a larger overlap in mass and radius between planets formed inside and outside the ice line for M dwarfs.

Nevertheless, the paucity of planets with intermediate values of water mass fraction (0 < f_H2O ≲ 0.2) persists for all stellar masses, yielding a depleted (M_⋆ ≤ 0.4 M_⊙) or empty (M_⋆ ≥ 0.7 M_⊙) diagonal band in the MR diagram of Fig. 1 between the rocky and condensed-water lines. This is due to pebble accretion, which proceeds quickly beyond the ice line, typically allowing the planets to reach the pebble isolation mass before crossing it. This means that solid accretion tends to stops before the icy planets have the possibility to accrete dry pebbles and mix their bulk composition (particularly true for M_⋆ ≳ 1.0 M_⊙). If this “diagonal valley” were to not be confirmed by future observations, then planetesimal-based (e.g. Brügger et al. 2020) or hybrid pebble-planetesimal models (Alibert et al. 2018; Guilera et al. 2020) would be favoured. Alternatively, it could indicate a slower pebble accretion resulting from a lower fragmentation velocities of icy pebbles. This “diagonal valley” in the M-R diagram appears more clearly when plotting, instead of radius, mean density (normalised by the density of an Earth-like composition, as in Luque & Pallé 2022, also see Appendix B.6). We show this in Fig. 4, including the compositional lines of Earth-like, pure-silicate composition, condensed-water, and steam worlds. The colour bars indicate the planetary equilibrium temperature and the different symbols give the different stellar masses (0.1–0.7 M_⊙ on the left panel, 1.0–1.5 M_⊙ on the right). A clear pattern emerging from Fig. 4 is a depletion (void) of planets with ρ ≈ 0.6 − 0.9 ρ_⊕, S around stars of 0.1–0.7 M_⊙ (1.0–1.5 M_⊙). Our model assumes an Earth-like composition for all the rocky components. In reality, different stellar abundances would yield rocky planets with different iron fractions, implying that rocky planets could be expected down to the pure silicate line (in orange) of Fig. 4. This means that, in reality, the density valley would be more blurred, but should still exist at ρ ≈ 0.6 − 0.8 ρ_⊕, S as a consequence of pebble accretion. This is overall in line with the conclusions of Luque & Pallé (2022), who found a valley at ρ ≈ 0.65 ρ_⊕, S for M-planets.

Fig. 4. Density-mass diagram (density normalised by an Earth-like composition) as a function of equilibrium temperature (colour bar). Left panel: Planets orbiting stars of 0.1 M⊙ (circles), 0.4 M⊙ (squares), and 0.7 M⊙ (diamonds). Right panel: Planets orbiting stars of 1.0 M⊙ (circles), 1.3 M⊙ (squares), and 1.5 M⊙ (diamonds). Note: the colour bars span different ranges in the two panels. For both panels, black-bordered symbols indicate fHHe < 1%, dark-grey 1 ≤ fHHe ≤ 10%, and light-grey fHHe > 10%. The dashed-lines show the compositional curves of Earth-like (brown), pure-silicate (orange), condensed water (blue), and steam worlds (light-blue). The grey-dotted lines correspond to radii of 1 and 4 R⊕.

We note as well that water worlds with “condensed” water exist mainly for M_⋆ = 0.1 and 0.4 M_⊙, where the equilibrium temperatures are below ∼400 K (Fig. 4, bottom left, and Fig. C.2). For 0.7 M_⊙, water-worlds are very common, but water is present as steam along all the planetary envelope. Thus, for planets with T_eq ≳ 400 K, we predict that water worlds should cluster around the line of 50% rock and 50% steam by mass (at approx. 0.2–0.4 ρ_⊕, S), as long as planetary cores are built predominately from pebbles. This steam line has planetary radii larger than the condensed-water line by 30–50% for planetary masses of 1 − 5 M_⊕, following the tendency reported by Turbet et al. (2019) and Aguichine et al. (2021). This means that the existence of steam-worlds helps to carve a deeper radius valley for T_eq ≳ 400 K (also reported by Burn et al. 2024, for planets around Sun-like stars, with T_eq ≳ 600 K). In this regard, our results do not entirely match those of Luque & Pallé (2022), who report water worlds of T_eq > 400 K falling on the condensed-water line. A possible explanation might be steam atmospheres with a reduced water mass fraction (see discussion in Appendix C.3). Increasing the sample of M planets will be crucial for testing our model predictions.

Finally, Figs. 1 and 4 also illustrate how mini-Neptunes range from being predominantly water-condensed worlds for M_⋆ = 0.1 and 0.4 M_⊙, to steam worlds for M_⋆ = 0.7 M_⊙, to steam-H/He planets for M_⋆ ≳ 1.0 M_⊙, with a H/He mass fraction below 10%. This prediction holds within the ranges in equilibrium temperature considered in this work (see Fig. 4). The case of M_⋆ = 1.5 M_⊙ appears somehow to be the exception, where several mini-Neptunes are dry due to the high loss of water during evolution (see Fig. 1 and Appendix B.7).

5. Conclusions

We studied the combined formation and evolution of planets around single stars in the mass range 0.1 ≤ M_⋆ ≤ 1.5 M_⊙, with the goal of characterising the radius valley for different stellar types. Overall, we find that the tendency (found in V20) of having larger water worlds over smaller rocky planets persists for all stellar ranges, but the radius valley separating them fades towards M dwarfs (in agreement with Luque & Pallé 2022; Bonfanti et al. 2024). This is a consequence of orbital migration, which allows small icy planets to reach the disc inner regions for low-mass stars, producing a larger overlap in mass and radius between the rocky and icy populations.

Despite of the “filling” of the radius valley towards low stellar masses, we find that when considering the full range of stellar masses, its dependence on stellar mass is in excellent agreement with observations, following $\mathrm{d}\log R_{\mathrm{valley}}/\mathrm{d}{M}_{\star }=0.{14}_{-0.01}^{+0.02}$. We also find a good agreement with observations regarding the location of the radius valley with orbital period (see Sect. 3.3), with a negative slope carved by photoevaporation for KGF-stars and a flat slope for M dwarfs.

Our end-to-end simulations show that super-Earths are bare rocky cores that emerge as evaporated worlds around all stars. On the other hand, mini-Neptunes, are, in their vast majority, water worlds that lost all or part of their primordial H/He due to photoevaporation. In addition, we confirm that the particular phase of water is key in shaping the radii of the exoplanet population. If the outer atmospheric temperatures are low enough for water to condense (T_eq ≲ 400 K), then a water-world will transition from a steam to an icy or liquid world, reducing its size by ∼15% at a planetary mass of 5 M_⊕.

Finally, within our pebble-based formation and evolution model, the radius valley emerges as a divider in planetary composition. Thus, we note that it is more visible in terms of mean density, where a clear valley occurs at normalised mean density of ρ/ρ_⊕, S ∼ 0.6 − 0.8 across all stellar types.

Acknowledgments

We thank the anonymous referee for insightful criticism. J.V. acknowledges support from the Swiss National Science Foundation (SNSF) under grant PZ00P2_208945. J.H. and C.M. acknowledge support from the Swiss National Science Foundation (SNSF) under grant 200021_204847 “PlanetsInTime”. M.P.R. is partially supported by PICT-2021-I-INVI-00161 from ANPCyT, Argentina. M.P.R., O.M.G. and M.M.B. are partially supported by PIP-2971 from CONICET (Argentina) and by PICT 2020-03316 from Agencia I+D+i (Argentina). M.P.R. and O.M.G. acknowledge support by ANID, – Millennium Science Initiative Program – NCN19_171. O.M.G. gratefully acknowledges the invitation and financial support from ISSI Bern in early 2024 through the Visiting Scientist Program. Software: For this publication the following software packages have been used: Python-matplotlib (https://matplotlib.org/) by Hunter (2007), Python-seaborn (https://seaborn.pydata.org/) by Waskom et al. (2020), Python-numpy (https://numpy.org/), Python-pandas (https://pandas.pydata.org/).

References

Appendix A: Methodology

A.1. Formation model: Initial conditions

We applied the global planet formation model PLANETALP (Ronco et al. 2017; Guilera & Sándor 2017; Guilera et al. 2019) with stellar masses of 0.1, 0.4, 0.7, 1.0, 1.3, and 1.5 M_⊙. As in V20, we used the 12 protoplanetary disc profiles from Andrews et al. (2010), with two uniform viscosity parameters (α = 10⁻⁴, 10⁻³), where the initial gas surface density is given by:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{g}}={\mathrm{\Sigma }}_{\mathrm{g},0}{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{-\gamma }{e}^{-{(r/r_{\mathrm{c}})}^{2-\gamma }},\end{array}$$(A.1)

where Σ_g, 0 is a normalisation parameter determined by the disc initial mass (M_d, 0), γ is the exponent that represents the surface density gradient, and r_c is the characteristic or cut-off radius of the disc. Following Burn et al. (2021), we scale the initial disc mass and characteristic radii as:

$$\begin{array}{cc}& {\mathrm{M}}_{\mathrm{d},0}=\left(\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_{\odot }}\right){\mathrm{M}}_{\mathrm{d},0}^{\mathrm{sun}},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& r_{\mathrm{c}}={\left(\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{1/1.6}{r}_{\mathrm{c}}^{\mathrm{sun}},\hfill \end{array}$$(A.3)

where ${\mathrm{M}}_{\mathrm{d},0}^{\mathrm{sun}}$ and $r_{\mathrm{c}}^{\mathrm{sun}}$ stem from Andrews et al. (2010) (displayed in Tab. B.1 of Venturini et al. 2020a). In Venturini et al. (2020a), we considered only solar-mass stars and we set the inner edge of the disc at 0.1 au, based on hydrodinamical simulations (Flock et al. 2019) and on the mean orbital period of the innermost planet of planetary systems (Mulders et al. 2018). Also, as in Burn et al. (2021), we scale the disc inner edge as:

$$\begin{array}{c}\hfill {\mathrm{r}}_{\mathrm{int}}={\left(\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{1/3}0.1\mathrm{au}.\end{array}$$(A.4)

This scaling is based on identifying the location of r_int as the point where the Keplerian orbital period matches the stellar rotation period (assumed to be the same for all stars as in Burn et al. 2021). However, we note that the location of the disc inner edge and its scaling with stellar mass is very uncertain.

Given the initial gas surface density by Eq. (A.1), the initial dust surface density is defined as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{d}}={\eta }_{\mathrm{ice}}{Z}_0{\mathrm{\Sigma }}_{\mathrm{g}},\end{array}$$(A.5)

where η_ice represents the sublimation/condensation of water-ice and adopts values of η_ice = 1/2 inside the ice line and η_ice = 1 outside of it (Lodders et al. 2009). Z₀ is the initial dust-to-gas ratio (adopting values of 0.0068, 0.0099, 0.0144, 0.0210, and 0.0305) to span the known ranges of stellar metallicities (as assumed in V20). Regarding the thermodynamical quantities of the discs around stars with different masses, we used the corresponding effective temperatures, T_eff, and stellar radii, R_⋆, from Baraffe et al. (2015) for the different protostars at 0.5 My to compute the disc vertical structure (as in Guilera et al. 2017, 2019).

As in Venturini et al. (2020a), we launch seven embryos per disc (one at a time), three of them are inside and four of them beyond the ice line. For the embryos that start within the ice line in the solar-case, the initial semi-major axis is defined with uniform log-spacing between r_int and r_ice − 0.1 au, and between r_ice + 0.1 au and 16 au for the icy embryos. For the non-solar cases, the corresponding initial semimajor axis scales as a_ini(M_⋆/M_⊙)^1/3 as the disc’s inner edge.

A.2. Gas accretion: Calculation and assumptions

Gas accretion is computed both in the attached and detached phases as in Venturini et al. (2020b). During the attached phase (where the planet’s envelope connects smoothly with the gaseous disc), gas accretion is calculated by solving the 1D spherically symmetric structure equations for the envelope before the planet reaches the pebble isolation mass (following the method of Alibert & Venturini (2019), which uses deep neural networks trained on pre-computed structure models for sub-critical core masses). When solving the structure equations, the envelope opacity is taken from Freedman et al. (2014) for the gas, and from Bell & Lin (1994) for the dust (see the discussion on the choice of opacities in Appendix C.1.3.) Once the pebble isolation mass is attained, solid accretion stalls and the core becomes critical. At this stage, we adopt the prescription of Ikoma et al. (2000), which is valid precisely in this regime when solid accretion is halted. After reaching the pebble isolation mass, at a certain point, the planet would accrete more gas than what the disc can provide due to its viscosity. Here, the planet detaches from the disc and accretes at a rate given by the disc viscous accretion (3πνΣ). In the detached phase, the gas accretion can be damped even further when the planet opens a gap. This is taken into account following Eqs. (36)-(39) of Tanigawa & Ikoma (2007).

A.3. Gas removal by giant impacts

Once the gas of the disc has completely dispersed, during the first million years of evolution, giant collisions between the remaining planets and embryos may become an important mechanism for efficiently removing the planet’s atmosphere (Ogihara & Hori 2020; Ogihara et al. 2020) before atmospheric mass loss due to photoevaporation takes place (Izidoro et al. 2017). Here, as in Venturini et al. (2020a), we formed seven embryos per disc, but only one at a time, which implies that N-body interactions between protoplanets embedded in the gaseous disc are not modelled. However, we can simply estimate the atmospheric mass-loss due to potential collisions between the formed planet and another one less massive with final period < 100 days, also formed in the same disc, but in a different simulation. For the sake of simplicity and in alignment with the findings of Ogihara & Hori (2020), who observed one or two giant impacts in their N-body analysis, we restricted each planet to a single collision. We follow the procedure developed in Ronco et al. (2017), where the core mass of the planet after the collision is the sum of the target’s and impactor’s cores, and where the final envelope mass (M_E) is estimated with the global atmospheric mass-loss fractions (X_loss) computed in Inamdar & Schlichting (2015) (their Fig.5) between super-Earths and mini-Neptunes. Thus, if $M_{\text{E}}^i$ and $M_{\text{E}}^j$ are the envelope masses of the target, i, and the impactor, j, the envelope mass of the target after the collision is given by $M_{\text{E}}^i=M_{\text{E}}^i(1-X_{\text{loss}}^i)+M_{\text{E}}^j(1-X_{\text{loss}}^j)$, where $X_{\text{loss}}^i$ and $X_{\text{loss}}^j$ are the atmospheric mass-loss fractions of the target and the impactor, respectively. For each family of collisions (i.e. for all the isolated collisions between a planet and each of the other less massive planets in the same disc), we kept the most destructive one, which removes the greatest amount of atmosphere.

The main caveat of our simple collision model is clearly the lack of N-body simulations since we can only handle one planet per simulation. Our approach to calculating potential collisions post-gas stage with other (less massive) planets formed under the same initial disc conditions in different simulations does not guarantee their simultaneous formation in a single simulation. On the other hand, considering only one collision per planet, particularly the most destructive one (the one that generates the maximum envelope mass-loss), could reduce the mass of the envelope between 16%-100%, with a mean value of 72%, as stated in Appendix D of V20. This consideration could be overestimating the envelope mass-loss rate since Ogihara & Hori (2020) showed that their mass-loss fraction ranges from a few percent to approximately 90%, but with a typical value around 20%. However, they also confirm that this percentage is higher if head-on collisions take place. On another note, Matsumoto et al. (2021), who conducted N-body simulations to trace planet sizes during the giant impact phase with envelope stripping through impact shocks, considering empirical envelope-loss rates obtained from smoothed particle hydrodynamics simulations by Kegerreis et al. (2020). These authors also mentioned that head-on collisions can lead to increased rates of envelope mass loss. Additionally, their findings indicate that protoplanets in inner orbits undergo a higher frequency of collisions, potentially leading them to become bare cores after giant impacts. This is an outcome we also found within our simulations. Thus, our model captures the overall outcomes observed in more detailed studies, despite its simplicity and reduced precision in modelling collisions.

A.4. Evolution model

After disc dispersal and the removal of mass from giant impact (for the ’collision case’), we computed the cooling and contraction of the planets’ envelopes, including the effect of mass-loss due to photoevaporation during 2 Gyr with the code COMPLETO (Mordasini 2020). We refer to this stage as the ’evolution’ phase. The choice of 2 Gyr is to ensure that all the stars considered in this study are still in the main sequence. The code solves the 1D equations for mass, momentum, and energy conservation. When calculating the luminosity evolution, the gravitational and internal energy of the core and envelope, as well as radiogenic heating were considered (Linder et al. 2019). The stellar X-ray flux evolution was taken from the model of McDonald et al. (2019). The envelopes are assumed to be composed of H, He, and H₂O (with the compounds uniformly mixed) and with the initial amounts stemming from the accretion process. The equation of state (EOS) of the H/He component corresponds to Chabrier et al. (2019) and for water, we use the AQUA EOS from Haldemann et al. (2020). The EOS for the iron-silicate core is the modified polytropic EOS of (Seager et al. 2007), which does not yield very precise planetary radii for sub-Earth size planets. For this reason, we re-computed the final radius at 2 Gyr with the more-up-to-date interior structure model of BICEPS (see next section). The atmospheric escape rates are taken from (Kubyshkina & Fossati 2021). The atmospheric escape follows a dependency of ${\text{Z}}_{\text{env}}^{-0.77}$ (Owen & Jackson 2012; Mordasini 2020), where Z_env is the water mass fraction of the envelope, namely, ${\text{Z}}_{\text{env}}=\frac{{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}}{{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}+{\text{M}}_{\text{HHe}}}$. We caution that the evaporation rates of envelopes enriched in heavy elements are still uncertain. It is imperative that future hydrodynamic escape calculations address this aspect.

A.5. Transit radius calculation and internal structure

At the end of the planetary evolution stage, for each modelled planet, we calculated its corresponding transit radius. For that purpose, we used the internal structure model of BICEPS by Haldemann et al. (2024). BICEPS solves the internal structure equations and calculates the planetary transit radius for a large range of planetary compositions. For this work, we assumed that all planets are structured in the following way: i) in the centre there is an iron core made out of pure iron, ii) surrounded by a mantle of rocky material (MgSiO₃), and iii) the outermost layer contains the volatile elements H, He, and H₂O. The composition of the volatile layer is given by M_HHe and M_H2O, which are provided for each planet by the evolution model. As in V20, the water is uniformly mixed with the H, and He; thus, Z_env is uniform throughout the volatile layer. For planets that lost all their volatiles during the evolutionary stage or did not accrete any volatiles, the structure model starts the integration at the mantle layer. We note that the mass ratio between iron core and rocky mantle is kept at Earth-like values of 1:2.

Depending on the planet’s internal luminosity, the orbital distance to its host star and the host stars effective temperature the thermal structure of the volatile layer will be different. To account for this effect, BICEPS uses the non-grey atmosphere model of Parmentier et al. (2015) to calculate the temperature profile of the outermost volatile layer. We further use the Schwarzschild criterion to determine if at a certain depth the volatile layer is stable against convection, namely, if the temperature gradient is radiative or adiabatic (Kippenhahn et al. 2012). At last the transit radius is calculated by determining the radius at the chord optical depth of τ_chord = 2/3 (Guillot 2010; Parmentier & Guillot 2014). The opacities used in this calculation are taken from Freedman et al. (2014). For more details on the internal structure model (especially the applied EOSs), we refer to Haldemann et al. (2024).

A.6. Model limitations

Several simplifications affect the results of our formation-evolution-structure calculations.

1. Solar elemental abundances:

   • (a)

     The rocky component of the planets is always assumed to be Earth-like in composition, namely, with 33% iron and 67% silicates by mass. In reality, a spread in iron and silicate abundances is expected for stars that have elemental abundances different than solar. The minimum iron mass fraction expected for stars in the thin disc (where the Sun and most planet-hosting stars are) is of ∼20%, (Michel et al. 2020, their Fig. A.1).

   • (b)

     The ice-to-rock ratio of 1:1 just beyond the ice line also stems from solar abundances. Stars with other elemental abundances are expected to exhibit different ratios. For thin-disc stars, the minimum water mass fraction of solid material beyond the ice line is of ∼30% (Michel et al. 2020, see their Fig. A.2).

   These two compositional assumptions imply that the "depleted diagonal band" found in the MR diagram of Fig.1 should still exist, but could be narrower, between lines of 20% iron mass fraction for rocky composition and of 30% water mass fraction. This would also translate in a persistent density valley for all stellar types between the same compositional lines.

2. Only one planetary seed per disc: N-body interactions are not yet included in PLANETALP. This implies that the planet typically reaches the disc inner edge when migrating inwards, and the possibility of being trapped in resonant chains with other planets at different semi-major axes cannot be modelled. This model limitation has a stark impact on the planets’ final orbital period, which should not be directly compared to observations.

3. Planetary seeds are all placed at the beginning of the simulations (time = 0). This also affects the final orbital period of wet versus rocky planets. Formation timescales are shorter for wet planets (as we showed in V20) and, as a consequence, they reach the inner disc at early times, when the zero torque location (the planet migration trap) is located near the inner edge of the disc. As time evolves, this migration trap location moves slightly outwards and dry planets tend to be trapped at larger periods. This is why in Fig. 3, we can see that most of the wet planets tend to be located at orbital periods of ∼12 days, while the dry ones are mostly located at orbital periods of ≳13 days.

4. Planetary cores are assumed to grow only due to pebbles. In reality, pebbles are converting into planetesimals in the disc (Drążkowska & Alibert 2017; Lenz et al. 2019), thus planetesimals could also contribute to the core growth. Except for very few cases (Alibert et al. 2018; Guilera et al. 2020; Kessler & Alibert 2023), formation models tend to assume only one of the two solid accretion types due to the uncertain pebble-to-planetesimal ratios along the disc. Planetesimal accretion would yield more mixed compositions (Brügger et al. 2020; Burn et al. 2021), blurring the density valley to some extent.

5. The atmospheric escape rates for atmospheres enriched in water are taken from empirical photoevaporation laws from protoplanetary discs (Ercolano & Clarke 2010) (Sect.A.4). Detailed hydrodynamical atmospheric escape simulations for mixtures of H, He and H₂O should be developed in the future.

6. Simplifications in the internal structure calculations: our upper irradiated atmospheres assume a solar composition (Parmentier et al. 2015). In addition, the gas opacities in the radiative zones are calculated with the model of Freedman et al. (2014), which accounts for certain degree of envelope enrichment, but are not specific to water.

Appendix B: Extended results

B.1. Bare cores after disc dissipation

The effect found in V20 is in principle expected to occur for any stellar mass because it is related to the change of the pebbles’ properties at the water ice line. This, in turns, produces more massive icy cores than rocky ones and this effects the core sizes as well.

In Fig. B.1, we show the histograms of the core sizes of the simulated planets at the end of formation. The blue bars indicate a core water mass fraction larger than f_H2O, core > 45%, the red ones, dry cores (f_H2O, core < 5%), and the green bars intermediate, mixed compositions (5%≤f_H2O, core ≤ 45%). We note that, as expected, icy cores tend to be bigger than rocky cores, for all stellar masses. A deficit of cores with sizes in the range of ∼1.5 − 2 R_⊕  is clear for M_⋆ ≥ 0.4 M_⊙. The valley separates dry from wet cores for M_⋆ ≥ 0.4 M_⊙. Interestingly, the lower the stellar mass, the larger the overlap in size between the dry and wet populations, and the larger the number of planets with intermediate compositions. In particular, for M_⋆ = 0.1 M_⊙ there is no second peak, as many icy cores are very small. This behaviour is due to migration. Type I migration occurs for lower planet mass at lower stellar masses (Paardekooper et al. 2011). This allows for smaller icy planets to reach the inner system for the low stellar mass cases. This feature was also reported by Burn et al. (2021).

Fig. B.1. Histograms of core sizes at the end of formation, for the different stellar masses, and for planets with final orbital period below 100 days. Red bars indicate fH2O, core < 5%, green bars 5%≤fH2O, core < 45%, and blue bars fH2O, core ≥ 45%, where fH2O, core is the core water mass fraction at the end of formation. The black lines show the overall core size distribution.

B.2. After 2 Gyr of evolution

Figure B.2 shows the histograms of the planet radii at 2 Gyr of evolution, for both the nominal and collisional cases. In agreement with Fig. B.1, the radius valley is non-existent for M_⋆ = 0.1 M_⊙, and starts to be visible for M_⋆ ≥ 0.4 M_⊙ for the collisional case (which resembles, rather, the bare core cases of Fig. B.1 since more atmospheric mass is removed when we allow for giant impacts). Interestingly, our model predicts that the case of M_⋆ = 0.7 M_⊙ is the one for which the radius valley is the most prominent, with a clear peak of water-rich mini-Neptunes for both nominal and collisional cases. This is in line with the findings of Kunimoto & Matthews (2020) which report the occurrence rate of mini-Neptunes increasing for K dwarfs compared to G dwarfs. It is important to clarify that these histograms should not be taken as absolute occurrence rates to compare directly with observations. Indeed, this study is not a population synthesis, in the sense that the initial conditions are not taken with the weights of the observed distributions of disc properties. We do consider the ranges of possible initial conditions stemming from observations, but not the weights of the distributions. Nevertheless, the fractions of certain types of planets (e.g. mini-Neptunes) can be compared relatively, between our different stellar masses. Our planet formation and evolution parameter study predicts that mini-Neptunes reach their peak of occurrence for stellar masses of M_⋆ ∼ 0.7 M_⊙ (strictly speaking, 0.4 < M_⋆ < 1.0 M_⊙ due to our binning in stellar mass). For the nominal cases of M_⋆ ≳ 1.0 M_⊙, the water-rich planets are no longer concentrated around a peak value, but they are, rather, spread over a large range of sizes. This is mainly the effect of gas accretion onto cores that are more massive compared to M_⋆ ≲ 1.0 M_⊙ and which can bind larger amounts of gas. When gas is removed by collisions, the mini-Neptune peak re-emerges, similarly to what we found in V20. Another effect that removes mini-Neptunes for the case of the most massive stars (M_⋆ ≳ 1.3 M_⊙) is the strong photoevaporation, which completely removes the H/He/H₂O envelopes for more planets compared to M_⋆ = 1.0 M_⊙. Hence, planets that were born beyond the ice line as ice-rich, can become bare rocky cores by evolution. In our simulations, this happens for 28% of the original water-rich mini-Neptunes for M_⋆ = 1.5 M_⊙ and for 14% of the original water-rich mini-Neptunes for M_⋆ = 1.3 M_⊙.

Fig. B.2. Histograms of planets’ radii after 2 Gyr of evolution, for different stellar masses, and for planets with final orbital period below 100 days. Red bars indicate $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }<5\%$, green bars $5\%\le f_{{\text{H}}_2\text{O},\text{core}}^{\prime }<45\%$, and blue bars $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }\ge 45\%$, where $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }$ is the core water mass fraction at 2 Gyr of evolution. The black lines show the overall planet size distribution. Left panels display results for the nominal set-up, and right panels for the collisional set-up (see Sect.A.3).

B.3. Water-worlds around low-mass stars migrated from beyond the ice line

The ice line of protoplanetary discs around low-mass stars is much closer to the central star than for more massive stars. For example, while for M_⋆ = 1.0 M_⊙, the ice line typically locates initially at ∼2-3 au, for M_⋆ = 0.1 M_⊙ it is at∼0.5 au at the beginning of the simulations. Because of this, the existence of water worlds around low-mass stars could in principle be explained without type I migration, in a scenario where wet-planets with orbital periods below 100 days simply formed in situ, with the ice line closer in. However, this never happens in our simulations, water worlds always accrete the bulk of their ices beyond the ice line and migrate inwards afterwards. To illustrate this, we show in the top panel of Fig. B.3 that the ice lines in discs around M_⋆ = 0.1 M_⊙ never manage to cross the 100-day orbit threshold (marked with the horizontal grey dashed line at ∼0.2 au) before photoevaporation opens a gap in the disc (time from which the solid lines become dashed lines). The bottom panel of Fig. B.3 shows that planets forming within the ice line for M_⋆ = 0.1 M_⊙ are always dry (planetary tracks represented with black curves only), while those that start their growth outside the ice line, but end inside it are water-rich (planetary tracks represented by turquoise curves that then turn into black once they cross the ice line) and can only reach the regions of orbital period below 100 days (marked by the dashed vertical line) through type I migration. We have checked that this scenario always takes place with all the icy planets around all the stars considered. For simplicity we only show it for some of them, selected randomly, around 0.1 M_⊙.

Fig. B.3. Time evolution of the ice lines for the discs around 0.1M⊙ (top) for the different protoplanetary discs considered in our simulations (see Sect.A.1, for α = 10−4). The ice lines turn from solid to dashed lines when the mid-gap opens in the protoplanetary disc (i.e. when the disc dissipates at the ice line’s location). The horizontal grey dashed line represents the orbit of ∼100 days. Growth tracks of planets forming around 0.1M⊙ that end within the 100-day period (bottom). Black lines indicate that the planet is located inside the ice line, while torquoise lines indicate the planet is beyond the iceline. Clearly, all planets that had accreted some ice started to form beyond a 100-day period.

B.4. Steam atmospheres

An interesting feature of Fig.1 is that for M_⋆ ≥ 0.7 M_⊙, most of the water-rich, H/He-poor planets (f_H2O > 0.45 and f_HHe < 10⁻², see Figs.1 and B.4) lie above the condensed-water line. To understand what is is behind the ’puffing up’ of the water-rich planets, we first need to isolate the effect of H/He from the effect of temperature on the planetary radii. For this, we plot in Fig.B.5 the M-R diagram as a function of equilibrium temperature for M_⋆ = 0.4.−0.7 M_⊙. The planets are also distinguished by a H/He mass fraction smaller or larger than 0.01 (we refer to them as H/He poor if f_HHe < 10⁻² and H/He-rich otherwise). Three things are clear from this diagram:

1. Water-rich and H/He-poor planets with the lowest equilibrium temperature (below 400 K, the ones around M_⋆ = 0.4 M_⊙) lie on the condensed-water line. These planets are cold enough for water to be in condensed form.

2. All of the water-rich, H/He-poor planets with T_eq ≥ 400 K lie on a different line (light blue in the figure). They are well described as R_P(M_P) = $a{M}_{\mathrm{P}}^b$, with a = 1.64183 and b = 0.180572 (see Sect.B.5).

3. All the planets with f_HHe ≳ 10⁻² are scattered either above the light-blue M-R curve for M_⋆ = 0.7 M_⊙ or above the condensed-water line for M_⋆ = 0.4 M_⊙.

Fig. B.4. Same as Fig.1 but the colour bar represents the planets’ mass fraction of H/He.

Fig. B.5. MR diagram of water-rich planets (fH2O ≥ 0.45) around M⋆ =  0.4 M⊙ (circles) and M⋆ =  0.7 M⊙ (diamonds). Black-bordered symbols indicate fHHe ≤ 10−2, while grey, fHHe > 10−2. The colour bar is the equilibrium temperature of the planets at zero albedo. The dashed-blue line shows the condensed-water line, the light-blue dashed a fitting MR relation to the steam worlds, described in Sect.B.5.

Indeed, for all the planets with f_HHe > 10⁻², the radius is strongly affected by the amount of primordial gas and by the equilibrium temperature. On the other hand, the planets that follow the light-blue line are H/He poor but their full atmospheres are hot enough for water to be present as vapour. This is better illustrated in Fig.C.2 where we show specific pressure-temperature atmospheric profiles (Appendix C.2).

B.5. MR relations: Analytical fits

In this section, we specify the mass-radius (MR) analytical curves deployed in Fig.1. Such MR relations were obtained by fitting our results for the rocky, condensed-water, and steam planets computed with BICEPS (Sect.A.5.)

For rocky planets, defined as those planets with final water contents lower than 5% in mass and H-He envelopes lower than 10⁻⁶ with respect to the total mass, we find the following relation:

$$\begin{array}{c}\hfill \frac{{\mathrm{R}}_{\mathrm{P}}}{{\mathrm{R}}_{\oplus }}=a_{\mathrm{rocky}}{\left(\frac{{\mathrm{M}}_{\mathrm{P}}}{{\mathrm{M}}_{\oplus }}\right)}^{b_{\mathrm{rocky}}},\end{array}$$(B.1)

with a_rocky = 0.999009 ± 0.0002802 and b_rocky = 0.279514 ± 0.0002377. Our fit is in very good agreement with the one found by Zeng et al. (2019), who found a_rocky = 1 and b_rocky = 1/3.7 ≈ 0.27027. We note that we do not have any bare rocky planet with mass larger than 10 M_⊕, reason for which this fit is valid for planetary masses below that.

For condensed-water planets, which are planets with a final water content greater than 45% in terms of mass (and less or equal than 50% due to the assumed ice-to-rock ratio at the ice line), H-He envelopes lower than 10⁻² by mass, and T_eq < 400 K, we find:

$$\begin{array}{c}\hfill \frac{{\mathrm{R}}_{\mathrm{P}}}{{\mathrm{R}}_{\oplus }}=a_{\mathrm{cw}}{\left(\frac{{\mathrm{M}}_{\mathrm{P}}}{{\mathrm{M}}_{\oplus }}\right)}^{b_{\mathrm{cw}}},\end{array}$$(B.2)

with a_cw = 1.24191 ± 0.001667, and b_cw = 0.267404 ± 0.0008558. The fit was performed for planets with 0.1 ≤ M_P ≤ 20 M_⊕ and therefore should only be used in this mass range. We note again that our fit is in very good agreement with the one found by Zeng et al. (2019) for those planets with 50% of water by mass, who found a_cw = 1.24, and b_cw = 1/3.7 ≈ 0.27027.

Finally, for the steam-worlds, defined as the condensed-water worlds except that T_eq  ≥  400 K, we obtain the following fit:

$$\begin{array}{c}\hfill \frac{{\mathrm{R}}_{\mathrm{P}}}{{\mathrm{R}}_{\oplus }}=a_{\mathrm{steam}}{\left(\frac{{\mathrm{M}}_{\mathrm{P}}}{{\mathrm{M}}_{\oplus }}\right)}^{b_{\mathrm{steam}}},\end{array}$$(B.3)

with a_steam = 1.64183 ± 0.008192, and b_steam = 0.180572 ± 0.001608. The fit was performed for planets with 1 ≤ M_P ≤ 40 M_⊕ and should only be used in this mass range. We note that our calculations indicate an increase of planet radius of 30%, 25%, and 10% between steam worlds and condensed-water worlds of masses of 1, 2, and 5 M_⊕, respectively (for planets with water mass fractions of 50%). As an interesting final remark, we note that the MR relation we find for steam worlds does not depend strongly on equilibrium temperature for different stellar masses.

B.6. Mean density: Histograms

Figure B.6 shows the histograms of the planetary mean densities for the different stellar masses. The left panels are for the planets around the M and K dwarfs, while the right panel for planets orbiting G-F and A stars. In the top panels, the density refers to the ‘normalised density’, namely, the mean density divided by the density that the planet would have if it had an Earth-like composition (i.e. a mass radius-relation given by Eq. B.1). The bottom panels display the histograms of the physical density in cgs. Clearly, the density valley is better visualised in terms of the normalised density, especially for the low-mass stars. The peak of planets at 0.5-0.6 ρ/ρ_⊕, S of the left-top panel corresponds to the condensed-water worlds, which do not exist in our results for high-mass stars due to the higher equilibrium temperatures (see Fig. 4). In terms of physical density (bottom panels), the valley occurs at ρ≈ 4-5 g/cm³ for M and K dwarfs, and for ρ≈ 3-4 g/cm³ for G-F-A spectral types.

Fig. B.6. Planetary mean density for the different stellar masses. Left (right) panels show values for the low- (high-)mass stars, with the each colour corresponding to one stellar mass as indicated in the legend. The top panels show the normalised density (see text), while the bottom panels show the density in cgs.

B.7. A few dry mini-Neptunes in the MR diagram

Figure 1 contains a few dry ‘outliers’ in the region of the MR diagram above the condensed-water line, where most planets are water-rich. These planets were actually formed beyond the ice line, just like the wet planets that surround them. The difference with the latter is that they lost most of their water (and H/He) due to photoevaporation during the evolutionary phase. To illustrate this, we show in Table B.1, the initial and final planet mass, water mass, and H/He mass for two planets that end up very close in the MR diagram of M_⋆ = 1.5 M_⊙. We also tabulate the rock mass (invariable during the evolution), and the initial and final water mass fraction (f_H2O), H/He mass fraction (f_HHe) and envelope metallicity (Z_env, defined as the mass of water divided by the mass of water plus H/He). One of the planets ends up water-rich and

the other water-poor after 2 Gyr of evolution (referred as ‘wet’ and ‘dry’ planets in the table, respectively). The reason behind the different evolutionary paths is the atmospheric escape rate going as $\sim {\text{Z}}_{\text{env}}^{-0.77}$ (with Z_env being the envelope water mass fraction, which remains constant throughout the evolution by hypothesis). This means that the smaller the initial Z_env, the faster the atmospheric mass loss. Indeed, the two planets have very different initial Z_env, despite of having similar initial f_H2O. The smaller Z_env of the dry planet is a consequence of having a much larger initial amount of H/He. Because of the smaller Z_env of the dry planet, it loses 99.45% of its initial envelope during 2 Gyr of evolution, while the ‘wet’ planet loses only 6.46% of it. This leaves the ‘dry’ planet with only f_H2O = 0.55% for the bulk water fraction at the end of the evolutionary phase, while the ‘wet’ planet retains f_H2O = 47.6%. Despite being highly ‘desiccated’, the dry planet would still show a high water signature in its atmospheric spectra, due to it retaining an atmospheric water mass fraction (or envelope metallicity) of Z_env ≈ 55%. We caution, however, that atmospheric escape rates of envelopes enriched with elements heavier that H/He are highly uncertain and that more theoretical research is needed to provide meaningful predictions of water mass fraction during the evolutionary phase.

Appendix C: Extended discussions

C.1. Dependence of the radius valley location on the formation model parameters

Planet formation models contain a set of free parameters, with a range of possible values inferred from either observations or lab experiments. These parameters can affect the outcome of planet formation simulations in regard to the predicted planet mass, planet radius, and orbital distance to the star, among others. The prime parameters that can result in very different outcomes in our pebble-based formation model are:

1. the α-viscosity parameter

2. the fragmentation threshold velocity of pebbles

3. the opacity of the gaseous envelope

While a deep analysis of the impact of these parameters on our results is beyond the scope of the present work, it is nevertheless pertinent to discuss the effect they could have on the location of the radius valley.

C.1.1. Dependence on the α-viscosity parameter

The α-viscosity parameter affects the growth of the planetary cores in two ways. On the one hand, it affects the pebbles’ sizes (the higher the α, the more turbulent is the disc, hence the more the pebbles collide and break, maintaining the pebbles at smaller sizes). The larger the pebbles, the faster the core grows (e.g. Lambrechts et al. 2014; Venturini et al. 2020b). On the other hand, the disc viscosity impacts the disc’s structure, including the disc aspect ratio. This means that the pebble isolation mass (the maximum mass a core can acquire by pebble accretion) is indirectly affected by α, because the pebble isolation mass depends strongly on the disc aspect ratio. In principle, the larger the α, the larger the pebble isolation mass, because the larger the viscosity, the more massive the planet needs to be to open a partial gap. We checked that the pebble isolation mass at the ice line is of about 8 M_⊕ for α = 10⁻⁵, compared to ∼15 M_⊕ for α = 10⁻⁴ for a standard disc around a Sun-like star. A planet of 8 M_⊕ with f_H2O = 0.5 and without H-He would have a size of 2.86 R_⊕ (from using the MR relation of Eq.B.3), which means that the mean of water-worlds forming at such low viscosity are still expected to contribute mainly to the peak of mini-Neptunes and not to fill the radius valley (at ≈1.8 − 1.9 R_⊕ for M_⋆ = 1.0 M_⊙).

In Venturini et al. (2020b), we discussed extensively the impact of α on the planetary growth by pebble accretion within the ice line, considering values of α = 10⁻⁵, 10⁻⁴and10⁻³ (which are the values inferred from disc observations; see e.g. Rosotti 2023). We found that for α = 10⁻³, the rocky pebbles were so small (∼10 μm) that planets were basically not growing inside the ice line. On the other extreme, for α = 10⁻⁵, we found that the pebbles reach much bigger sizes (∼1 cm) in the dry inner regions, and planetary growth proceeds extremely fast (in the order of 10⁴ years). We also found that for the formation of rocky cores, considering α = 10⁻⁵ increases the maximum bare rocky core mass from 5 M_⊕ to 7 M_⊕ with respect to α = 10⁻⁴. This would still yield a maximum radius of bare rocky cores of 1.69-1.95 R_⊕ (from using the MR relation of (Zeng et al. 2019) and of Eq.B.1), that is, at the location of the valley.

It is also important to mention that the planet migration prescriptions adopted in our model are not really valid for α ∼ 10⁻⁵. For such small values of α, even a low-mass planet could open a partial gap. McNally et al. (2019) showed that in this case, vortices form at the edges of the partial gap, changing the migration regime. Thus, based on our previous results, in this work we decided to keep the same approach as in V20, using only α = 10⁻⁴ and α = 10⁻³. We analysed the effect of this two values of the α-viscosity parameter on the location of the radius valley on Fig.C.1. The plot shows the MR diagrams for the synthetic planets, distinguishing according to the cases where α = 10⁻⁴ (orange circles) or α = 10⁻³ (green triangles). Overall, we notice that the void of the ‘diagonal band’ delimited by the Earth-like composition and the condensed-water lines (which sets the location of the radius and density valley; see main text) is preserved for both values of α. For α = 10⁻³, we have very few rocky planets, which is expected due to the extremely low pebble sizes, resulting in negligible core growth inside the ice line (Venturini et al. 2020b). Overall, the jump in pebble size at the ice line is the underlying reason for the dichotomy in core size and core composition in our pebble-based formation model (Venturini et al. 2020a). This jump exists at any value of the α-viscosity parameter (Fig. 4, top panel, from Venturini et al. 2020b).

Fig. C.1. MR diagrams according to the adopted values of the α-viscosity parameter. Green triangles: α = 10−3. Orange circles:α = 10−4. Both nominal and collisional cases are displayed.

C.1.2. Dependence on the fragmentation threshold velocity of pebbles

Regarding the impact of the fragmentation threshold velocities on the results, in V20 we showed that for Sun-like stars, considering v_th= 1 m/s for silicates and v_th= 5 m/s for icy pebbles gives very similar results to our nominal assumption of v_th= 1 m/s for silicates and v_th= 10 m/s for icy pebbles. If instead v_th = 2 m/s is chosen for the icy pebbles, fewer planets with the maximum core ice mass fraction of 0.5 form and more planets with intermediate ice compositions emerge. This would contribute to the filling of the radius valley. Finally, we tested as well a very low fragmentation threshold velocity of v_th= 1 m/s, constant along the disc. In that case, the pebbles remain so small along the disc, that planet formation by pebble accretion is very inefficient: no planet forms with mass above Earth and no planet migrates within 100 days of orbital period.

The values of the fragmentation threshold velocity suited for protoplanetary disc conditions have been a topic of intense debate in the past years. The values adopted in this work (10 m/s beyond the water-ice line and 1m/s within the ice line), are the standard values used in dust growth and evolution studies (Drążkowska et al. 2016; Drążkowska & Alibert 2017), and on planet formation models, (e.g. Drążkowska et al. 2021; Schneider & Bitsch 2021; Savvidou & Bitsch 2023). These values are also supported by many experimental works (e.g. Gundlach & Blum 2015; Musiolik et al. 2016). While recent laboratory experiments found that amorphous water-ice particles could be as fragile as silicates (e.g. Gundlach et al. 2018; Musiolik & Wurm 2019), more recent works, such as those of Nietiadi et al. (2020) and Musiolik (2021), revealed that ices of various volatile molecules on the surface of dust aggregates act like liquids during collisions. This behaviour increases the sticking properties and enables higher velocity collisions that do not lead to the fragmentation of ice-coated dust aggregates. Thus, given the current knowledge on the topic, we consider that a value of 10 m/s for the fragmentation threshold velocity of ice-rich grains is a sufficiently reasonable assumption. In addition, fragmentation threshold velocities of 1 m/s for silicate and 10 m/s for icy pebbles were recently adopted by Cañas et al. (2024) to explain the density dichotomy in the Kuiper Belt Objects. For our results, the location of the radius valley is invariant as long as the fragmentation threshold velocity of icy pebbles is in the range of 5-10 m/s.

C.1.3. Dependence on the envelope’s opacity

The opacities in the gaseous envelope control its cooling and contraction and, hence, the gas accretion rate as long as the planet’s envelope connects smoothly with the gaseous disc (attached phase). The largest uncertainty in the opacity values stem from the poorly constrained grain sizes within a planetary envelope. In principle, grains can settle and grow compared to ISM sizes, which goes in the direction of reducing the dust opacities compared to the ISM. However, the level of turbulence in planetary envelopes is hard to predict, the growth of grains depends upon their size and structure (e.g. fractal versus non-fractal, which is unknown), as well on the magnitude of the envelope ‘recycling’. Indeed, recent 3D hydrodynamical simulations have shown that as long as a planet is located at short orbital distances (within approximately 1 au from the central star Ormel et al. 2015; Moldenhauer et al. 2022; Wang et al. 2023), the gas that is accreted by the planet can flow back to the protoplanetary disc. This means that for planets building the envelopes at short orbital distances (as in our scenario, where planets migrate fast towards the disc inner edge and continue accreting gas at the stranded locations), the small dust grains would be constantly replenishing the envelopes as fresh gas flows in. This is why we chose the high grain opacity values from Bell & Lin (1994), which correspond to the ISM size grains (and whose values vary with pressure and temperature).

In any case, due to the uncertainty on the dust opacity values of the planetary envelopes during planet formation, it is worth discussing its impact on the location of the radius valley. The study of Mordasini (2020) studies the emergence of the radius valley taking initial conditions from population synthesis studies that assume (contrary to our assumptions) a very low dust opacity during the formation phase (Bell & Lin 1994, reduced by a factor of 0.003). The location of the valley, corresponding to the maximum mass of bare rocky cores, is at 1.67 R_⊕ at 0.1 au for a Sun-like star Mordasini (2020, as seen in his Fig.4). We obtained the exact same value in our simulations with high dust opacity (nominal case of M_⋆ = 1.0 M_⊙ displayed in Fig.B.2). Indeed, Mordasini (2020) discussed the practically negligible effect of the opacities on the radius valley location in Sect.2.3.7 of that paper. This is because at the end of disc lifetime, the lower the envelope dust opacity, the more massive the gaseous envelope for a given core mass. A more massive envelope is more extended than a thin one, which increases the atmospheric escape rate (which depends on the planet radius to the cube Owen & Wu 2017; Mordasini 2020). Thus, even if the planet started with a thicker envelope, photoevaporation would remove it completely for the ranges of core masses we are discussing here (below approximately 10 M_⊕). We found this exact same behaviour in our study of Venturini et al. (2020b) when we tested a set of simulations with 100× reduced Bell & Lin (1994) opacities, namely: the location of the radius valley was unchanged.

C.2. Condensed-water or steam worlds: Dependence on equilibrium temperature and model assumptions

It is well known that the atmospheric composition, the distribution of heavy elements and the choice of irradiation model affect the calculation of the planetary radii (Guillot 2010). In V20, we argued that assuming the water as mixed with H/He was the most physically motivated assumption since for Sun-like stars all our synthetic planets had temperatures high enough for water to be in the form of steam throughout all the atmospheres. In this work, we made the same assumption of water mixed with H/He during the evolution simulations and the computation of the planetary radii. However, we noted that our synthetic planets can now have temperatures low enough for water to condense, particularly around M dwarfs (our M_⋆ = 0.1, 0.4 M_⊙ cases, Appendix B.4). The fate of condensed water in H₂O-H-He atmosphere is not trivial to predict. In principle, condensed water on the top layers would settle into deeper regions and if there is atmospheric circulation, would rise up again. Thus, the assumption of uniform envelope metallicity might still hold. For the case of pure water envelopes, a better physical treatment would be to allow water to sink to the next layer, as long as the condensation conditions in pressure and temperature hold. Unfortunately, our structure code cannot handle yet compositional gradients properly. However, we caution that changing the distribution of water will probably change the value of planet radius in cases where water condensation takes place. Not only does the distribution of water affect the resulting transiting radius, but also other model assumptions, such as the choice of opacities and the irradiation model. We note that, compared to our results, the models of Aguichine et al. (2021) predicted a substantially larger planet radius for planets composed of half-rocks-half-water by mass when T_eq ≥ 400K (their Fig. 5). The key difference with our model is their assumption of the core-envelope boundary to be extremely hot (in a region of the phase diagram of water where water is in ’supercritical’ state). This assumption is based on evolution calculations accounting for the ‘runaway greenhouse effect of water’ (e.g. Turbet et al. 2019). More recent work points out that such effect might have been overestimated in the past due to the assumption of purely convective envelopes (Selsis et al. 2023).

In our calculations, we did not impose any condition on the pressure-temperature of the envelope-core boundary, we simply integrated the structure equations inwards from the top of the atmosphere (with the outer boundary conditions depending on the equilibrium temperature and on the luminosity that the planet irradiates according to its cooling history; see Appendix A.4). The types of P-T profiles that we find for planets with pure water envelopes are illustrated in Fig.C.2. In the figure, we show the profiles of two planets, one of the two has a fully condensed water envelope while the other a full steam one. To better capture the transition between the two, we ran the case of the steam world again, but reducing its equilibrium temperature (red and green curves of Fig.C.2). We note that the transition between the condensed and steam worlds happens for 383 < T_eq < 396K; in our model, this translates into a temperature at the top of the atmosphere of 268 < T_out < 276 K, respectively. The reason for T_out being smaller that T_eq by ∼80 − 120 K is the non-grey atmospheric model of Parmentier et al. (2015). The previous grey atmospheric model of Guillot (2010) also yields T_out smaller than T_eq, but only by ∼20 K. The nearly 100 K difference in T_out between the two atmospheric models makes the water of the planets orbiting M_⋆ = 0.4 M_⊙ take a condensed form when using the non-grey irradiated model of Parmentier et al. (2015) (used in COMPLETO) – as compared to the vapour form when using the one of Guillot (2010) (used in COMPLETO). This is shown in Fig. C.3. To summarise, the radii of water-worlds near the condensation line are very sensitive to model assumptions. More efforts are needed in modelling the atmospheres of water worlds to compute accurate transit radii.

Fig. C.2. Pressure-temperature profiles of two synthetic mini-Neptunes with fH2O = 0.5 and equilibrium temperatures of 466 K (planet 503, blue line) and 326 K (planet 506, orange line). The black solid lines in the background are the phase boundaries of pure H2O given by the AQUA EOS (Haldemann et al. 2020). Both planets have a similar mass of approx. 1.5 M⊕, but different sizes due to the different physical state of water (planet 503 is a steam-world, with RP = 1.690 R⊕ , and planet 506 is a condensed-water world with RP = 1.657 R⊕ ). The red and green curves correspond to planet 503, but the equilibrium temperature was reduced to see where the transition from steam to condensed-water world occurs.

Fig. C.3. MR diagrams comparing planetary radii obtained with BICEPS (purple triangles) versus COMPLETO (orange circles). For both set-ups water is assumed to be mixed uniformly with the H/He.

C.3. Locating the steam worlds

A controversial aspect of the Luque & Pallé (2022) study is the conclusion that the planets that fall on the condensed-water line in the MR diagram are indeed 50% rock+50% condensed water by mass, because most of those exoplanets are too hot for water to be condensed. Our models suggest that an upper equilibrium temperature for the transition would be 400 K, and other works point to even lower temperatures (e.g. Aguichine et al. 2021). The bulk of planets in the sample from Luque & Pallé (2022) have equilibrium temperatures higher than this. To analyse this in more detail, comparing the results of our calculations with observations, we plot on the top panels of Fig. C.4, all real exoplanets with T_eq ≤ 400 K for 0.25 ≤ M_⋆ < 0.55 M_⊙ (M dwarfs) and T_eq ≤ 600 K for 0.55 ≤ M_⋆ < 0.85 M_⊙ (K dwarfs).

Fig. C.4. Comparison of real exoplanets for M⋆ =  0.4 M⊙ (left panels) and M⋆ =  0.7 M⊙ (right panels). The limit in equilibrium temperature is for proper comparison with our simulations. The brown- and blue-dashed lines correspond to Earth-like and 50% water+50% Earth-like compositions, with dark-blue indicating condensed-water and light-blue, steam-worlds (as defined in Appendix B.5). The green-dotted curves correspond to 50% steam+50% Earth-like from Aguichine et al. (2021) (Teq = 400 K for dark-green and Teq = 600 K for light-green). The real exoplanets are those with a relative error in radius less than 25% and in mass of less than 70%, with the data taken form the NASA Exoplanet Archive on 08.09.23. For the bottom panels, the real exoplanets are shown in the background in grey.

For the M dwarfs, we note that only six real exoplanets have T_eq ≤ 400 K and fairly good mass and radius measurements. Three of those have mass and radius compatible with the condensed-water line (within errors). We note that for the planets that have sizes above the condensed-water line, we actually have synthetic planets with that mass and radius (Fig.C.4, bottom-left panel). These planets, despite of having water in condensed form, have some remnant H/He that increases the planet radii compared to the condensed-water line.

For the K dwarfs (Fig.C.4, right panels), the coldest real planet (T_eq ≈ 350 K) falls exactly on the condensed-water line, as our calculations predict. Another five real planets fall in between the steam- and condensed-water lines, where our model predicts that the atmospheres should be in the form of steam. Nevertheless, we do have synthetic planets in that region of the parameter space (bottom-right panel). These planets have actually steam atmospheres but the total water content is of f_H2O ≈ 20 − 40% (Fig.1). Thus, this could be a plausible explanation for the real exoplanets falling in between the steam- and condensed-water lines for T_eq ≳ 400 K (composition that could also explain the exoplanets analysed by Luque & Pallé 2022).

We also note that while our model yields some steam-worlds with masses between 1 and 3 M_⊕ and T_eq ≈ 500 − 600 K for M_⋆ = 0.7 M_⊙, no real exoplanets have yet been confirmed in that part of the parameter space (see top-right panel of Fig. C.4). Nevertheless, the synthetic planets falling in that part of the MR diagram represent only 5% of the synthetic water worlds produced for that stellar mass, so we expect those planets to be rare. In addition, we note that the number of sub-Neptunes around K dwarfs with T_eq ≤ 600 K is still low to draw any strong conclusions.

Another clear difference between the synthetic and real exoplanets observed in this figure is the lack of real rocky exoplanets. This is because of the imposed maximum equilibrium temperature. Real rocky exoplanets exist but typically at higher equilibrium temperatures. Improved modelling on the disc inner edge and the addition of N-body interactions and orbital evolution due to tides is needed to better reproduce the observed distribution of planetary periods.

All Tables

All Figures

Fig. 1. MR diagram for different stellar masses. The colour bar represents the total water mass fraction in the planets. Dots correspond to the nominal output, diamonds to collisions. The brown- and blue-dashed lines correspond to a theoretical composition of Earth-like and 50% water + 50% Earth-like, as defined in Appendix B.5, with dark-blue indicating condensed-water worlds and light-blue indicating steam worlds. The green-dotted curves correspond to 50% steam + 50% Earth-like from Aguichine et al. (2021) (Teq = 400 K for dark-green and Teq = 600 K for light-green). Grey dots are real exoplanets with mass and radius measurements better than 75 and 20%, respectively (data taken from the NASA Exoplanet Archive on 08.09.23). The real exoplanets in each panel correspond to the stellar masses of 0–0.25 M⊙ for M⋆ = 0.1 M⊙, 0.25–0.55 M⊙ for M⋆ = 0.4 M⊙, 0.55–0.85 M⊙ for M⋆ = 0.7 M⊙, 0.85–1.15 M⊙ for M⋆ = 1.0 M⊙, 1.15–1.4 M⊙ for M⋆ = 1.3 M⊙, and 1.4–1.7 M⊙ for M⋆ = 1.5 M⊙.

In the text

Fig. 2. Radius valley fit as a function of the stellar mass (blue line with errors in light-blue) for all our synthetic planets (colour-coded as a function of the water mass fraction). Dots correspond to the nominal output, diamonds to collisions. The green dashed line with the light-green shaded area show the best fit found by Petigura et al. (2022), the violet dashed line and the lilac shaded area represent the results found by Berger et al. (2020), and the mustard dashed line and the light-mustard shaded area represents the best fit by Luque & Pallé (2022). Our best fit was found using https://github.com/parkus/gapfit (Loyd et al. 2020).

In the text

Fig. 3. Radius valley fit as a function of orbital period (blue line with errors in light-blue). The colour bar indicates total water mass fraction. Top panel: Planets around KGF-type stars. The green line and shade show the best fit found by Petigura et al. (2022) for these stellar types. Bottom panel: Planets around M dwarfs. The green line and shade show the best fit found by Luque & Pallé (2022). The particular clustering in the orbital period in both panels is due to the model limitations (Appendix A.6).

In the text

Fig. 4. Density-mass diagram (density normalised by an Earth-like composition) as a function of equilibrium temperature (colour bar). Left panel: Planets orbiting stars of 0.1 M⊙ (circles), 0.4 M⊙ (squares), and 0.7 M⊙ (diamonds). Right panel: Planets orbiting stars of 1.0 M⊙ (circles), 1.3 M⊙ (squares), and 1.5 M⊙ (diamonds). Note: the colour bars span different ranges in the two panels. For both panels, black-bordered symbols indicate fHHe < 1%, dark-grey 1 ≤ fHHe ≤ 10%, and light-grey fHHe > 10%. The dashed-lines show the compositional curves of Earth-like (brown), pure-silicate (orange), condensed water (blue), and steam worlds (light-blue). The grey-dotted lines correspond to radii of 1 and 4 R⊕.

In the text

	Fig. B.1. Histograms of core sizes at the end of formation, for the different stellar masses, and for planets with final orbital period below 100 days. Red bars indicate fH2O, core < 5%, green bars 5%≤fH2O, core < 45%, and blue bars fH2O, core ≥ 45%, where fH2O, core is the core water mass fraction at the end of formation. The black lines show the overall core size distribution.
In the text

Fig. B.2. Histograms of planets’ radii after 2 Gyr of evolution, for different stellar masses, and for planets with final orbital period below 100 days. Red bars indicate $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }<5\%$, green bars $5\%\le f_{{\text{H}}_2\text{O},\text{core}}^{\prime }<45\%$, and blue bars $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }\ge 45\%$, where $f_{{\text{H}}_2\text{O},\text{core}}^{\prime }$ is the core water mass fraction at 2 Gyr of evolution. The black lines show the overall planet size distribution. Left panels display results for the nominal set-up, and right panels for the collisional set-up (see Sect.A.3).

In the text

Fig. B.3. Time evolution of the ice lines for the discs around 0.1M⊙ (top) for the different protoplanetary discs considered in our simulations (see Sect.A.1, for α = 10−4). The ice lines turn from solid to dashed lines when the mid-gap opens in the protoplanetary disc (i.e. when the disc dissipates at the ice line’s location). The horizontal grey dashed line represents the orbit of ∼100 days. Growth tracks of planets forming around 0.1M⊙ that end within the 100-day period (bottom). Black lines indicate that the planet is located inside the ice line, while torquoise lines indicate the planet is beyond the iceline. Clearly, all planets that had accreted some ice started to form beyond a 100-day period.

In the text

	Fig. B.4. Same as Fig.1 but the colour bar represents the planets’ mass fraction of H/He.
In the text

Fig. B.5. MR diagram of water-rich planets (fH2O ≥ 0.45) around M⋆ =  0.4 M⊙ (circles) and M⋆ =  0.7 M⊙ (diamonds). Black-bordered symbols indicate fHHe ≤ 10−2, while grey, fHHe > 10−2. The colour bar is the equilibrium temperature of the planets at zero albedo. The dashed-blue line shows the condensed-water line, the light-blue dashed a fitting MR relation to the steam worlds, described in Sect.B.5.

In the text

	Fig. B.6. Planetary mean density for the different stellar masses. Left (right) panels show values for the low- (high-)mass stars, with the each colour corresponding to one stellar mass as indicated in the legend. The top panels show the normalised density (see text), while the bottom panels show the density in cgs.
In the text

	Fig. C.1. MR diagrams according to the adopted values of the α-viscosity parameter. Green triangles: α = 10−3. Orange circles:α = 10−4. Both nominal and collisional cases are displayed.
In the text

Fig. C.2. Pressure-temperature profiles of two synthetic mini-Neptunes with fH2O = 0.5 and equilibrium temperatures of 466 K (planet 503, blue line) and 326 K (planet 506, orange line). The black solid lines in the background are the phase boundaries of pure H2O given by the AQUA EOS (Haldemann et al. 2020). Both planets have a similar mass of approx. 1.5 M⊕, but different sizes due to the different physical state of water (planet 503 is a steam-world, with RP = 1.690 R⊕ , and planet 506 is a condensed-water world with RP = 1.657 R⊕ ). The red and green curves correspond to planet 503, but the equilibrium temperature was reduced to see where the transition from steam to condensed-water world occurs.

In the text

	Fig. C.3. MR diagrams comparing planetary radii obtained with BICEPS (purple triangles) versus COMPLETO (orange circles). For both set-ups water is assumed to be mixed uniformly with the H/He.
In the text

Fig. C.4. Comparison of real exoplanets for M⋆ =  0.4 M⊙ (left panels) and M⋆ =  0.7 M⊙ (right panels). The limit in equilibrium temperature is for proper comparison with our simulations. The brown- and blue-dashed lines correspond to Earth-like and 50% water+50% Earth-like compositions, with dark-blue indicating condensed-water and light-blue, steam-worlds (as defined in Appendix B.5). The green-dotted curves correspond to 50% steam+50% Earth-like from Aguichine et al. (2021) (Teq = 400 K for dark-green and Teq = 600 K for light-green). The real exoplanets are those with a relative error in radius less than 25% and in mass of less than 70%, with the data taken form the NASA Exoplanet Archive on 08.09.23. For the bottom panels, the real exoplanets are shown in the background in grey.

In the text