MNRAS 000, 1–12 (2025) Preprint 24 March 2025 Compiled using MNRAS LATEX style file v3.

On the road to the radius valley: distinguishing between gas dwarfs and
water worlds with young transiting exoplanets
James G. Rogers1★
1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom

Accepted XXX. Received YYY; in original form ZZZ

arXiv:2503.17364v1 [astro-ph.EP] 21 Mar 2025

ABSTRACT
The detection of young transiting exoplanets represents a new frontier in our understanding of planet formation and evolution.
For the population of observed close-in sub-Neptunes, two proposed formation pathways can reproduce their observed masses
and radii at ∼ Gyr ages: the “gas dwarf” hypothesis and the “water world” hypothesis. We show that a sub-Neptune’s size at
early ages ≲ 100 Myrs is strongly dependent on the bulk mean molecular weight within its envelope. As a result, gas dwarfs
and water worlds should diverge in size at early ages since the mean molecular weight of gas dwarf envelopes is predicted to be
smaller than that of water worlds. We construct population models under both scenarios that reproduce Kepler demographics in
the age range ∼ 1 − 10 Gyrs. We find tentative evidence that the gas dwarf model is more consistent with the small population of
young exoplanets < 40 Myrs from TESS. We show that planet radius is relatively insensitive to planet mass for young, puffy sub-
Neptunes, meaning that well-characterised masses are not necessarily required to exploit the effects of mean molecular weight
at the population level. We confirm the predicted difference in planet size between the models is also true under mixed-envelope
scenarios, in which envelopes consist of mixtures of hydrogen and steam. We highlight that transit surveys of young exoplanets
should target the youngest observable stellar clusters to exploit the effects of mean molecular weight.
Key words: planets and satellites: atmospheres - planets and satellites: formation

1 INTRODUCTION dominated atmospheres accreted from the nascent protoplanetary

disc (e.g. Lee et al. 2014; Ginzburg et al. 2016). Then, through at-
Twenty years on from the first sub-Neptune discoveries (Santos et al.
mospheric escape processes, including boil-off (Owen & Wu 2016;
2004; Rivera et al. 2005), understanding their nature still remains a
Ginzburg et al. 2016), core-powered mass loss (Ginzburg et al. 2018)
mystery (e.g. Bean et al. 2021). With typical sizes between ∼ 2−4𝑅 ⊕ ,
and XUV photoevaporation (Owen & Wu 2013; Lopez & Fortney
and orbital periods ≲ 100 days, sub-Neptunes orbit ∼ 50% of Sun-
2013), planets that are smaller in mass and closer to their host stars
like stars (Petigura et al. 2013; Batalha et al. 2013; Fressin et al. 2013).
were stripped of their envelopes, leaving behind bare cores. These
Sub-Neptunes sit above the “radius valley”, an observed paucity of
rocky cores thus formed the population of super-Earths, whereas
planets at ∼ 1.8𝑅 ⊕ , although its exact location varies with orbital
the planets that could retain their H2 /He envelopes formed the sub-
period and stellar mass (e.g. Fulton et al. 2017; Petigura et al. 2018;
Neptune population. This evolutionary model is supported by di-
Van Eylen et al. 2018; Petigura et al. 2022; Ho & Van Eylen 2023).
rect observations of escaping hydrogen and helium (e.g. Dos Santos
Their counterparts, referred to as super-Earths, sit below the valley.
2023).
Crucially, mass measurements have shown that the bulk density for
planets above and below the valley are markedly different (e.g. Wu On the other hand, the “water world” hypothesis states that the
& Lithwick 2013; Weiss & Marcy 2014; Rogers 2015; Wolfgang radius valley represents the end-product of two distinct planet for-
et al. 2016; Chen & Kipping 2017). Super-Earths are consistent mation pathways (e.g. Mordasini et al. 2009; Venturini et al. 2016;
with a rocky, Earth-like composition. Conversely, sub-Neptunes have Raymond et al. 2018; Zeng et al. 2019; Mousis et al. 2020; Burn et al.
much lower bulk densities, suggesting the presence of low-density 2024a). Rocky super-Earths were formed in the inner regions of pro-
compositional components within these planets. The radius valley toplanetary discs, where the major solid components are silicates.
can thus also be interpreted as an approximate delineation between Sub-Neptunes formed beyond the water-ice line, where the available
planet types of different bulk densities (Luque & Pallé 2022). Despite solid mass increases due to the condensation of volatile species. They
the discovery of thousands of sub-Neptunes, their compositions and accrete large quantities of ice/water content (typically ∼ 50% of their
formation pathways are still highly debated. mass) and then migrate inwards to the orbital locations we observe
Two major hypotheses exist to explain these observations. In the them at today.
“gas dwarf” scenario, super-Earths and sub-Neptunes formed from a Seeking evidence in favour of either model at the population level
single population of progenitor planets with thick primordial H2 /He- is a complex challenge, largely due to the compositional degeneracy
when only measuring a planet’s bulk density (e.g. Rogers 2015; Bean
et al. 2021). Models for both mechanisms can reproduce demographic
★ E-mail: jr2011@cam.ac.uk observations in various parameter spaces, such as orbital period,

© 2025 The Authors

2 J. G. Rogers
planet radius, planet mass, and stellar mass spaces (e.g. Gupta & of 𝑇eq . The sound speed in Equation 1 is evaluated at the radiative-
Schlichting 2019; Rogers & Owen 2021; Burn et al. 2024a). Water convective boundary, 𝑐2s ≡ 𝑘 B𝑇eq /𝜇, where 𝜇 is the mean molecular
worlds are predicted to be more prominent around M-dwarfs since weight of the envelope, and is a crucial variable in this study. We
the water-ice line may be closer to the host star, providing a larger highlight here that, for simplicity, we assume the mean molecular
reservoir of water-rich material closer to the location where sub- weight to be constant in the envelope as a function of radius. We will
Neptunes are observed. However, observational evidence for this frequently refer to 𝜇 as the “bulk” mean molecular weight, which can
remains unclear (Rogers et al. 2023b; Ho et al. 2024). be interpreted as a mass-averaged value throughout the envelope. The
In this paper, we present a new approach to distinguishing between density profile in this radiative layer is given by:
gas dwarfs and water worlds, which relies on the characterisation of   
𝐺 𝑀c 1 1
planet demographics at early stellar ages. We exploit the difference in 𝜌 = 𝜌rcb exp − . (2)
predicted mean molecular weights in the envelopes of such planets, 𝑐2s 𝑟 𝑅rcb
which produce very different evolutionary behaviour at young ages. To evolve the model, we numerically solve a set of coupled ordinary
Gas dwarf envelopes are expected to be H2 /He-dominated, yielding differential equations:
a relatively low mean-molecular weight and, thus, large scale height.
𝑑𝑀env ¤
The thickness of this envelope is sensitive to the thermal state of = − 𝑀, (3)
the planet, which will be larger when the planet is young. Thus, 𝑑𝑡
gas dwarfs are predicted to be puffy at young ages. Water worlds,
on the other hand, are expected to have volatile-rich envelopes with 𝑑𝐸 p
= −𝐿 rad . (4)
high mean molecular weights and small scale heights. The thick- 𝑑𝑡
ness of such envelopes will respond weakly to changes in thermal Here, Equation 3 represents the conservation of mass, where 𝑀env is
state, meaning that they will not undergo significant contraction with the envelope mass and 𝑀¤ is the mass loss rate. Equation 4 represents
time (see also Lopez 2017; Aguichine et al. 2024). This provides the conservation of planet energy, 𝐸 p .1 Here, 𝐿 rad is the radiative
an observational difference between the models, which can be ex- luminosity of the planet, which, in radiative-convective equilibrium,
ploited. Recent atmospheric characterisation of evolved, temperate is given by:
sub-Neptunes from JWST have revealed a wide range in upper at-
4
64𝜋∇ad 𝐺 𝑀p 𝜎𝑇eq
mospheric mean molecular weights (e.g. Madhusudhan et al. 2023;
𝐿 rad = , (5)
Benneke et al. 2024; Piaulet-Ghorayeb et al. 2024; Davenport et al. 3 𝑐2s 𝜅 rcb 𝜌rcb
2025). Meanwhile, JWST/HST observations of young sub-Neptunes
where 𝜅rcb is the gas Rosseland mean opacity evaluated at the
have suggested low mean molecular weights (e.g. Barat et al. 2024;
radiative-convective boundary.
Thao et al. 2024). However, these observations are sparse, biased
Calculating the envelope mass 𝑀env and planetary energy 𝐸 p re-
towards M-dwarfs, and only probe the upper atmospheres of sub-
quires integrating through the planetary structure. For envelope mass:
Neptunes, preventing robust constraints from being placed on forma-
tion pathways at the population level. ∫ 𝑅out
Our paper is laid out as follows: in Section 2, we present a semi-
𝑀env = 4𝜋𝑟 2 𝜌(𝑟)𝑑𝑟, (6)
analytical model for gas dwarf and water world structure and evolu- 𝑅c
tion, and justify that their radii differ greatly at young ages in Section
where 𝑅c is the core radius and 𝜌 is given by Equations 1 and 2. The
3. In Section 4, we compare population-level models with observed
outer boundary, 𝑅out is given by:
planets from Kepler and TESS. Discussion and conclusions are pro-
vided in Sections 5 and 6, respectively. 𝑅out = min(𝑅B , 𝑅H ), (7)
where 𝑅B = 𝐺 𝑀c /2𝑐2s is the Bondi radius and 𝑅H = 𝑎(𝑀c /3𝑀∗ ) 1/3
is the Hill radius. Planetary energy is calculated by:
2 A SIMPLE MODEL FOR GAS DWARFS AND WATER ∫ 𝑅out 
𝑘 B𝑇 (𝑟) 𝐺 𝑀c

WORLDS 𝐸p = 4𝜋𝑟 2 𝜌(𝑟) − 𝑑𝑟
𝑅c 𝜇 (𝛾 − 1) 𝑟
(8)
The envelopes of close-in, highly irradiated sub-Neptunes are ex- 𝑀c 𝑘 B𝑇c
+ ,
pected to be in radiative-convective equilibrium. To investigate their 𝜇c (𝛾c − 1)
evolution, we construct semi-analytic models based on previous
where the first and second terms represent energy stored in the en-
works for H2 /He-dominated envelopes (e.g. Owen & Wu 2017;
velope and core, respectively. For the core, we follow Misener &
Gupta & Schlichting 2019; Misener & Schlichting 2021) and steam-
Schlichting (2021) and assume an incompressible and isothermal
dominated envelopes (Rogers et al. 2025). Assuming an ideal gas
core at a temperature of 𝑇c with mean molecular weight 𝜇 𝑐 and
equation of state, the density profile in the convective envelope is
adiabatic index 𝛾c = 4/3 from the Dulong-Petit law.
assumed to be adiabatic and given by (e.g. Ginzburg et al. 2016):
To solve Equations 3 and 4, we use a RK45 adaptive-step integrator
   1
𝛾−1
with a numerical tolerance of 10 −6 . At each timestep, the updated
𝐺 𝑀c 𝑅rcb
𝜌(𝑟) = 𝜌rcb 1 + ∇ad 2 −1 , (1) envelope mass, 𝑀env , and planet energy, 𝐸 p , are used to solve for the
𝑐 s 𝑅rcb 𝑟
new radiative-convective boundary density, 𝜌rcb , and radius, 𝑅rcb .
where the subscript ‘rcb’ refers to quantities evaluated at the
radiative-convective boundary. Here, 𝜌 is the gas density, 𝛾 is the 1 Strictly speaking, Equation 4 should have an extra term proportional to
adiabatic index, ∇ad ≡ (𝛾 − 1)/𝛾 is the adiabatic temperature gra- ¤ since escaping material also removes energy from the
the mass loss rate, 𝑀,
dient, 𝑀c is the planet’s core mass and 𝑅rcb is the position of the planet. This is important for extreme mass loss processes such as boil-off (e.g.
radiative-convective boundary. Above the convective region sits a ra- Owen & Wu 2016). However, our evolution models begin after this phase.
diative layer, assumed to be isothermal at an equilibrium temperature Hence we can safely ignore this term.

MNRAS 000, 1–12 (2025)

 On the road to the radius valley 3
We also calculate the position of the photospheric radius, 𝑅ph , and weight of 𝜇c = 60, motivated by Earth (Szurgot 2015). To calculate
transit radius, 𝑅tr , via RK45 numerical integration (e.g. Guillot 2010). the core radius, 𝑅c , we adopt the mass-radius relations of Fortney
This is required since the difference between photospheric and transit et al. (2007), assuming an Earth-like interior composition of 25%
radii of young planets with large scale heights can become significant iron, 75% rock, as inferred from Rogers & Owen (2021); Rogers
(≳ 10%). et al. (2023a).
Since we are interested in young planets and their host stars, we Hydrogen-dominated gas is susceptible to atmospheric escape,
include pre-main sequence stellar evolution from Johnstone et al. which can include mechanisms such as boil-off during disc dispersal
(2021). The stellar luminosity, 𝐿 ∗ (𝑡), for a given stellar mass and age, (Owen & Wu 2016; Ginzburg et al. 2016; Rogers et al. 2024a), core-
is used to calculate planetary equilibrium temperature as follows: powered mass loss (e.g. Ginzburg et al. 2018; Gupta & Schlichting
 1 2019; Misener et al. 2024) and XUV photoevaporation (e.g. Owen &
𝐿∗ 4
𝑇eq = , (9) Jackson 2012; Kubyshkina & Fossati 2022; Schulik & Booth 2022;
16𝜎𝜋𝑎 2 Caldiroli et al. 2022). For simplicity, we choose to model planets af-
where we have assumed a zero Bond albedo. ter the protoplanetary disc has fully dissipated. As shown in Rogers
For initial conditions, we take a simple, agnostic approach in set- et al. (2024a), boil-off and core-powered mass loss should dominate
ting each planet’s cooling timescale, 𝜏cool , equal to the planet’s age, during disc dispersal, then photoevaporation dictates the planetary
𝑡 pl . The cooling timescale, which is closely related to the Kelvin- outflow once the gas disc is optically thin to XUV photons through
Helmholtz timescale used in other works (e.g. Owen & Wu 2017; its midplane. Nevertheless, photoevaporation and core-powered mass
Owen 2020), is given by: loss are similar hydrodynamic escape processes (Owen & Schlicht-
ing 2024), and thus, this simplifying assumption has little effect on
𝐸 c − 𝐸 env
𝜏cool = , (10) our results. We calculate photoevaporative mass loss rates with the
𝐿 rad following (Lecavelier Des Etangs 2007; Erkaev et al. 2007):
where the energy of the envelope 𝐸 env is always negative if the en- 3
𝐿 𝜋𝑅ph
velope is gravitationally bound. This initial condition of 𝜏cool = 𝑡pl 𝑀¤ = 𝜂 XUV , (11)
enforces an initial planet entropy appropriate for its age. In the early 4𝜋𝑎 2 𝐺 𝑀p
stages of planetary evolution, the planet’s thermal state and cooling where 𝜂 is the mass loss efficiency, calculated by interpolating
rate are closely tied to its formation process and contraction. Assum- over the grid of hydrodynamic models from Owen & Jackson
ing the cooling timescale is comparable to the planet’s age ensures (2012). Typical efficiency values for sub-Neptunes can range be-
that the initial model reflects a system for which cooling balances tween ∼ 10 −3 − 10 −1 , depending on planet mass and radius. The
with the internal thermal content over the planet’s lifetime. In other XUV luminosity from the host star, 𝐿 XUV , is taken from the pre-
words, if 𝜏cool ≪ 𝑡pl , the planet will start very inflated and contract main sequence stellar evolution models from Johnstone et al. (2021).
unphysically quickly. On the other hand, 𝜏cool ≫ 𝑡pl would imply a
very contracted initial size and very little thermal evolution with time,
which is likely also unphysical for a recently formed planet. We note 2.2 Water Worlds
here that extreme mass loss events, such as boil-off during protoplan-
For our simplest water world models, we assume a pure steam enve-
etary disc dispersal, may increase a planet’s initial cooling timescale
lope with a mean molecular weight of 𝜇 = 18 amu, a gas adiabatic
since high entropy material is removed from the planet very quickly,
index of 𝛾 = 4/3, and Rosseland mean opacities from Kempton et al.
thus introducing premature cooling (Owen & Wu 2016; Ginzburg
(2023) for pure H2 O gas. For the core, we assume an interior compo-
et al. 2016; Owen 2020; Rogers et al. 2024a; Tang et al. 2024).
sition of 50% rock, 50% water, as predicted by population synthesis
To account for this, we also enforce that a planet’s initial radiative-
models (e.g. Mordasini et al. 2009; Raymond et al. 2018; Zeng et al.
convective boundary cannot exceed 0.1𝑅B , where 𝑅B = 𝐺 𝑀p /2𝑐2s is
2019; Venturini et al. 2016; Mousis et al. 2020; Burn et al. 2024a).
the Bondi radius. This constraint is appropriate for a planet in which
This yields a core mean molecular weight of 𝜇c = 27 amu. We use
a bolometrically-driven transonic wind, triggered by disc dispersal,
this composition in the mass-radius relations of Fortney et al. (2007)
has largely shut off (see Owen & Wu 2016). These initial conditions
to calculate a core radius, 𝑅c . Similar to the steam models of Rogers
are thus appropriate for a planet for which its nascent protoplanetary
et al. (2025), we choose not to consider the atmospheric escape of
disc has dispersed and any consequential mass loss has concluded.
steam-dominated (or, more generally, high mean molecular weight)
In addition to the initial conditions, various model parameters,
envelopes due to the many uncertainties of this mechanism. Little
including the core’s mean molecular weight, 𝜇c , gas mean molecular
is currently known about how processes such as molecular disso-
weight, 𝜇, gas opacity, 𝜅, and gas adiabatic index, 𝛾, depend on the
ciation/ionisation, drag, diffusion, and line-cooling may affect the
assumed formation model. These are discussed in Sections 2.1 and
escape of high mean molecular weight envelopes, and thus, we leave
2.2.
this for future work. Nevertheless, studies have shown that water-rich
sub-Neptunes should store most of their water in their rocky interiors
2.1 Gas Dwarfs (Dorn & Lichtenberg 2021; Luo et al. 2024), meaning that we need
only consider small envelope mass fractions, which may either suffer
For our simplest gas dwarf models, we assume a gas mean molecular from inefficient atmospheric escape and/or may be replenished from
weight of 𝜇 = 2.35 atomic mass units (amu), appropriate for Solar- their outgassing interiors.
metallicity gas (e.g. Anders & Grevesse 1989), a gas adiabatic index
of 𝛾 = 7/5 for molecular hydrogen, and interpolate the Rosseland
mean opacity tables from Kempton et al. (2023) for Solar-metallicity
3 THE IMPACT OF ENVELOPE MEAN MOLECULAR
gas. For the cores of gas dwarfs, multiple studies have shown that
WEIGHT ON PLANET SIZE
their interiors are consistent with an Earth-like composition (e.g.
Owen & Wu 2017; Gupta & Schlichting 2020; Rogers & Owen The mean molecular weight of a gaseous envelope directly controls
2021; Rogers et al. 2023a). Thus, we assume a core mean molecular its scale height, 𝐻. In the simplistic case of an isothermal ideal gas,

MNRAS 000, 1–12 (2025)

4 J. G. Rogers
atmospheric bulk mean molecular weights, which are dictated by
planet formation and evolution.

3.1 Gas dwarf vs. water world evolution

We now consider the evolution of gas dwarfs and water worlds. Figure
2 demonstrates their evolution in orange and blue, respectively. We
show planet models with core masses of 4𝑀⊕ , 6𝑀⊕ and 8𝑀⊕ in
solid, dashed and dotted lines, respectively. Models are initialised at
𝑡 = 5 Myrs. In the left-hand panel, we construct “minimal” models,
as described in Sections 2.1 and 2.2, which represent the simplest
physical construction of each formation scenario. Under this set of
models, gas dwarfs and water worlds have bulk mean molecular
weights of 𝜇 = 2.35 amu and 𝜇 = 18.0 amu, respectively, initial
envelope mass fractions of 𝑋 = 0.05 and 𝑋 = 0.01, respectively, and
orbital periods of 10 days. Gas dwarfs and water worlds converge
in size after ∼ 1 Gyr due to their internal composition producing
the same radius: water worlds with water-rich, low-density cores and
compact envelopes, gas dwarfs with rocky/iron-rich, high-density
cores and extended envelopes. However, their evolution is divergent
at early ages. The low bulk mean molecular weight of gas dwarf
envelopes means that the energy of formation, which is abundant at
early times, puffs the planet up significantly when compared to water
Figure 1. The transit radius evolution of sub-Neptune envelopes for different
bulk mean molecular weights of 𝜇 = 2.35, 5.0, 10 and 18.0 amu in blue,
worlds. This is particularly apparent for planet ages ≲ 100 Myrs.
green, orange and red, respectively. All planets have a 5𝑀⊕ core, with an In the right-hand panel of Figure 2, we consider more complex
Earth-like core composition of 25% iron, 75% silicate. Each model has an models of gas dwarfs and water worlds, informed by recent theoret-
envelope mass fraction of 𝑋 ≡ 𝑀env /𝑀c = 0.01 and a constant equilibrium ical and observational studies. In the case of gas dwarfs, multiple
temperature of 𝑇eq = 800 K and does not undergo any mass loss. works have shown that the chemical interaction of Solar-metallicity
gas with a magma ocean at high temperatures will result in an in-
crease in envelope mean molecular weight (Kite & Schaefer 2021;
the scale height is: Schlichting & Young 2022; Misener & Schlichting 2022; Misener
𝑐2 𝑘 𝑇 𝑟2 1 et al. 2023; Rogers et al. 2024b; Young et al. 2024). Most notably,
𝐻 (𝑟) = s = B ∝ , (12) chemical reactions between hydrogen gas and molten silicate can pro-
𝑔 𝜇 𝐺 𝑀p 𝜇
duce large quantities of vapour species such as SiO, MgO, FeO and
where 𝑔 is the gravitational acceleration at radius 𝑟. The scale height H2 O, increasing the mean molecular weight above Solar values. As
is inversely proportional to the mean molecular weight for a given an example, recent JWST observations of TOI-270 d, an archetypal
temperature (or, equivalently, thermal state of the envelope). The temperate sub-Neptune, have been interpreted with a mean molecu-
mean molecular weight also impacts the scale height for a non- lar weight of 𝜇 ∼ 5 amu (Benneke et al. 2024), consistent with the
isothermal envelope, such as in the convective interior. From Equa- global chemical equilibrium models of Schlichting & Young (2022);
tion 1, one can show that the pressure gradient, which is analogous Rogers et al. (2024b). To capture this effect in our simple evolu-
to the scale height, |𝑑𝑃/𝑑𝑟 | ∝ 𝜇1/(𝛾−1) . In other words, the enve- tionary formalism, we produce gas dwarf models with bulk mean
lope pressure in a convective interior falls off more steeply for a high molecular weights of 𝜇 = 5.0 amu, which is indicative of mixed
mean molecular weight composition, implying a less puffy envelope miscible envelopes. Since H2 O is a dominant opacity carrier, we
for a given mass and entropy. In summary, two envelopes of the same use the Rosseland mean opacities from Kempton et al. (2023) for a
mass and entropy but with different mean molecular weights will, H2 /He/H2 O mixture with a H2 O mass fraction of 60% to reproduce
therefore, be different sizes. this bulk mean molecular weight. We continue to assume atmospheric
As a sketch of the idea, Figure 1 demonstrates the evolution of a escape according to Equation 11 with the same efficiencies as in the
sub-Neptune for different bulk mean molecular weights, ranging from minimal gas dwarf model, but highlight that molecular cooling and
𝜇 = 2.35 amu (appropriate for Solar metallicity H2 /He-dominated drag should reduce such mass loss rates. This model thus provides
gas) to 𝜇 = 18.0 amu (appropriate for a pure H2 O gas). All planets an upper limit in terms of contraction rate as a test of whether the
have a 5𝑀⊕ core, with an interior Earth-like core composition. Each mixed-envelope models can still be distinguished at early ages. In or-
model has an envelope mass fraction of 𝑋 ≡ 𝑀env /𝑀c = 0.01. der for this model to reproduce the typical sizes of sub-Neptunes, we
For simplicity, we assume a constant equilibrium temperature of must increase its envelope mass fraction, as highlighted in Thorngren
𝑇eq = 800 K (ignoring the pre-main sequence evolution of the host & Fortney (2019); Misener & Schlichting (2022). For the right-hand
star for now), a constant gas opacity of 𝜅 = 0.1 cm2 g −1 and do not panel of Figure 2, gas dwarfs have an initial envelope mass fraction
include any mass loss. Whereas the size of planets converge after of 𝑋 = 0.3. Note that it is beyond the scope of this work to model the
Gyrs of evolution, they differ significantly at early times. The high convection-inhibition brought about from mean molecular weight
temperatures from formation increase the size of planets with lower gradients in the lower envelope, as also highlighted in Misener &
mean molecular weight more substantially, causing them to be more Schlichting (2022); Misener et al. (2023). Our goal is to demonstrate
inflated when young. In this work, we highlight that a population the population-level impacts of bulk mean molecular weight on the
of sufficiently young transiting exoplanets would allow one to probe sizes of young sub-Neptunes. Nevertheless, convection inhibition can

MNRAS 000, 1–12 (2025)

 On the road to the radius valley 5

Figure 2. The radius evolution for different sub-Neptune models: gas dwarfs with Earth-like interiors and H2 /He-dominated envelopes (orange) and water
worlds with water-rich interiors and steam-dominated atmospheres (blue). Planet size is shown for core masses of 4𝑀⊕ , 6𝑀⊕ and 8𝑀⊕ in solid, dashed and
dotted, respectively. The left-hand panel shows “minimal” models, representing the simplest construction of each planet formation hypothesis. Gas dwarfs have
water-poor, Earth-like interiors, host Solar-metallicity envelopes with bulk mean molecular weights of 𝜇 = 2.35 amu, initial 5% envelope mass fractions and
undergo photoevaporative mass loss. Water worlds have water-rich interiors, host pure steam envelopes with bulk mean molecular weights of 𝜇 = 18 amu
with envelope mass fractions of 1% and do not undergo mass loss. The right-hand panel shows mixed envelope models, for which gas dwarfs have bulk mean
molecular weights of 𝜇 = 5 amu and initial envelope mass fractions of 30%. Water worlds have bulk mean molecular weights of 𝜇 = 6 amu with envelope mass
fractions of 1%. Despite converging to the same radius after ∼ Gyrs of evolution under various assumptions, gas dwarfs and water worlds behave very differently
at early times.

alter the thermal evolution of planets, and should be considered in and assumed an adiabatic equation of state for gas pressure, 𝑃 ∝ 𝜌 𝛾 .
future work. The internal energy for an ideal gas envelope scales as 𝐸 i ∝ 𝑃¯env𝑉env ,
Similar to our mixed-envelope models of gas dwarfs, the right- where 𝑉env is the volume of the atmosphere and 𝑃¯env is a mass-
hand panel of Figure 2 also shows water worlds with H2 /He/H2 O averaged envelope pressure (e.g. Kippenhahn et al. 2012). In the
mixed envelopes. Recent population synthesis models from Burn case of gas dwarfs with non-self-gravitating envelopes, particularly
et al. (2024a,b) have considered the accretion of solar-metallicity at young ages, the volume is dominated by the radius of the envelope,
gas onto water-rich sub-Neptunes. They predict envelope water mass hence 𝑉env ∝ 𝑅p3 and the pressure can be approximated as 𝑃¯env ∝
fractions as low as ∼ 70%, equating to a mean molecular weight (𝑀env /𝑅p3 ) 𝛾 from the adiabatic equation of state. Combining, and
of 𝜇 ≈ 6 amu, which we adopt in our mixed-envelope models of recalling that 𝑀env ≡ 𝑋 𝑀c , one finds that:
water worlds. We use the Rosseland mean opacities from Kempton
2 −1
et al. (2023) for a H2 /He/H2 O mixture with a H2 O mass fraction 𝑅p ∝ 𝑋 3 𝑀c 3 . (14)
of 70%. We retain an atmospheric mass fraction of 𝑋 = 0.01 from
our “minimal” models, motivated by theoretical studies of water- In other words, the radius of a planet is insensitive to its mass for a
rich planets hosting most of their H2 O in the interior (Luo et al. given atmospheric mass fraction and age.2 This is well-documented
2024). Recent JWST upper-atmospheric characterisation of GJ 9827 in other works (see Figure 1 from Lopez & Fortney 2014).
d reveals an abundance of H2 O and a mean molecular weight between
𝜇 ∼ 10 and 18 amu (Piaulet-Ghorayeb et al. 2024). This planet is
thus consistent with some predicted water world scenarios. 4 YOUNG TRANSITING EXOPLANETS AS A PROBE OF
Regardless of the underlying assumptions, Figure 2 highlights that PLANET FORMATION AT THE POPULATION LEVEL
the radius evolution of gas dwarfs and water worlds are dramatically
different at early ages. One can also see that planet size is relatively In Section 3 we have demonstrated that the radius evolution of sub-
insensitive to planet mass as a function of time. This can be un- Neptune envelopes can vary substantially with different bulk mean
derstood from a simple virial theorem argument, which states that, molecular weights. We now build population-level models under the
for a spherically-symmetric envelope in hydrostatic equilibrium, its two categories, gas dwarfs vs. water worlds, to showcase how such
gravitational potential energy, 𝑈 = −𝐺 𝑀c 𝑀env /𝑅p , is related to its
internal energy, 𝐸 i : 2 Equation 14 is derived assuming an adiabatic index of 𝛾 = 5/3, since
𝑈 = −3(𝛾 − 1)𝐸 i , (13) 𝛾 = 4/3 provides an analytic singularity. When the full semi-analytic structure
equations are numerically solved, as in Figure 2, one recovers an insensitive
where we have safely ignored vanishing boundary terms at large radii planet size as a function of planet mass for a fixed envelope mass fraction.

MNRAS 000, 1–12 (2025)

6 J. G. Rogers
scenarios would present at early times. For simplicity, we only con- in the survey. For each sub-Neptune formation model, we repeat this
struct population models from the “minimal” formation scenarios process to produce 1000 synthetic surveys. We then use a kernel
since these are better studied when compared to mixed-envelope density estimator (KDE) to calculate a probability density function
models as described in Sections 2.1 and 2.2. As such, we assume (PDF) for planet detection in period-radius space from the combined
that gas dwarfs and water worlds have bulk envelope mean molec- surveys. This PDF is then scaled to the integrated expectation value
ular weights of 𝜇 = 2.35 amu and 18.0 amu, respectively, for the for survey yield, i.e. the number of detected planets in each modelled
remainder of this paper. Nevertheless, Figure 2 demonstrates that survey. An underlying formation model is a good fit to the data if
more complex mixed-envelope models also show similar differences this expectation value approximately matches the actual yield of 926
in radius evolution at early times. planets in the CKS catalogue, as well as the distribution shape of
observed planets in period-radius space.
Given a population-level model that reproduces the late-time de-
4.1 Modelling synthetic transit surveys mographics, the next step is to inspect the model predictions at early
The prerequisite for our population models is that they reproduce the times. We choose to compare these models with a small sample of
observed exoplanet population at late times ≳ 1 Gyr. To do so, we observed transiting exoplanets from TESS (Vach et al. 2024). Figure
choose to compare the models at late times with the California Ke- 2 highlights that the gas dwarf and water world models diverge the
pler Survey (CKS) from Petigura et al. (2022). We seek to construct most at the earliest times, hence we only include planets within this
underlying models that, when injected through the CKS pipeline, re- sample with host stellar ages < 40 Myrs.3 We evaluate each mod-
produce the observed (biased) population. For each model, we evolve elled planet radius at a random time, uniformly drawn up to 40 Myrs
106 planets, with stellar masses drawn from a Gaussian distribution and again add a radius measurement uncertainty of 5%. In the TESS
with a mean of 1𝑀⊙ and a standard deviation of 0.3𝑀⊙ (Rogers survey of Vach et al. (2024), the stellar sample for stars < 40 Myrs
et al. 2021). We also randomly draw a disc dispersal time, which contains ≈ 2500 stars. In order to produce synthetic surveys, we again
defines the initial time for each planet model, uniformly drawn be- assume the same occurrence rate of 0.7 planets per star, and draw
tween 3 Myrs and 10 Myrs (e.g. Kenyon & Hartmann 1995; Ercolano 1750 planets from our underlying population of young planets. Since
et al. 2011; Koepferl et al. 2013). The typical ages of host stars in Kepler and TESS have different completeness and noise properties,
the CKS catalogue are in the range of ≈ 1 − 10 Gyrs (e.g. Berger we approximate the probability of detection for TESS with the same
et al. 2020). Hence, we terminate each planet model at a random functional form as in Equation 16 with 𝐴 = 0.85, 𝑘 = 5.0, 𝑙 = 1.00
time, uniformly drawn between these limits. Finally, we statistically and 𝜃 = 2.5, which best reproduces Figure 8 from Vach et al. (2024)
model observational uncertainty in planetary radii by introducing a for young FGK stars. Similarly, the signal-to-noise ratio for TESS
random Gaussian perturbation to each planet radius with zero mean transits is calculated with Equation 17 with 𝑇obs = 60 days. Al-
and a standard deviation of 5% the planet’s true size. though TESS does not quantify photometric precision with a CDPP,
We produce synthetic CKS surveys at late times by randomly draw- we find an effective log(CDPP) randomly drawn from a Gaussian
ing 25, 200 planets from the underlying population of 106 modelled distribution with a mean and standard deviation of −1.5 and 0.25, re-
planets. This value comes from an assumed occurrence rate of 0.7 spectively, best reproduces completeness maps from Fernandes et al.
planets per star in the CKS stellar catalogue of 36, 000 stars (Fulton (2022); Vach et al. (2024). As with the CKS synthetic surveys, we
et al. 2017). Thus, the expectation value for planets orbiting the sur- repeat this entire procedure 1000 times and produce a combined
veyed stars is 0.7 × 36, 000 = 25, 200. For each planet, we calculate PDF for planet detection, scaled to the integrated expectation value
the probability of transit and detection following Fulton et al. (2017); of synthetic survey yield. This is compared to the 14 planets with
Rogers & Owen (2021); Petigura et al. (2022). The probability of hosts stars < 40 Myrs in the Vach et al. (2024) survey.
transit is given by: We stress that comparing models at early evolutionary times is
only valid if the population models reproduce the observations at
𝑅∗
𝑝 tr = , (15) late times. In order to reproduce the CKS sample at ∼ 1 − 10 Gyrs,
𝑎 both population-level models require different planetary initial con-
where 𝑅∗ is the stellar radius and 𝑎 is the semi-major axis. The proba- ditions, including initial envelope mass fraction distributions, core
bility of detection, 𝑝 det , takes the form of a Γ cumulative distribution mass distributions, and orbital period distributions. We discuss these
function of the form: in the following sections.
∫ 𝑚−𝑙
𝜃
𝑝 det = 𝐴 Γ(𝑘) 𝑡 𝑘−1 𝑒 −𝑡 d𝑡, (16)
0
4.2 Gas dwarf population model
where 𝐴 = 1.0, 𝑘 = 17.56, 𝑙 = 1.00 and 𝜃 = −0.49 following Fulton
et al. (2017). The signal-to-noise ratio, 𝑚, for a transiting planet is: For the gas dwarf case, both super-Earths and sub-Neptunes are
  2 √︂ drawn from the same underlying population, and atmospheric es-
𝑅p 𝑇obs 1 cape is responsible for separating the two sub-populations over time
𝑚= , (17)
𝑅∗ 𝑃 CDPP and forming the radius valley. We use the inferred underlying dis-
where 𝑇obs = 4 years is the observation time for Kepler, 𝑃 is the tributions from Rogers & Owen (2021), who fitted the XUV pho-
orbital period, and the combined differential photometric precision toevaporation “minimal” model from Section 2.1 to the CKS data
(CDPP) is a measure of noise for the host star, measured in ppm. of Fulton & Petigura (2018) using a Bayesian hierarchical model.
To match the results from Koch et al. (2010), we randomly draw Both the core mass and initial envelope mass fraction distributions
log(CDPP) as a Gaussian distribution with a mean and standard are parameterised using Bernstein polynomials of 5th order. These
deviation of −4.0 and 0.25, respectively. The product of 𝑝 tr × 𝑝 det
is the probability that a given planet will be observed in a transit
survey. For each planet, we draw a random number between 0 and 1. 3 We also update this sample to include the recently detected 3 Myr old
If this number is less than 𝑝 tr × 𝑝 det , then this planet is “observed” planet from Barber et al. (2024).

MNRAS 000, 1–12 (2025)

 On the road to the radius valley 7
are peaked distributions, with centres at ∼ 4𝑀⊕ and ∼ 2%, respec- One can see the emergence of the radius valley in the gas dwarf
tively. The orbital period distribution is parameterised with a smooth case, with the ratio of super-Earths to sub-Neptunes increasing with
broken power-law as follows: time as atmospheric escape strips primordial sub-Neptunes into bare
𝑑𝑁 1 rocky super-Earths. The typical sizes of gas dwarf sub-Neptunes
∝ 𝑃 , (18) also significantly decrease with time as a result of their low bulk
𝑑𝑃 ( 𝑃0 ) −𝑘 1 + ( 𝑃𝑃0 ) −𝑘2 mean molecular weights. The water world scenario, on the other
where 𝑃0 = 5.5 days, 𝑘 1 = 2.3, and 𝑘 2 = 0.0. We assume all cores hand, shows very little evolution since the super-Earths and sub-
have a composition of 25% iron, 75% rock (Rogers & Owen 2021; Neptunes formed separately in the disc phase, and the water world
Rogers et al. 2023a) with Solar-metallicity H2 /He envelopes with sub-Neptunes have high bulk mean molecular weight, meaning very
mean-molecular weights of 𝜇 = 2.35 amu. Again, we stress that little thermal contraction.
models that include self-consistent interior-envelope chemical cou- Figure 4, on the other hand, shows the predicted biased survey
pling are still in their infancy (e.g. Kite & Schaefer 2021; Schlichting yields at late and early times, as described in Section 4.1. These are
& Young 2022; Misener et al. 2023). They are also computationally compared to the CKS exoplanet sample from Petigura et al. (2022),
expensive, meaning that population modelling is currently challeng- and the young TESS sample from Vach et al. (2024). In this case,
ing. While we leave these considerations for future work, Figure 2 colour bars represent the expectation value for integrated survey
highlights that the difference in early evolution between gas dwarfs yield. Both models reproduce the observations at late times, with two
and water worlds under this scenario is still significant, meaning our distinct sub-populations of super-Earths and sub-Neptunes, separated
results will remain largely the same. by the radius valley. The expectation value for survey yield is 910+20
−16
for gas dwarfs, and 931+24
−14
for water worlds, which are both consistent
with the real CKS survey yield of 926.
4.3 Water world population model
The nature of TESS’ lower photometric precision and the higher
The fundamental prediction of the water world hypothesis is that activity of young stars means that the expected yields of such surveys
super-Earths and sub-Neptunes formed through different channels. are considerably lower than Kepler. The survey completeness is sig-
Therefore, we require different underlying distributions for the two nificantly lower for planets smaller than ∼ 3𝑅 ⊕ , meaning that super-
classes of planets. For simplicity, we choose not to model these plan- Earths are very unlikely to be detected under either model. As an
ets during the disc phase, as done in population-synthesis modelling example, a planet of radius 3𝑅 ⊕ at 0.1 AU orbiting a main-sequence
(e.g. Burn et al. 2024a), which includes large-scale disc migration 1𝑀⊙ star is predicted to have a signal-to-noise ratio of 𝑚 ∼ 100 from
to bring water worlds interior to the water-ice line. Instead, we pro- Equation 17 as detected by Kepler, but 𝑚 ∼ 5 if orbiting a young
duce appropriate underlying distributions for the post-disc evolution active star detected by TESS. This equates to a probability of detec-
phase that reproduce the CKS data at ∼ Gyr times. To match the ratio tion from Equation 16 of 𝑝 det ∼ 1 and 𝑝 det ∼ 0.07, respectively. As
of super-Earths and sub-Neptunes in these observations, 35% of the a result of this systematic reduction in detectability for young active
underlying planet population are modelled as super-Earths, 65% as stars observed by TESS, the expectation value for survey yield under
sub-Neptunes. the gas dwarf model is significantly reduced to 15+3 planets. The
−2
For the core mass distribution, we adopt a log-Gaussian distribu- actual yield from TESS is 14 planets < 40 Myrs in Vach et al. (2024).
tion: For the water world model, the the expected yield is 3+3 planets,
−1
𝑑𝑁

(log 𝑀c − log 𝜇M ) 2
 which is inconsistent with TESS at the ∼ 2 − 3𝜎 level. The gas dwarf
∝ exp − 2
. (19) model is better at reproducing the observations since these planets
𝑑 log 𝑀c 2𝜎M
are larger at younger ages due to their lower bulk mean molecular
For super-Earths, we set the mean and standard deviation to 𝜇M = weights, meaning they are more likely to be detected. One can see
3𝑀⊕ and 𝜎M = 0.35, respectively. For sub-Neptunes, we set 𝜇M = that the gas dwarf model predicts planet detections at significantly
10𝑀⊕ and 𝜎M = 0.30, respectively. larger radii at early times, as seen in the TESS data.
For the orbital period distributions, we use Equation 18 with 𝑃0 = Although the young transiting exoplanet sample from TESS is
4.5 days, 𝑘 1 = 2.0, and 𝑘 2 = −0.5 for super-Earths, and 𝑃0 = visually more consistent with the gas dwarf model, it is interesting
8.5 days, 𝑘 1 = 2.5, and 𝑘 2 = 0.1 for sub-Neptunes. These values are to note that the predicted orbital period distribution is noticeably
all consistent with the CKS study of Petigura et al. (2022). different to the data. Young TESS planets pile up at ≳ 10 days, despite
As discussed, the water world hypothesis predicts that super-Earths completeness being significantly higher at shorter orbital periods.
and sub-Neptunes have different core compositions due to their dif- This specific realisation of the gas dwarf model thus over-predicts the
ferent formation channels. We assume Earth-like compositions of number of sub-Neptunes at short orbital periods at young ages. The
25% iron, 75% rock for super-Earths with no atmospheres (Rogers & majority of these predicted close-in sub-Neptunes would eventually
Owen 2021; Rogers et al. 2025). For sub-Neptunes, we adopt a core lose their atmospheres to become super-Earths on ∼ 100 Myr to
composition of 50% water, 50% rock and draw a random steam en- ∼ Gyr timescales. The discrepancy could thus be caused by an under-
velope mass fraction log-uniformly between 0.001% and 10% (Luo prediction of planets that were “born rocky”, i.e. began life as super-
et al. 2024). Earths in the post-disc epoch. Alternatively, it could be due to post-
disc inward orbital migration that changes the period distribution
with time. We discuss this in Section 5.2. Although not as statistically
4.4 Population-level results
significant, it is also worth noting that the gas dwarf model slightly
The underlying unbiased planet populations for gas dwarf and water under-predicts the number of large sub-Neptunes ∼ 10R ⊕ at young
worlds are presented in the left and right-hand panels of Figure 3, ages. This could be explained by higher initial entropies of formation
respectively. These are presented at late times (1 − 10 Gyrs) and or perhaps delayed radial contraction due to convection inhibition
early times (< 40 Myrs) in the top and bottom panels, respectively. in the envelope of such planets (e.g. Misener & Schlichting 2022),
Colour bars represent the relative probability of planet occurrence. which we discuss in Section 5.4.

MNRAS 000, 1–12 (2025)

8 J. G. Rogers

Figure 3. Population-level predictions of gas dwarf and water world hypotheses for sub-Neptune evolution in the left and right-hand panels, respectively. The
populations are shown at times of 1 − 10 Gyrs and < 40 Myrs in the top and bottom panels, respectively. Colour bars represent relative probability of planet
occurrence.

A full statistical comparison between models and data is beyond and an observational and theoretical route to determine the nature of
the scope of this work, and requires Bayesian model comparison sub-Neptunes at the population level.
that takes into account the many free parameters of each model. In
other words, there are likely realisations of both models that repro-
duce the data marginally better than those presented in this work, 5.1 Gas opacities and cooling timescales
hence a comparison is not fair unless these realisations are fully ex-
An important envelope property that controls the evolution of a planet
plored. Model comparison involves calculating a Bayes factor via the
is gas opacity. Whereas the mean molecular weight influences the in-
posterior distributions of both models, which is a computationally
stantaneous size of an envelope (via the scale height), the gas opacity
expensive exercise and hence left for future work when the young
also affects the instantaneous rate of change of envelope size. From
TESS sample has increased in size. We highlight that the sample of
Equation 5, one can see that the radiative luminosity, which controls
young transiting exoplanets is small, and thus current inferences are
the cooling and contraction of a planet, is inversely proportional to
likely dominated by small number statistics. Nevertheless, various
the opacity. Thus, higher gas opacities yield slower thermal con-
studies have shown a tentative overabundance of large sub-Neptunes
traction. As a demonstration of this, Figure 5 shows the evolution
(≳ 4𝑅 ⊕ ) orbiting young stars, which may hint at low mean-molecular
of the same planet from Figure 1, however, here, we vary the gas
weight, gas dwarf evolution (e.g. Fernandes et al. 2022; Christiansen
Rosseland mean opacity from 0.01 cm2 g −1 to 10.0 cm2 g −1 with a
et al. 2023; Vach et al. 2024; Fernandes et al. 2025).
mean molecular weight of 𝜇 = 5.0 amu. To highlight the difference
in evolution, we set the initial radiative-convective boundary radii to
𝑅rcb = 5𝑅 ⊕ . Again, planets converge to the same size after several
Gyrs, but their initial cooling rate is altered by varying opacities. One
can also see that the planets with the lowest opacities contract with a
5 DISCUSSION
timescale much shorter than their age (seen by near-vertical lines in
We have shown that the bulk mean molecular weight of a small ex- Figure 5). This suggests that the initial condition of 𝑅rcb = 5𝑅 ⊕ was
oplanet envelope directly affects its thermal contraction with time. likely inappropriate for these envelopes. The probability of detect-
At early times this can result in dramatically different evolution de- ing such a planet at its initially large size is thus vanishingly small,
pending on the assumed formation model. In particular, we have given its rapid contraction. Since we have observed several young
focussed on the gas dwarf and water world hypotheses regarding the sub-Neptunes at young ages with large radii ∼ 10 𝑅 ⊕ , it follows from
formation of sub-Neptunes. We now discuss caveats to this strategy, a statistical argument that their cooling timescales must be at least

MNRAS 000, 1–12 (2025)

 On the road to the radius valley 9

Figure 4. Same as Figure 3 but now accounting for survey biases and uncertainty. Colour bars represent the integrated expectation value for survey yield. We
perform synthetic California Kepler Surveys on each model in the top panels, with the real observed population shown as grey points (Petigura et al. 2022).
Similarly, we perform synthetic TESS surveys for host stars < 40 Myrs in the bottom panels, with real observed planets shown as grey points (Vach et al. 2024).

equal to their ages. In our work, we chose to adopt agnostic initial at the youngest possible age. This minimises the number of cooling
conditions in which the cooling timescale of the planet is equal to its timescales a planet has evolved over, and thus captures the state of
age, which in turn sets the initial size of the planet at the point of disc each planet closer to its initial conditions. The major reason that we
dispersal. More work is needed to investigate the initial conditions are unable to infer the underlying nature of sub-Neptunes at ∼ Gyr
of planets under different scenarios, as discussed in Section 5.4. timescales is that they have evolved for many cooling timescales,
To avoid the potential degeneracy between gas opacity and mean destroying the information of their individual formation pathway.
molecular weight when comparing models with observations of Recent detections of very young transiting exoplanets < 10 Myrs
young transiting planets, one should use self-consistent models that (e.g. Barber et al. 2024) have shown that these detections are indeed
couple the two properties. For example, a steam-dominated enve- possible and offer exciting glimpses into newly-formed planets.
lope will have larger mean opacities than solar metallicity H2 /He-
dominated gas. This coupling is also true for the adiabatic index,
5.2 On the benefits of probing the young super-Earth
𝛾, which in this work parametrises the equation of state within a
population
convective envelope (see Equation 1) and can also vary depending
on the species present in an atmosphere. The adiabatic index is im- In Figure 4, we see that the gas dwarf formation model is visually
portant since it controls the size of a planet’s convective envelope better at reproducing the young planet population than the water
for a given entropy. For example, a lower adiabatic index, implying world model. However, one can also see that the gas dwarf model
more degrees of freedom for a species, reduces the size of a convec- over-predicts the number of detected young sub-Neptunes at short
tive envelope for a fixed entropy. This effect, however, is smaller in orbital periods ≲ 10 days when compared to the TESS data. One
magnitude than that of changing the bulk mean molecular weight in plausible explanation for this discrepancy is an under-prediction of
the envelope. For our simplistic models, we have used the opacities planets that are “born rocky” within the gas dwarf model. Rogers &
for the appropriate mean molecular weight and H2 /He/H2 O mixture, Owen (2021) used CKS data to infer that ≲ 20% of planets should
thus maintaining self-consistency within the model. Nevertheless, be born as super-Earths by the end of disc dispersal, which may
future studies should include full opacity and equation of state ta- be explained with a combination of processes, including inefficient
bles, which account for intricate changes in chemical properties as gas accretion, boil-off and core-powered mass loss (e.g. Lee et al.
the planet evolves. Improvement to model sophistication in this man- 2022; Owen & Schlichting 2024; Rogers et al. 2024a). Explicitly
ner will help in refining the difference in the size of gas dwarfs and modelling these processes will likely reduce the number of predicted
water worlds at young ages. sub-Neptunes at short orbital periods at early times, and thus allow
If one wishes to exploit the effects of bulk mean molecular the model to better reproduce the observations. We leave this increase
weight, an important observational strategy is to observe planets in model complexity for future work. An alternative explanation for

MNRAS 000, 1–12 (2025)

10 J. G. Rogers
gas dwarf models. These studies have been limited, however, by
various observational effects, including the heightened activity of
pre-main sequence host stars and the requirement of high-resolution
imaging to resolve stars within their dense clusters. Recall from
Section 5.1 that the most informative planet detections are those
from the youngest stellar clusters, which are limited in number and
magnitude for missions such as TESS.
Mission concepts, such as the Early eVolution Explorer (EVE), aim
to perform dedicated transit surveys of young stellar clusters (Howard
et al. 2024). With the increase in planet yield from this style of survey,
one could also investigate mixture models in which planets from
both populations are present. This statistical analysis lends itself to
a Bayesian framework, allowing one to infer the appropriate fraction
of each model that best fits the data.
As shown in Section 3, the radius of an inflated planet is insensi-
tive to planet mass, which conveniently means that well-characterised
planet masses are not necessary to make inferences from population
models. Recent work has also shown that most young transiting
exoplanets are very unlikely to be progenitor hot-Jupiters (Karalis
et al. 2024). Nevertheless, transit spectroscopy of young planet at-
mospheres from JWST/HST can constrain a planet’s mass from its
scale height (e.g. Barat et al. 2024; Thao et al. 2024). Atmospheric
characterisation may, in fact, be a more straightforward route to a
young planet’s mass than conventional methods, such as radial ve-
Figure 5. The radius evolution of exoplanet atmospheres for different Rosse- locity measurements, which are notoriously challenging to obtain
land mean gas opacities of 𝜅 = 0.01, 0.1, 1.0 and 10.0 cm2 g −1 in blue, green, around young active stars. This also comes with the obvious benefit
orange and red, respectively. All planets have a 5𝑀⊕ core, with an Earth-like
of constraints on atmospheric chemistry at early times. EVE thus
core composition of 25% iron, 75% silicate. Each model has an atmospheric
mass fraction of 𝑋 ≡ 𝑀atm /𝑀c = 0.01 and a constant equilibrium temper-
also serves as a path-finder for suitable young planet targets for at-
ature of 𝑇eq = 800 K. All planets are initialised with a radiative-convective mospheric characterisation with JWST and future missions such as
boundary set at 5𝑅 ⊕ . ARIEL and HWO.

the lack of detected TESS planets ≲ 10 days is that the orbital period 5.4 Theoretical uncertainties
distribution of sub-Neptunes evolves during the post-disc era. Inward In this work, we have utilised simple radiative-convective equilibrium
migration of sub-Neptunes and super-Earths could be caused by the models to demonstrate that the bulk mean molecular weight of an
interaction of their escaping atmospheres with stellar winds, resulting envelope alters its size most noticeably at young ages. However, these
in a net angular momentum loss for the planet (Hanf et al. 2025), or models have neglected some important physical effects, which should
due to tidal migration (e.g. Lee & Owen 2025). be closely examined in future work.
One way to determine which of these mechanisms is responsible The models in this study have assumed a constant mean molecular
for the lack of observed young sub-Neptunes at short orbital periods is weight as a function of envelope radius. As discussed in Section 3,
to also measure the orbital period distribution of super-Earths. If there recent theoretical work has shown that planet envelopes chemically
is an abundance of super-Earths at short orbital periods at early times, react with their interiors, which in turn changes the mean molecular
it is likely that inefficient gas accretion and/or atmospheric escape weight within the atmosphere (e.g. Kite & Schaefer 2021; Schlichting
produced this population during the disc phase. Characterising the & Young 2022). Mean molecular weight gradients can form, which
young super-Earth population can also be used to measure the ratio may inhibit convection in the deep interior and thus change the en-
of super-Earths to sub-Neptunes as a function of time. One clear velope structure and evolution (e.g. Leconte et al. 2017; Misener
prediction of the gas dwarf model is that this ratio should increase & Schlichting 2022). The effects of phase equilibria and miscibility
with time since sub-Neptunes are being stripped of their atmospheres between various species such as H2 , H2 O, and silicate melt also re-
to become super-Earths (e.g. Gupta & Schlichting 2020; Rogers & main largely unexplored (e.g. Mousis et al. 2020; Vazan et al. 2022;
Owen 2021). However, probing this area of parameter space (𝑅p ≲ Pierrehumbert 2023; Gupta et al. 2024; Young et al. 2024). Cou-
2𝑅 ⊕ ) around young stars is technically challenging for instruments pled models of planet interiors and envelopes must be considered,
such as TESS, and will require next-generation transit observations. particularly for young planets due to their higher entropies and thus
active geochemistry. As discussed in Section 4, the gas dwarf popula-
tion model presented in this work slightly under-predicts the number
5.3 An observational route forward
of large sub-Neptunes ∼ 10 R ⊕ at young ages. More sophisticated
More detections of young transiting exoplanets are required to answer evolution modelling of such effects may help in interpreting this
science questions regarding planet formation and evolution at the under-prediction.
population level. Previous studies have shown that the young planet Beyond the interior and envelope structure of a planet, more atten-
population is inconsistent with the evolved Kepler population at the tion should be given to understanding the initial conditions of planet
∼ 1 − 3𝜎 level (Fernandes et al. 2022; Christiansen et al. 2023; evolution. For gas dwarf models (as discussed in Section 2.1), studies
Vach et al. 2024; Fernandes et al. 2025), with a strong indication have shown that the rapid disc dispersal process in the inner regions
that young sub-Neptunes are more inflated, and more consistent with ≲ 1 AU induces an extreme mass loss event referred to as “boil-off”,

MNRAS 000, 1–12 (2025)

 On the road to the radius valley 11
in which a planet may lose ∼ 90 % of its atmospheric mass (e.g. Owen spectroscopy may provide a clearer path to determining the masses
& Wu 2016; Ginzburg et al. 2016; Rogers et al. 2024a; Tang et al. of young planets as shown in recent works (e.g. Thao et al. 2024).
2024). These studies have shown that the initial cooling timescale
• There is tentative evidence that the population of observed tran-
of a gas dwarf post-disc dispersal may be longer than its age, which
siting planets from TESS (Vach et al. 2024) is inconsistent with a
is caused by extreme cooling of the envelope during boil-off (Owen
simple water world model (Figure 4). However, this sample is small
2020). In this study, we accounted for boil-off by enforcing initial
and, thus, the results may be driven by small-number statistics. We
planet sizes not exceeding 0.1𝑅B (see Owen & Wu 2016). Neverthe-
encourage future statistical comparisons to place emphasis on ac-
less, in future studies, the initial conditions of sub-Neptunes under
counting for this effect, as well as other survey biases in forward
various formation and composition models should be investigated.
models, for a fair comparison.
In addition, more work is needed to understand the mass loss of high
mean molecular weight species in exoplanet atmospheres. This is Inferring the composition of sub-Neptunes at the population level
a complex challenge due to various species-dependant heating and will inevitably rely on theoretical models. Future studies should place
cooling processes, as well as drag and diffusion. emphasis on understanding the implications of chemical interactions
Finally, the key to statistically inferring properties from a popula- between the interior and envelopes of such planets, as well as possible
tion of young transiting exoplanets is to use self-consistent structure initial conditions in their post-disc phase under a range of formation
and evolution models that can reproduce the properties of evolved scenarios.
exoplanets at ∼ Gyr epochs. These must also be computationally
fast in order to run large quantities of population-level models. The
balance between model complexity and speed may be aided with
ACKNOWLEDGEMENTS
machine-learnt emulators (Rogers et al. 2023a).
I would like to thank the anonymous referee for valuable comments
that helped improve this paper. I would also like to thank Sydney Vach
for sharing the TESS planet sample used in this study, as well as Jen-
6 CONCLUSIONS
nifer Burt, Neal Turner, Andrew Mann, Madyson Barber, George
In this paper, we have discussed the effects of bulk envelope mean Zhou, Eric Gaidos, Hilke Schlichting and James Owen for discus-
molecular weight on the early evolution of close-in exoplanets. The sions that helped improve this paper. JGR gratefully acknowledges
lower the bulk mean molecular weight, the larger the radius of a support from the Kavli Foundation.
planet at early times, which provides an observational window into
their composition and, thus, formation pathway. We have focussed
on sub-Neptunes, which can currently be explained with the “gas DATA AVAILABILITY
dwarf” hypothesis (rocky interiors and extended H2 /He dominated
envelopes) or the “water world” hypothesis (rocky interiors with All models will be made available upon reasonable request.
several tens of per cent in water mass fraction). A population of young
transiting exoplanets will allow for model comparison to determine
which of these models best describes reality. Our conclusions are as REFERENCES
follows:
Aguichine A., Batalha N., Fortney J. J., Nettelmann N., Owen J. E., Kempton
• The observed size of a young sub-Neptune is sensitive to the E. M. R., 2024, arXiv e-prints, p. arXiv:2412.17945
bulk mean molecular weight of its envelope. The lower the mean Anders E., Grevesse N., 1989, Geochimica Cosmochimica Acta, 53, 197
molecular weight, the larger the scale height, and thus the larger the Barat S., et al., 2024, Nature Astronomy, 8, 899
radius inflation. However, as planets cool, they contract. As a result, Barber M. G., et al., 2024, Nature, 533, 221
Batalha N. M., et al., 2013, ApJS, 204, 24
the range in planet radii as a function of mean molecular weight is
Bean J. L., Raymond S. N., Owen J. E., 2021, Journal of Geophysical Research
greatest at earliest times (Figure 1).
(Planets), 126, e06639
• The radii of evolved sub-Neptunes can be explained by two Benneke B., et al., 2024, arXiv e-prints, p. arXiv:2403.03325
interior compositions: "gas dwarfs" with low global water mass frac- Berger T. A., Huber D., Gaidos E., van Saders J. L., Weiss L. M., 2020, arXiv
tions and low envelope mean molecular weights in the approximate e-prints, p. arXiv:2005.14671
Burn R., Mordasini C., Mishra L., Haldemann J., Venturini J., Emsenhuber
range of 𝜇 ∼ 2.35 − 5 amu, or "water worlds" with high water mass
A., Henning T., 2024a, arXiv e-prints, p. arXiv:2401.04380
fractions and high envelope mean molecular weights in the range of Burn R., Bali K., Dorn C., Luque R., Grimm S. L., 2024b, arXiv e-prints, p.
𝜇 ∼ 6 − 18 amu. The radii of sub-Neptunes under these formation arXiv:2411.16879
pathways converge after ∼ Gyrs of evolution. However, their radii Caldiroli A., Haardt F., Gallo E., Spinelli R., Malsky I., Rauscher E., 2022,
diverge significantly at early ages ≲ 100 Myrs (Figure 2). A&A, 663, A122
Chen J., Kipping D., 2017, ApJ, 834, 17
• The prediction of different planet radii for the gas dwarf and
Christiansen J. L., et al., 2023, AJ, 166, 248
water world hypotheses can be exploited with a population of young Davenport B., et al., 2025, arXiv e-prints, p. arXiv:2501.01498
transiting exoplanets in order to probe mean molecular weight and, Dorn C., Lichtenberg T., 2021, ApJ, 922, L4
hence, sub-Neptune composition. We show that surveys of young Dos Santos L. A., 2023, IAU Symposium, 370, 56
planets should target the youngest observable stellar clusters to max- Ercolano B., Clarke C. J., Hall A. C., 2011, MNRAS, 410, 671
imise this difference to distinguish between formation channels. Erkaev N. V., Kulikov Y. N., Lammer H., Selsis F., Langmayr D., Jaritz G. F.,
Biernat H. K., 2007, A&A, 472, 329
• For young, puffy sub-Neptunes, planet size is relatively insen- Fernandes R. B., et al., 2022, arXiv e-prints, p. arXiv:2206.03989
sitive to planet mass. This removes the necessity of characterising Fernandes R. B., et al., 2025, arXiv e-prints, p. arXiv:2503.10856
planet masses to infer bulk composition, which is notoriously difficult Fortney J. J., Marley M. S., Barnes J. W., 2007, ApJ, 659, 1661
for planets orbiting young, active stars. Nevertheless, atmospheric Fressin F., et al., 2013, ApJ, 766, 81

MNRAS 000, 1–12 (2025)

12 J. G. Rogers
Fulton B. J., Petigura E. A., 2018, AJ, 156, 264 Rogers J. G., Gupta A., Owen J. E., Schlichting H. E., 2021, MNRAS, 508,
Fulton B. J., et al., 2017, AJ, 154, 109 5886
Ginzburg S., Schlichting H. E., Sari R., 2016, ApJ, 825, 29 Rogers J. G., Janó Muñoz C., Owen J. E., Makinen T. L., 2023a, MNRAS,
Ginzburg S., Schlichting H. E., Sari R., 2018, MNRAS, 476, 759 519, 6028
Guillot T., 2010, A&A, 520, A27 Rogers J. G., Schlichting H. E., Owen J. E., 2023b, ApJ, 947, L19
Gupta A., Schlichting H. E., 2019, MNRAS, 487, 24 Rogers J. G., Owen J. E., Schlichting H. E., 2024a, MNRAS, 529, 2716
Gupta A., Schlichting H. E., 2020, MNRAS, 493, 792 Rogers J. G., Schlichting H. E., Young E. D., 2024b, ApJ, 970, 47
Gupta A., Stixrude L., Schlichting H. E., 2024, arXiv e-prints, p. Rogers J. G., Dorn C., Aditya Raj V., Schlichting H. E., Young E. D., 2025,
arXiv:2407.04685 ApJ, 979, 79
Hanf B., Kincaid W., Schlichting H., Cappiello L., Tamayo D., 2025, AJ, 169, Santos N. C., et al., 2004, A&A, 426, L19
19 Schlichting H. E., Young E. D., 2022, PSJ, 3, 127
Ho C. S. K., Van Eylen V., 2023, MNRAS, 519, 4056 Schulik M., Booth R., 2022, arXiv e-prints, p. arXiv:2207.07144
Ho C. S. K., Rogers J. G., Van Eylen V., Owen J. E., Schlichting H. E., 2024, Szurgot M., 2015, in 46th Lunar and Planetary Science Conference.
arXiv e-prints, p. arXiv:2401.12378 Tang Y., Fortney J. J., Murray-Clay R., 2024, arXiv e-prints, p.
Howard W. S., et al., 2024, arXiv e-prints, p. arXiv:2411.08092 arXiv:2410.08577
Johnstone C. P., Bartel M., Güdel M., 2021, A&A, 649, A96 Thao P. C., et al., 2024, AJ, 168, 297
Karalis A., Lee E. J., Thorngren D. P., 2024, arXiv e-prints, p. Thorngren D., Fortney J. J., 2019, ApJ, 874, L31
arXiv:2408.16793 Vach S., et al., 2024, AJ, 167, 210
Kempton E. M. R., Lessard M., Malik M., Rogers L. A., Futrowsky K. E., Ih Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., Fulton
J., Marounina N., Muñoz-Romero C. E., 2023, ApJ, 953, 57 B. J., Petigura E., Snellen I., 2018, MNRAS, 479, 4786
Kenyon S. J., Hartmann L., 1995, ApJS, 101, 117 Vazan A., Sari R., Kessel R., 2022, ApJ, 926, 150
Kippenhahn R., Weigert A., Weiss A., 2012, Stellar Structure and Evolu- Venturini J., Alibert Y., Benz W., 2016, A&A, 596, A90
tion. Astronomy and Astrophysics Library, Springer Berlin Heidelberg, Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6
https://books.google.co.uk/books?id=wdSFB4B_pMUC Wolfgang A., Rogers L. A., Ford E. B., 2016, ApJ, 825, 19
Wu Y., Lithwick Y., 2013, ApJ, 772, 74
Kite E. S., Schaefer L., 2021, ApJ, 909, L22
Young E. D., Stixrude L., Rogers J. G., Schlichting H. E., Marcum S. P., 2024,
Koch D. G., et al., 2010, ApJ, 713, L79
arXiv e-prints, p. arXiv:2408.11321
Koepferl C. M., Ercolano B., Dale J., Teixeira P. S., Ratzka T., Spezzi L.,
Zeng L., et al., 2019, Proceedings of the National Academy of Science, 116,
2013, MNRAS, 428, 3327
9723
Kubyshkina D., Fossati L., 2022, arXiv e-prints, p. arXiv:2211.10166
Lecavelier Des Etangs A., 2007, A&A, 461, 1185
Leconte J., Selsis F., Hersant F., Guillot T., 2017, A&A, 598, A98 This paper has been typeset from a TEX/LATEX file prepared by the author.
Lee E. J., Owen J. E., 2025, ApJ, 980, L40
Lee E. J., Chiang E., Ormel C. W., 2014, ApJ, 797, 95
Lee E. J., Karalis A., Thorngren D. P., 2022, arXiv e-prints, p.
arXiv:2201.09898
Lopez E. D., 2017, MNRAS, 472, 245
Lopez E. D., Fortney J. J., 2013, ApJ, 776, 2
Lopez E. D., Fortney J. J., 2014, ApJ, 792, 1
Luo H., Dorn C., Deng J., 2024, arXiv e-prints, p. arXiv:2401.16394
Luque R., Pallé E., 2022, Science, 377, 1211
Madhusudhan N., Sarkar S., Constantinou S., Holmberg M., Piette A. A. A.,
Moses J. I., 2023, ApJ, 956, L13
Misener W., Schlichting H. E., 2021, MNRAS, 503, 5658
Misener W., Schlichting H. E., 2022, MNRAS, 514, 6025
Misener W., Schlichting H. E., Young E. D., 2023, MNRAS, 524, 981
Misener W., Schulik M., Schlichting H. E., Owen J. E., 2024, arXiv e-prints,
p. arXiv:2405.15221
Mordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139
Mousis O., Deleuil M., Aguichine A., Marcq E., Naar J., Aguirre L. A.,
Brugger B., Gonçalves T., 2020, ApJ, 896, L22
Owen J. E., 2020, MNRAS, 498, 5030
Owen J. E., Jackson A. P., 2012, MNRAS, 425, 2931
Owen J. E., Schlichting H. E., 2024, MNRAS, 528, 1615
Owen J. E., Wu Y., 2013, ApJ, 775, 105
Owen J. E., Wu Y., 2016, ApJ, 817, 107
Owen J. E., Wu Y., 2017, ApJ, 847, 29
Petigura E. A., Howard A. W., Marcy G. W., 2013, Proceedings of the National
Academy of Science, 110, 19273
Petigura E. A., et al., 2018, AJ, 155, 89
Petigura E. A., et al., 2022, AJ, 163, 179
Piaulet-Ghorayeb C., et al., 2024, ApJ, 974, L10
Pierrehumbert R. T., 2023, ApJ, 944, 20
Raymond S. N., Boulet T., Izidoro A., Esteves L., Bitsch B., 2018, MNRAS,
479, L81
Rivera E. J., et al., 2005, ApJ, 634, 625
Rogers L. A., 2015, ApJ, 801, 41
Rogers J. G., Owen J. E., 2021, MNRAS, 503, 1526

MNRAS 000, 1–12 (2025)