Assessing Core-powered Mass Loss in the Context of Early Boil-off: Minimal Long-lived Mass Loss for the Sub-Neptune Population

Planet formation is an energetic process, leading to a hot initial condition and subsequent cooling of the planetary interior. The cooling of the H/He envelope allows gravitational binding energy and internal energy to be released as thermal energy, and as a result, the planet's interior contracts. Such thermal contraction is critical for modeling atmospheric escape, as the planet radius controls the potential well that H/He material is lifted from. Similar to previous planet evolution codes (J. J. Fortney et al. 2007; E. D. Lopez & J. J. Fortney 2014), the planet's gaseous H/He is assumed to consist of an envelope that is adiabatic and isentropic owing to the short mixing timescale of convection, with a radiative atmosphere on top. Given the entropy and the envelope mass for the envelope, the interior thermal state of a planet is completely defined by the equation of state, which we take from G. Chabrier et al. (2019), and hydrostatic equilibrium.

Between each time step, to evolve the interior, we track the energy fluxes that heat the envelope from below and cool the envelope at the top, which change the interior thermal states. The isentropic envelope links the change of specific entropy Δs of the envelope adiabat to the net heat transfer ΔQ = LΔt = Δs∫Tdm to the envelope within the time interval Δt, which gives

$$\begin{eqnarray}\begin{array}{rcl}-L & = & \displaystyle \frac{{ds}}{{dt}}{\int }_{{M}_{c}}^{{M}_{c}+{M}_{\mathrm{env}}}T\,{dm}\\ & = & -{L}_{\mathrm{env}}+{L}_{\mathrm{core}}\\ & = & -\max ({L}_{\mathrm{int}},{L}_{\mathrm{ad}})+{L}_{\mathrm{radio}}-{c}_{v}{M}_{c}\displaystyle \frac{{{dT}}_{c}}{{dt}},\end{array}\end{eqnarray} \tag{ 1 }$$

where M_c and M_env are the core mass and the envelope mass, respectively, and T is the temperature of each mass shell dm. On the right-hand side, we sum up each luminosity component to get the total envelope luminosity L, where the envelope cooling L_env at the top is set by the greater of either the radiative intrinsic luminosity L_int or advective luminosity L_ad, the energy transport from the bulk flow driven by atmospheric escape. The envelope heating from below is due to the total core luminosity L_core from both radiogenic heating L_radio and the core heat capacity.

The intrinsic luminosity L_int, the net cooling rate, is the total energy that a planet's interior radiates per unit time through its radiative atmosphere to space. To evaluate the intrinsic luminosity, we utilize a grid of one-dimensional atmosphere models for solar metallicity, as described in J. J. Fortney et al. (2007), over a range of surface gravities, incident bolometric fluxes, and interior temperatures. As a planet loses mass through atmospheric escape, the hydrodynamic wind is capable of transporting internal thermal energy out of the interior. This effect is generally minor for a long-lived mass loss, i.e., photoevaporation and core-powered mass loss. However, when a planet is rapidly blowing off its H/He envelope in the boil-off phase, the mass-loss rate is a few orders of magnitude greater than that at a later age, making such an energy advection, L_ad, potentially important. J. E. Owen & Y. Wu (2016) argue that the energy advection from boil-off leads to an interior cooling that is much more significant than the radiative cooling, with the advective cooling rate L_ad estimated by

$$\begin{eqnarray}&&{L}_{\mathrm{ad}}=\displaystyle \frac{\gamma }{\gamma -1}\dot{M}{c}_{s}^{2},\end{eqnarray} \tag{ 2 }$$

where γ is the adiabatic index, $\dot{M}$ is the mass-loss rate, and c_s is the sound speed. This cooling effect is evaluated in this work.

For simplicity, the core is assumed to be isothermal with the temperature found at the bottom of the envelope, essentially meaning that it is an efficient conductor. As a planet cools off, the core temperature correspondingly drops, which simultaneously releases thermal energy via the core's heat capacity. For the specific heat capacity at a constant volume of the core, we use c_v = 0.75 J K⁻¹ g⁻¹. To stabilize our numerical calculation, we assess the time derivative of core temperature dT_c /dt based on the change of core temperature over the last five time steps. For the same reason, the planets are not allowed to increase specific entropy over time, such that ${L}_{\mathrm{core}}\leqslant {L}_{\mathrm{env}}$. The consequence of this simplification is evaluated in Sections 3.1.5 and 6.2. In terms of L_radio from radioactivity, the dominant decaying isotopes are ²³⁵U, ²³⁸U, ⁴⁰K, and ²³²Th. The abundances and radioactive powers of these isotopes at early ages are derived based on their half-lives and the meteoritic abundances at the current solar age (E. Anders & N. Grevesse 1989; N. Nettelmann et al. 2011).

Based on the $\tfrac{{ds}}{{dt}}$ calculated from Equation (1) at each time step, we evolve the envelope specific entropy using a fifth-order ordinary differential equation (ODE) solver. The resulting specific entropy at a new time step automatically defines the amount of the total energy (gravitational energy and internal energy) that is allowed to be released through interior cooling and therefore the rate of thermal contraction. The envelope mass M_env also controls both the cooling process and the hydrostatic interior structure, as seen on the left-hand side of the equation. In the presence of mass loss, given the mass-loss rate $\dot{M}$, we evolve the total planetary mass M_p in a similar manner simultaneously, as discussed below.

The radiative atmosphere is a static layer whose radial structure passively evolves with time on top of the adiabatic interior as a planet contracts. We focus on planets located close enough to their host stars that the radiative−convective boundary (RCB) separating the adiabatic interior from the radiative atmosphere is set by the incident bolometric flux from the star. At the RCB, the planet's intrinsic luminosity is equal to the stellar energy deposited per time (i.e., stellar flux reduced by inward diffusion of photons since the optical depth at the RCB is generally large). Such radiative atmospheres typically exhibit two roughly isothermal regions—an outer region optically thin to both incident optical radiation and outgoing infrared radiation, and a deeper region that is optically thick to outgoing infrared (and, in its deepest parts, to incoming optical radiation as well; T. Guillot 2010). The radiative atmosphere's temperature–pressure (T−P) profile is thus primarily set by stellar heating. Stellar heating, especially that from M dwarfs and in XUV wavelengths, can decrease over time at young ages. However, the variability in bolometric flux from a Sun-like star is generally modest at the ages relevant to the boil-off phase, ∼3–10 Myr after the formation (I. Baraffe et al. 2015). Since we only focus on the history of mass loss and thermal evolution resulting from stellar bolometric flux, independent of stellar spectral types and XUV-driven escape, the time variability of stellar heating is ignored. This assumption only quantitatively affects the time-integrated mass loss by a small fraction.

For gas opacity to thermal radiation κ_th, the transition from optically thin to thick occurs where optical depth τ ∼ κ_th ρ H ∼ 1, for atmospheric density ρ and scale height H. In a hydrostatic atmosphere, ρ H ∝ Pr², so that the pressure where optical depth transitions occur depends only on radius r, which, after the initial blow-off period, typically changes modestly over the planet's evolution, leading to a roughly constant T−P profile in the radiative atmosphere. To capture this behavior, we do not perform a full radiative transfer calculation but rather generate a T−P profile using the widely used analytical two-stream radiative transfer model of T. Guillot (2010). The thermal opacity κ_th is chosen to be a constant 0.02 cm² g⁻¹, representative of cool H/He solar-composition atmospheres at P ∼ 1 bar (R. S. Freedman et al. 2014). In terms of the opacity to visible photons κ_ν , we compare the T−P profile of the analytical two-stream model to that from full radiative−convective equilibrium transfer models (J. J. Fortney et al. 2008, 2013) to find the best-fit opacity values (in the range of 0.0001−0.006 cm² g⁻¹) as a function of incident bolometric flux in order to well approximate the T−P profile from the full calculation. The values are found to be only weakly dependent on the surface gravity and the intrinsic luminosity from the interior. Beginning from the outer edge of the envelope, we integrate using hydrostatic equilibrium, taking into account the variation of gravity with radius.

The thickness of the radiative atmosphere is not negligible for low-mass planets because of their low gravity, which decreases with altitude and leads to an outwardly increasing scale height. Interestingly, the radiative atmosphere can contribute to a major fraction of the optical radius of a sub-Neptune, making it important for population studies. Our planetary radius is defined at 20 mbar, a typical pressure level for optical transits.

The bottom of the radiative atmosphere is separated from what we term the H/He "envelope" (Section 2.1) at the RCB. Our numerical model estimates the location of the RCB by finding the intersection between the P−T profiles of the adiabatic envelope and the radiative atmosphere. Compared to an RCB pressure estimated from full radiative transfer calculations, our approach only quantitatively shifts the RCB location in radius by a negligible fraction. The RCB pressure for a sub-Neptune varies dramatically with age, from typically at ∼10 bars (young ages) to 10 kbar (Gyr+ ages). As a sub-Neptune evolves, its interior cools off and the intrinsic luminosity declines. Consequently, the mismatch between the large stellar bolometric flux and the weaker intrinsic flux becomes ever larger, leading to a deeper RCB. Figure 1 illustrates this evolution of the RCB location (crosses) over time.

Zoom In Zoom Out Reset image size

Furthermore, we calculate the total atmospheric mass in the static radiative layer at each time step by integrating the mass shells in the radiative layer. The radiative atmospheric mass does not participate in thermal evolution and is therefore excluded from the H/He envelope mass to better assess the thermal evolution process. We find that the radiative atmospheric mass fraction is typically at most a few hundredths of a percent of the planetary mass during most of the evolution time, contributing a minor difference between the thermal evolution with and without it. However, an exception happens at the early stage of boil-off. The radiative atmospheric mass can constitute up to 20%–50% of the convective envelope mass for a young and inflated planet, thus playing a role in prolonging the thermal contraction timescale. This substantial radiative atmospheric mass results from the slowly varying density–radius profile ($\rho \sim {\rho }_{\mathrm{rcb}}\exp (-r/H)$ for a roughly isothermal atmosphere), which remains high (∼ρ_rcb, the density at the RCB) in the radiative atmosphere owing to the large scale height H that is comparable to the planetary radius r above the RCB. Note that the quantitative details remain uncertain here, especially if additional factors such as the opacity contribution from dust in both the atmosphere and the debris disk are considered.

Another exception occurs when the planetary envelope is about to be completely stripped. When a planet loses mass from the deepest part of the radiative atmosphere, right above the RCB, the mass there has to be steadily resupplied from the envelope to sustain an outflow (S. Ginzburg et al. 2016). By the time that the convective zone is depleted by the mass loss, the radiative atmospheric mass (∼0.01% of the planetary mass) becomes most of the total H/He mass. In this case, thermal evolution of the envelope starts to cease and the planet transitions into an envelope-free super-Earth. We terminate thermal evolution and mass loss once the envelope mass becomes negligible compared to the radiative atmospheric mass, and after that, we treat the planet as a bare core super-Earth, assuming that its radiative atmosphere is then completely lost in a sufficiently short period.

We evolve the total planetary mass given the mass-loss rate by solving the ODE problem ${{dM}}_{p}/{dt}=\dot{M}$. At each time step, we estimate the envelope mass M_env by subtracting the constant core mass M_c and the time-dependent radiative atmospheric mass M_atm from the total mass M_p :

$$\begin{eqnarray}&&{M}_{\mathrm{env}}={M}_{p}-{M}_{c}-{M}_{\mathrm{atm}}.\end{eqnarray} \tag{ 3 }$$

The atmospheric mass, M_atm, implicitly depends on M_env because the envelope mass M_env provides the lower boundary condition for the radiative atmospheric mass calculation. Since the atmospheric mass changes slowly compared to the mass-loss rate, to improve the computation efficiency we take the value of the atmospheric mass from the last converged solution and treat it as a constant within the next time step to avoid needing to solve Equation (3) self-consistently using an iterative method at each structure calculation (there are many within one time step). A full self-consistent calculation of M_atm is employed at the first time step of the evolution.

2.3.1. Parker Wind

Physically speaking, there is no fundamental difference between the nature of boil-off and core-powered escape. Both processes are directly caused by the large-radius H/He envelope and atmosphere in the absence of sufficient ambient confinement pressure, which cannot maintain a hydrostatic equilibrium and therefore leads to a Parker wind outflow. Similar to Bondi accretion in reverse, the Parker wind is a hydrodynamic process that occurs when the protoplanetary disk has dissipated. As long as sufficient energy is available to resupply the H/He material in the atmosphere from below (we evaluate this assumption in Section 3.1.2), the stellar radiation dictates the energy budget for the outflow and therefore the temperature profile in the radiative atmosphere, where the wind advection occurs. For a planet irradiated by a star with an incident flux F, its radiative atmosphere is roughly isothermal with the temperature set by the equilibrium temperature ${T}_{\mathrm{eq}}={\left(\tfrac{F}{4\sigma }\right)}^{1/4}$, where σ is the Stefan–Boltzmann constant. For simplicity, we model the hydrodynamic wind using an isothermal Parker wind (E. N. Parker 1958) beyond the RCB to assess the mass-loss rate. This avoids needing an implementation of a numerical hydrodynamic and radiative transfer calculation. The advective−convective boundary is chosen to be the same as the RCB, which is justified later in Section 3.1.2. The wind structure calculation is in parallel with the radiative atmosphere calculation, where the hydrostatic calculation is utilized to assess the optical transit radius (which is not an important variable until boil-off ends), and the steady-state Parker wind determines the mass-loss rate.

The physical structure of an isothermal Parker wind is given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{{c}_{s}^{2}}{{v}^{2}}\right)\displaystyle \frac{d}{{dr}}({v}^{2})=-\displaystyle \frac{{{GM}}_{p}}{{r}^{2}}\left(1-\displaystyle \frac{2{c}_{s}^{2}r}{{{GM}}_{p}}\right),\end{eqnarray} \tag{ 4 }$$

where ${c}_{s}=\sqrt{{{kT}}_{\mathrm{eq}}/\mu }$ is the isothermal sound speed, k is the Boltzmann constant, μ = 2.35m_H is the mean molecular weight of the wind, m_H is the mass of atomic hydrogen, v is the wind velocity, and r is the radius. At a radius known as the sonic point, a singularity occurs where the wind velocity reaches the sound speed and both sides of the equation vanish:

$$\begin{eqnarray}&&{R}_{s}=\displaystyle \frac{{{GM}}_{p}}{2{c}_{s}^{2}}.\end{eqnarray} \tag{ 5 }$$

At the wind base, density and pressure transition continuously between the two layers of the steady-state hydrodynamic wind and the hydrostatic adiabatic envelope. This gives us the mass-loss rate:

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{parker}}}=4\pi {R}_{{\rm{rcb}}}^{2}{\rho }_{{\rm{rcb}}}{v}_{{\rm{rcb}}}=4\pi {r}^{2}\rho v,\end{eqnarray} \tag{ 6 }$$

from which the density profile ρ for the wind as a function of radius r is obtained, where R_rcb, ρ_rcb, and v_rcb are the planet radius, density, and wind velocity at the RCB, respectively.

2.3.2. Transition Out of Boil-off

Boil-off is a short-lived hydrodynamic process, hence the name, with a mass-loss timescale shorter than 1 Myr (J. E. Owen & Y. Wu 2016) owing to the sudden change of the ambient confinement pressure from a protoplanetary disk. The mass-loss timescale is given by

$$\begin{eqnarray}&&{t}_{\dot{M}}=\displaystyle \frac{{M}_{\mathrm{env}}}{\dot{M}},\end{eqnarray} \tag{ 7 }$$

where $\dot{M}$ is the mass-loss rate of the isothermal Parker wind and M_env is the envelope mass. If the thermal contraction decreases the planetary size faster than the mass loss can substantially alter the envelope's mass, the wind velocity v at the RCB in Equations (4) and (6) declines quickly over time compared to the rate of mass evolution, and consequently, the planet has to stop blowing off. Therefore, we define the end of boil-off at the time when the mass-loss timescale becomes comparable to the Kelvin–Helmholtz contraction timescale t_cool:

$$\begin{eqnarray}&&{t}_{\dot{M}}={t}_{\mathrm{cool}},\end{eqnarray} \tag{ 8 }$$

where the Kelvin–Helmholtz contraction timescale is assessed with the total envelope luminosity L in Equation (1) that incorporates both the envelope cooling L_env and the core luminosity L_core:

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{{{GM}}_{c}{M}_{\mathrm{env}}}{\alpha {R}_{\mathrm{rcb}}L},\end{eqnarray} \tag{ 9 }$$

where α ≤ 1 is a dimensionless factor related to the mass concentration of the envelope. With our initial setup of our model, we find that the H/He envelope mass and therefore gravitational binding energy are very slightly inwardly concentrated at the very beginning of boil-off, yielding α = 0.4. The mass concentration gradually shifts outward to α = 0.8 for an old planet that has a lower H/He mass and specific entropy. During most of the boil-off phase, the mass and energy exhibit a central-to-outward distribution with α ≥ 0.5. In boil-off, the mass loss and thermal contraction are decoupled owing to the large timescale difference ${t}_{\dot{M}}\ll {t}_{\mathrm{cool}}$, leading to negligible thermal evolution and hence a minor decrease in entropy, compared to that in the later evolution. Therefore, the shrinkage in planetary size is due to the mass loss rather than thermal contraction.

After the transition time defined in Equations (8) and (9), the mass-loss rate declines and thermal contraction becomes the dominant physical effect in controlling the planetary size. The isothermal Parker wind transitions into the post-boil-off long-lived mass-loss phase. The long-lived mass-loss history is greatly dependent on the thermal contraction of the planet, in which core luminosity can potentially play an important role. In this case, the mass loss and thermal evolution are coupled. A representative combined evolution of planet radii, temperatures, and H/He mass fractions from our numerical sub-Neptune model is shown in Figure 1, with the boil-off phase shown with dashed lines and the long-lived mass-loss phase shown with solid lines.

2.3.3. Transition to XUV-driven Escape

As a planet loses mass to the Parker wind mass loss, its physical radius shrinks and the upper atmosphere becomes transparent to XUV photons. Therefore, the high-energy photons are able to penetrate deeply enough into the atmosphere to drive a more efficient wind with a higher temperature. Atmospheric escape starts to transition to XUV-driven when the optical depth to XUV photons, evaluated at the sonic point, is equal to unity (J. E. Owen & H. E. Schlichting 2024):

$$\begin{eqnarray}&&{\tau }_{{\rm{s}}}={\sigma }_{\nu 0}{{Hn}}_{0}\sim 1,\end{eqnarray} \tag{ 10 }$$

where σ_ν0 is the cross section for the absorption of XUV photons for hydrogen, H is the scale height of the upper atmosphere, and n₀ is the neutral hydrogen number density. Based on our numerical model, we find that such a transition typically happens when boil-off is about to be over, at around 1 Myr. For the purpose of our study, XUV-driven escape is not directly included in the mass loss from our numerical model.

The planetary radius dictates the intensity of boil-off. The radiative atmosphere of an inflated planet, with its large scale height, is substantially less bound to the interior, making the planet lose mass readily. Since planets with hotter internal thermal states possess larger radii at a given envelope mass fraction, we argue that the strength of boil-off is directly controlled by the envelope entropy. Physically speaking, because the final entropy dictates both the thermal contraction and mass-loss timescales in Equation (8), which are equal at the transition time, the final mass fraction is completely determined by the final entropy s_t . Given the negligible thermal evolution compared to mass loss during boil-off in Equation (8), the specific entropy remains nearly unchanged during boil-off. Therefore, the initial entropy s_i ∼ s_t largely determines the physical states at the end of boil-off and the beginning of the subsequent evolution.

However, what entropy a planet is born with is not well determined from previous work. To constrain the initial entropy at the beginning of boil-off, J. E. Owen & Y. Wu (2016) suggest that a model planet should not be allowed to thermally contract (due to cooling) faster than the disk dispersal time ∼10⁵ yr, the characteristic timescale for the formation of the gas-depleted inner hole; otherwise, it is always able to adjust to a new hydrostatic equilibrium during the disk dispersal, leading to no efficient boil-off. On the other hand, the planetary initial contraction time should be no faster than the disk lifetime 3–10 Myr (E. E. Mamajek 2009; U. Gorti et al. 2016), as in the opposite case the planet would shrink rapidly enough while the disk is present to allow for more gas accretion, and consequently the increased envelope mass would eventually slow down the contraction timescale.

For these reasons, in this work the initial entropy and radius are chosen such that the initial Kelvin–Helmholtz thermal contraction timescale in Equation (9) is comparable to the disk lifetime, chosen to be 5 Myr. Note that the core luminosity and advective cooling are initially 0, such that the total envelope luminosity is set by the intrinsic luminosity L = L_int. As the disk dispersal time of ∼10⁵ yr is observed to be substantially shorter than the disk lifetime, our physical choice ensures that the planet's contraction timescale at the onset of boil-off is long compared to the boil-off time and hence the process of boil-off is not sensitive to direct modeling of the disk dispersal. Our initial setup results in an initial RCB radius that is 10%–20% of the sonic radius defined in Equation (5) and an initial entropy of (9–11)k_b /atom depending on the bolometric flux and core mass.

Most of the previous core-powered escape work (S. Ginzburg et al. 2016, 2018; A. Gupta & H. E. Schlichting 2019; W. Misener & H. E. Schlichting 2021) originates from a single analytical model with similar mass-loss treatments. According to these authors, although the nature of both boil-off and core-powered mass loss is a Parker wind, core-powered mass loss differs from boil-off in the following three aspects: the energy source for mass re-equilibration, the mass-loss timescale (short-lived or long-lived), and whether or not the mass-loss rate is enhanced by core luminosity. We focus on assessing each of these three aspects in this section.

A primary challenge for modeling these processes regards the energy supply that overcomes the gravitational force, to continuously drive such winds. If the energy input is inadequate, the wind region rapidly cools off when the PdV work drains energy from the internal energy of the outflow, leading to a low wind temperature in Equation (4), and therefore the outflow slows down until the energy supply is sufficient to maintain the outflow. This is known as an energy-limited wind. Two layers where energy supply may potentially limit the mass-loss rate are invoked: the radiative atmosphere and adiabatic envelope, shown in Figure 2. The energy injection in the radiative zone is primarily from the bolometric luminosity and is assumed in these studies to be sufficient to sustain an outflow of the equilibrium temperature ∼T_eq. In Section 3.1.3, we show that this assumption is reasonable for most, though not all, planets. However, if mass cannot be replenished fast enough from within the convective zone, H/He material will eventually be depleted, yielding a low-density layer at the RCB, even if sufficient energy to power the wind is available in the radiative zone.

Zoom In Zoom Out Reset image size

Envelope re-equilibration was considered as the bottleneck for boil-off and core-powered mass loss in S. Ginzburg et al. (2016, 2018, hereafter GSS16 and GSS18). GSS18 assume that most of the H/He mass that replenishes the wind at the RCB starts to expand from as deep as the core−envelope surface, with the energy needed for the re-equilibration of the envelope completely from the interior cooling energy L_int. This requires a large amount of cooling energy to sustain the outflow, and thus the outflow can be largely energy limited. The total energy needed per second is estimated by the amount of gravitational energy to overcome:

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{loss},\mathrm{core}}=\displaystyle \frac{{{GM}}_{c}\dot{M}}{{R}_{c}},\end{eqnarray} \tag{ 11 }$$

where R_c is the core radius. Based on ${\dot{E}}_{\mathrm{loss},\mathrm{core}}$, these authors estimate an energy-limited mass-loss rate if the mass re-equilibration is insufficient:

$$\begin{eqnarray}&&{\dot{M}}_{e-\mathrm{lim},\mathrm{core}}=\displaystyle \frac{{L}_{{int}}}{{{GM}}_{c}/{R}_{c}}.\end{eqnarray} \tag{ 12 }$$

Once the intrinsic luminosity becomes more than the energy needed in Equation (11), they treat the wind as non-energy-limited, as in Equation (6), which is the maximum possible mass-loss rate given the sufficient energy and mass supply, known as what they call a Bondi-limited regime. They define core-powered mass loss as the later stage of Parker wind mass loss when core cooling L_core constitutes a large fraction of the envelope cooling, so ${L}_{\mathrm{int}}\sim {L}_{\mathrm{core}}$. In this case, the energy that liberates H/He mass from the interior is ultimately from the core, rather than the envelope itself. This happens to planets that have a larger core thermal energy reservoir than the gravitational binding energy. On the other hand, they refer to boil-off (what they call spontaneous mass loss) as the earlier stage when the envelope energy dominates over the core thermal energy. The above discussion provides a key justification for their mass-loss treatment and the necessity of distinguishing between boil-off and core-powered mass loss.

However, an important assumption made in their envelope re-equilibration argument is that the energy released from the envelope internal energy when H/He mass is lifted and cooled from the hot deep interior to the colder RCB is ignored in Equation (11). We suggest that their argument in Equations (11) and (12) would hold true only if the envelope were isothermal (in which case, the contribution from the internal energy can be ignored, corresponding to a steady-state isothermal envelope), but not for the adiabatic envelope assumed in their model. We revisit the envelope re-equilibration problem with an adiabatic envelope considered in the energy calculation, with particular attention to the internal energy exchange in Section 3.1.2. We identify a different and much physically narrower (thus requiring less energy to overcome) bottleneck region due to the deficiency of L_int in Section 3.1.3 (Figure 2).

3.1.1. A Fast Re-equilibration of Energy and Mass via Envelope Convection

We begin by verifying that the envelope remains adiabatic and quasi-hydrostatic during an energetic boil-off mass-loss phase. We calculate the re-equilibration timescales and the mass transport and internal energy transport propagating via convection within the envelope, shown in Figure 3. The top panel shows three timescales. These timescales are evaluated at the envelope surface (RCB), as we find that they are either greater than the timescales at the bottom of the envelope or comparable. For any perturbation at the envelope surface, from either removing envelope mass or thermal cooling, the pressure and density of the deep envelope correspondingly react and readjust to a new hydrostatic equilibrium on the dynamical timescale (dashed line with t_dyn = R_rcb/c_s ), which corresponds to the time that sound waves take to traverse the entire planetary radius (since the atmosphere is close to being in hydrostatic equilibrium). Under our initial physical conditions, the wind crossing timescale t_cross = R_rcb/v_rcb, where v_rcb is the wind speed at the surface, is always much longer than the dynamical timescale, given the small Mach number M = v_rcb/c_s ≪ 1 at the RCB, so that the density, pressure, and therefore radius of the envelope presumably respond instantaneously to the mass loss. During this process, the RCB pressure and temperature (and therefore density) physically remain unchanged, as both the outgoing intrinsic luminosity (primarily set by s) and the incident bolometric luminosity, which dictates the RCB conditions, are not affected owing to their long timescale to change. Consequently, the RCB has to penetrate to a deeper location in radius to react to the pressure decline at the original location.

Zoom In Zoom Out Reset image size

However, energy transport is limited by convection. It mixes H/He with different thermal states at different depths, which ultimately determines the temperature (entropy) profile and therefore the mass distribution. In boil-off, as the density and pressure of the envelope decrease, the envelope temperature also has to decrease in order to adjust to the lowered pressure, as shown in the dashed curve from the top-right panel of Figure 1. Note that the envelope specific entropy s barely changes owing to ${t}_{\dot{M}}\ll {t}_{\mathrm{cool}}$. This process eventually releases a large amount of the internal energy of the envelope, providing an important amount of energy to the re-equilibration process. This temperature decline is distinct from the interior thermal cooling that drives the thermal contraction and the change of specific entropy s, as during this process a negligible amount of net heat transfer happens to the envelope.

Such a temperature change and internal energy transport are via convection on a timescale set by the eddy diffusion timescale (dotted line). We estimate this diffusion timescale t_conv as the eddy turnover time (solid) H/v_conv from mixing-length theory times the number of turnovers for the energy and mass to be transported to the envelope surface, which we take to be ∼(ΔR/H)², where ${\rm{\Delta }}R={R}_{\mathrm{RCB}}-{R}_{\mathrm{core}}$, H is the pressure scale height of the envelope, and v_conv is the convective velocity. We calculate v_conv using Equation 3.7 from T. Guillot et al. (2004), with the mixing-length parameter α_m = 1 and the specific heat capacity at constant pressure c_p for a diatomic ideal gas. The convective energy flux F_conv is chosen to be the intrinsic cooling flux ${L}_{\mathrm{int}}/(4\pi {R}_{\mathrm{rcb}}^{2})$ evaluated at the RCB, and the profile parameter $\delta =-{(\partial \mathrm{ln}\rho /\partial \mathrm{ln}T)}_{P}$ is estimated using an adiabatic equation of state. The other physical quantities are computed from our numerical model.

The above discussion therefore gives ${t}_{\mathrm{dyn}}\ll {t}_{\mathrm{conv}}\ll {t}_{\dot{M}}\ll {t}_{\mathrm{cool}}$ even at a very early age shown in Figure 3 (top). Note that the shortest mass-loss timescale is found at the beginning of boil-off, ∼1 yr. Therefore, from the above discussion, we see that an adiabatic and hydrostatic envelope that is widely assumed in giant planet evolution models and sub-Neptune models with mild atmospheric escape including photoevaporation and core-powered mass loss still holds true even in the presence of vigorous boil-off, due to the fast mixing effect from envelope convection. This justifies our assumption made in Equation (1).

3.1.2. Lifting Mass from Convective Envelope

In this section we focus on whether there is enough internal energy transport to lift H/He mass from the deep envelope to supply the outflow at the RCB. Intuitively, since we have shown that the envelope remains adiabatic through the mass-loss process, it does not require energy input to sustain gas expansion within the envelope. Instead, the required energy for re-equilibration can be completely taken from the internal energy of the gas. The only additional energy required is the energy input to supply the kinetic energy of the wind, which we find is generally negligible in the envelope for an envelope radius R_rcb a few times smaller than the sonic radius R_s . Because it has important implications for the rate of core-powered mass loss, we spend the remainder of this section validating this intuition that after mass loss the convective envelope does not require external energy input to re-equilibrate on a new adiabat, even when that re-equilibration involves expansion of the envelope in the planet's potential well.

Within the envelope, convection is capable of transporting internal energy at a rate even greater than the wind advection is. To demonstrate this, in the bottom panel of Figure 3 we show that envelope convection advects an amount of internal energy (black; this does not equal the actual net amount of internal energy advected by the wind) much greater than the energy needed to overcome the gravitational force (gray), evaluated per radial distance at the RCB:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{conv}}\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{{kT}}{(\gamma -1)\mu }\right]\gg \displaystyle \frac{{{GM}}_{p}\dot{M}}{{r}^{2}}\ .\end{eqnarray} \tag{ 13 }$$

We find that the convective mass flux ${\dot{M}}_{\mathrm{conv}}=4\pi {r}^{2}\rho {v}_{\mathrm{conv}}$, evaluated at the RCB (black), corresponding to the minimal amount of mass flux that is required to transport the intrinsic radiation out of the envelope, is always a few orders of magnitude greater than the mass-loss rate $\dot{M}$ (gray) in the middle panel. It indicates that the redistribution of energy and mass through envelope convection is always sufficient to adjust the envelope adiabatically and hydrostatically against the perturbation generated from the mass loss, without needing an additional energy source, i.e., intrinsic luminosity, as opposed to the suggestion from GSS18. Additionally, this indicates that the wind advection only dominates the radiative zone rather than the envelope, as boil-off requires a mass re-equilibration much less than convection could yield, so we suggest that the advective−convective boundary coincides with the RCB.

So far, our re-equilibration arguments are valid for a steady-state convective mass/energy flux in the adiabatic envelope, but in reality the H/He mass has to be ultimately taken from each layer of the envelope (i.e., ∂ρ/∂t ≠ 0 at any radius in the envelope). In a more general case, we have to assess whether enough internal energy can be generated from the re-equilibration compared to the energy needed to lift the H/He material in each layer of the envelope. This involves the Euler energy equation for a fluid:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{t}}{\partial t}+\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}({h}_{t}+\phi )={h}_{t}\displaystyle \frac{\partial \rho }{\partial t}\ ,\end{eqnarray} \tag{ 14 }$$

where we have dropped the external heating and cooling terms because they are negligible in the convective envelope. The combined energy per volume E_t = ρ e + ρ v²/2 is the sum of kinetic energy and internal thermal energy, where e = c_v T is specific internal energy and c_v , ρ, v , and T are the specific heat capacity, density, velocity, and temperature, respectively. We have combined the specific enthalpy of the envelope, h = e + p/ρ, with the specific kinetic energy in h_t ≡ h + v²/2, and ϕ = − GM/r is the gravitational potential. Equation (14) may be reformulated in terms of the energy per volume E_tot ≡ ρ e + ρ ϕ, so that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{{tot}}}{\partial t}=-\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}B+B\displaystyle \frac{\partial \rho }{\partial t}\ ,\end{eqnarray} \tag{ 15 }$$

where we have dropped a term proportional to ∂ϕ/∂t since the gravitational potential is primarily determined by the core mass M_c and does not vary significantly within a mass-loss time step. Equation (15) implies the classic result that the Bernoulli constant B ≡ h_t + ϕ is conservative along streamlines (which in our problem are radial) for a steady-state flow (∂/∂t = 0). While our problem is not in steady state, the mass-loss timescale is substantially longer than the timescale over which the convective envelope re-equilibrates into a pseudo-steady state (Figure 3). We demonstrate in Appendix A that, as a result, the first term on the right-hand side of Equation (15) is small compared to the other two terms in the equation.

Equation (15) may be further simplified by noting that the Mach number of the fluid motion in the convective envelope ${ \mathcal M }\ll 1$, so that in the envelope the kinetic energy $\rho {v}^{2}/2\sim {{ \mathcal M }}^{2}\rho {{kT}}_{\mathrm{eq}}/\mu$ is negligible compared to the fluid's internal thermal energy. Consequently, Equation (15) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{tot}}}{\partial t}=(h+\phi )\displaystyle \frac{\partial \rho }{\partial t},\end{eqnarray} \tag{ 16 }$$

where E_tot = ρ e + ρ ϕ incorporates the internal thermal energy density and gravitational potential energy density. The first and second terms in parentheses on the right-hand side of Equation (16) correspond to the internal energy advection and change of gravitational energy. We integrate Equation (16) over the volume of the envelope and obtain the change of the total energy ΔE over a time interval Δt longer than the convection timescale:

$$\begin{eqnarray}&&{\rm{\Delta }}E=(h+\phi ){\rm{\Delta }}M,\end{eqnarray} \tag{ 17 }$$

where ΔE > 0 is the total energy gained by the system in order to expand, ΔM < 0 is the mass loss, and −(h + ϕ) is the minimal amount of energy needed per mass to lift material to infinity with temperature at infinity ${T}_{\inf }=0$. Recall that h + ϕ is approximately constant with radius and so can be evaluated at any point in the convective envelope.

To account for the outflow (Equation (4)) in the radiative atmosphere atop the convective zone, we must add several terms to Equation (17). The H/He mass taken from the envelope should be at least heated and inflated to T_eq and only needs to be lifted to the sonic radius, beyond which it is considered to be completely lost from the planet. The kinetic energy at the sonic radius is also considered here for a more realistic energy budget. This gives the total energy needed per second for the entire system to power the hydrodynamic outflow:

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{loss},\mathrm{tot}}=({c}_{p}{T}_{\mathrm{eq}}+\displaystyle \frac{1}{2}{c}_{s}^{2}-\displaystyle \frac{{{GM}}_{p}}{{R}_{s}}-h-\phi )\dot{M}=\left(\displaystyle \frac{{{GM}}_{p}}{{R}_{\mathrm{rcb}}}-\displaystyle \frac{3{{kT}}_{\mathrm{eq}}}{2\mu }\right)\dot{M},\end{eqnarray} \tag{ 18 }$$

where we have evaluated $h={c}_{p}{T}_{\mathrm{eq}}=(7/2){c}_{s}^{2}$ (with γ = 7/5 for molecular gas) and ϕ = −GM_p /R_rcb at the RCB and plugged in Equation (5). This amount of energy equals the energy needed by the non-energy-limited isothermal Parker wind in the radiative atmosphere:

$$\begin{eqnarray}\begin{array}{rcl}{\dot{E}}_{\mathrm{loss},\mathrm{rad}} & = & \left(\displaystyle \frac{1}{2}{c}_{s}^{2}+\displaystyle \frac{{{GM}}_{p}}{{R}_{\mathrm{rcb}}}-\displaystyle \frac{{{GM}}_{p}}{{R}_{s}}\right)\dot{M}\\ & = & {\dot{E}}_{\mathrm{loss},\mathrm{tot}}\approx \displaystyle \frac{{{GM}}_{p}\dot{M}}{{R}_{\mathrm{rcb}}}.\end{array}\end{eqnarray} \tag{ 19 }$$

For the last equality, as an approximation, we have ignored the kinetic energy $\tfrac{1}{2}{c}_{s}^{2}$ and the gravitational potential energy at the sonic point GM_p /R_s , as they are typically small (∼kT_eq/μ) compared to the gravitational potential energy at the RCB (GM_p /R_rcb) for planets with our initial RCB radius, which is 10%–20% of the sonic radius R_s (Equation (5)).

The equality of Equations (18) and (19) indicates that the energy for the whole system is due to the outflow above the RCB instead of the envelope, with the change of gravitational energy when mass is transported in the envelope completely from the advection of internal energy. Since the total energy change per mass lost h + ϕ during the re-equilibration can be approximated by the specific total energy e + ϕ at the RCB, the total energy per envelope mass is nearly conservative during boil-off. Therefore, we argue that the H/He mass can be treated as being taken from the envelope surface between each mass-loss time step. The outflow is not energy limited by the envelope as opposed to the GSS18 assumption in Equation (12), because the mass redistribution process in the envelope happens on a short convection timescale and releases the exact amount of internal energy required for envelope expansion as the required energy source, leading to a negligible change in the envelope energy budget (specific entropy). Note that the discussion above does not depend on the adiabatic index γ and where the envelope mass is concentrated.

To verify this behavior numerically, we set up hydrostatic and adiabatic (with constant adiabatic index) envelopes with different envelope masses but the same outer boundary conditions (as we find that the RCB pressure and temperature barely change during the boil-off) assessed without the self-gravity. The exact amounts from the left-hand side (the total energy difference) and the right-hand side (the analytical estimate we argued) of Equation (17) are found to be equal, as shown in Appendix A.

We display the results of a similar calculation based on our time-evolving numerical model in Figure 4 (solid gray and black lines). The analytical expression that estimates the energy needed per unit time ${\dot{E}}_{\mathrm{loss},\mathrm{rad}}$ (as an approximation we used the second equality in Equation (19); gray solid line) for the wind to transport mass in the radiative atmosphere matches the total energy gained in order to lose mass (black solid line).

Zoom In Zoom Out Reset image size

We compute the black curve in Figure 4 as follows. Between each numerical time step, we calculate the difference of the total energy ΔE of the envelope induced by the mass loss alone excluding any consequence from the thermal contraction and the mass transfer into the radiative region:

$$\begin{eqnarray}&&{\rm{\Delta }}E={\rm{\Delta }}U+{\rm{\Delta }}{E}_{{\rm{th}}}+L{\rm{\Delta }}t+{\rm{\Delta }}{E}_{{\rm{rcb}}},\end{eqnarray} \tag{ 20 }$$

where ΔU is the change of the gravitational binding energy, ΔE_th is the change of thermal energy (internal energy), LΔt is the energy lost due to thermal cooling represented by the total luminosity L out of the envelope (see Equation (1)) over a time interval Δt, and ΔE_rcb accounts for the total energy lost by the mass transfer into the radiative zone at the RCB. ΔE_rcb is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{rcb}}}=\left(-\displaystyle \frac{{{\rm{GM}}}_{p}}{{R}_{{\rm{rcb}}}}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{{{kT}}_{{\rm{eq}}}}{\mu }\right){\rm{\Delta }}{M}_{{\rm{atm}}},\end{eqnarray} \tag{ 21 }$$

where ΔM_atm is the mass transfer into the radiative atmosphere (note that Equation (21) is not the same as Equation (17), but rather is simply the gravitational and internal energy of the material removed from the convective zone and added to the radiative zone). The results are shown in Figure 4. The close overlap between the gray (analytical) and black (numerical) lines shows that the escaping mass can be treated as entirely lost from the RCB without any other energetic consequence for the adiabatic envelope.

From Equation (17), we derive the rate of change in total energy of the envelope as $\dot{E}=\phi \dot{M}+h\dot{M}$. The first term represents the change in gravitational binding energy, and the second term corresponds to the internal energy that is advected out of the envelope per unit time through the bulk motion, the same as in Equation (2). However, we argue that it does not decrease the specific entropy of the envelope as an envelope cooling, in opposition to the argument in J. E. Owen & Y. Wu (2016). Instead, it is a consequence of mass loss that removes the energy from the internal thermal energy reservoir by decreasing the system mass. The advective energy flux ${L}_{\mathrm{ad}}=h\dot{M}$ then becomes a constant throughout the assumed isothermal atmosphere, leading to no internal energy exchange in the radiative atmosphere, so it is inefficient in powering the outflow under our assumption (see discussion in Section 6.4).

3.1.3. Expansion of Radiative Atmosphere: The Bottleneck in the Deep Atmosphere

We next focus on the energy supply that sets the wind temperature and powers the outflow in the radiative atmosphere. In Figure 5, there are two possible energy sources: the stellar heating L_bol from above, and the planet's own cooling energy L_int from below. The stellar heating is assessed at radius R_h , where stellar photons entering the atmosphere reach a critical optical depth τ_h , defined to be the depth where external bolometric heating from the star is locally comparable to the energy needed to overcome PdV work. Using this evaluation radius, ${L}_{\mathrm{bol}}=\pi {R}_{h}^{2}{F}_{\mathrm{bol}}$, where F_bol is the stellar bolometric flux incident on the planet. We note that the local flux at this radius is ${F}_{\mathrm{bol}}{e}^{-{\tau }_{\nu }}$, where τ_ν is the optical depth to visible photons, but we define L_bol at R_h to avoid overestimating the energy available to the wind from stellar photons.

Zoom In Zoom Out Reset image size

As a comparison, the total energy needed to overcome PdV work (Equation (19)) is shown with the dashed red line. The figure is based on the evolution of a 3.6 M_⊕ planet with the non-energy-limited isothermal Parker wind mass loss in Equation (6). The intrinsic luminosity is dwarfed by the stellar bolometric luminosity over the entire evolution, such that the intrinsic luminosity is a secondary energy source in the radiative atmosphere. This finding is independent of the stellar flux we chose. We find that boil-off can be energy limited owing to the deficiency of the stellar bolometric energy deposit in the early stage (the transition is marked by circles). Based on Equation (19), the mass-loss rate is given by

$$\begin{eqnarray}&&{\dot{M}}_{e-\mathrm{lim},\mathrm{bol}}=\displaystyle \frac{{L}_{\mathrm{bol}}}{{{GM}}_{p}/{R}_{\mathrm{rcb}}},\end{eqnarray} \tag{ 22 }$$

which we term "bolometric-limited." As a planet shrinks as a result of the mass loss, the energy demand declines exponentially and the stellar bolometric energy that is directly deposited in the radiative atmosphere becomes sufficient to power escape of the atmosphere after ∼1 kyr.

However, a bottleneck region exists in the deeper part of the radiative atmosphere (Figure 2), where the stellar bolometric flux becomes exponentially more diffuse, turning into less useful PdV work. We find that even with sufficient total stellar energy deposited in the atmosphere, the inadequacy of energy injection for the outflow can still happen in a layer that is dark to visible radiation, above the RCB and below a radius that is optically thick to visible photons. In Figure 6 we calculate the ratio between the bolometric absorption and energy needed for each layer of the radiative atmosphere at an age after the bolometric-limited phase for the same model planet:

$$\begin{eqnarray}&&\displaystyle \frac{{PdV}}{\mathrm{absorption}}=\displaystyle \frac{{{GM}}_{p}\dot{M}{dr}/{r}^{2}}{\pi {r}^{2}\sigma {T}_{{\rm{eq}}}^{4}{e}^{-{\tau }_{\nu }}d{\tau }_{\nu }}=\displaystyle \frac{4{{vge}}^{{\tau }_{\nu }}}{\sigma {T}_{{\rm{eq}}}^{4}{\kappa }_{\nu }},\end{eqnarray} \tag{ 23 }$$

where g = GM_p /r² is the local gravity, v is the wind velocity, and κ_ν is the opacity to visible photons. To compute the second equality in Equation (23), we plugged in Equation (6), and τ_ν is defined as

$$\begin{eqnarray}&&{\tau }_{\nu }=\int \rho {\kappa }_{\nu }{dr}.\end{eqnarray} \tag{ 24 }$$

Above the critical optical depth to visible photons τ_h , each layer of the radiative atmosphere gets more energy from the direct absorption of the bolometric radiation than the amount needed to expand, whereas the deep atmosphere is inefficient in absorbing stellar energy. We find that the critical optical depth τ_h is weakly dependent on planetary radius (and therefore age), core mass, and stellar bolometric flux. τ_h = 10 typically gives a good estimate (see Figure 6 for an example). Note that the reemission process of the radiative atmosphere at infrared wavelengths is ignored in this calculation, which might assist in the energy absorption in the deep atmosphere and push the bottleneck to a deeper location.

Zoom In Zoom Out Reset image size

To break through the bottleneck, the only possible energy source is from the radiative interior cooling. In Figure 5, we show that the radiative intrinsic luminosity L_int (solid black) starts to be ample enough to expand the bottleneck region after 10⁴ yr, before which the mass loss should be energy limited owing to the bottleneck effect. For ease of comparison with other regimes, we name this process "intrinsic-limited" escape. As a comparison, the energy needed for the bottleneck region ${\dot{E}}_{\mathrm{loss},{\tau }_{h}}$ below τ_h is shown in solid red. The intrinsic-limited regime typically lasts longer and is more energy limited (with a lower mass-loss rate in Figure 7) than the bolometric-limited regime as long as the τ_h radius occurs above the RCB (circle). This makes it the common limiting effect for boil-off (in the example in Figure 7, the bolometric-limited regime is not relevant). We only find the bolometric-limited regime to be relevant for very-low-mass (≤1 M_⊕) planets at young ages, whose τ_h occurs physically below the RCB.

Zoom In Zoom Out Reset image size

Analytically, ${\dot{E}}_{{\rm{loss}},{\tau }_{h}}$ can be estimated by assuming a steady-state outflow:

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{loss},{\tau }_{h}}=\displaystyle \frac{{{GM}}_{p}\dot{M}}{{R}_{{rcb}}}-\displaystyle \frac{{{GM}}_{p}\dot{M}}{{R}_{{\tau }_{h}}},\end{eqnarray} \tag{ 25 }$$

where ${R}_{{\tau }_{h}}$ is the radius of the τ_h surface. The validity of this assumption can be verified by our numerical model in Figure 4 (dotted lines). Similar to the procedure for the convective envelope, the change of total energy between time steps is

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\tau }_{h}}={\rm{\Delta }}{U}^{{\prime} }-{\rm{\Delta }}{E}_{{rcb}}^{{\prime} }+{\rm{\Delta }}{E}_{{\tau }_{h}}^{{\prime} },\end{eqnarray} \tag{ 26 }$$

where ${\rm{\Delta }}{E}_{\mathrm{rcb}}^{{\prime} }$ and ${\rm{\Delta }}{E}_{{\tau }_{h}}^{{\prime} }$ are the gravitational binding energy fluxes carried by the bulk flux through the RCB and τ_h surfaces, respectively. (The primes indicate that these quantities are for the deep radiative region, while in Equation (20) they were for the convective envelope.) They consist of the mass transfer through the corresponding surfaces due to both the mass redistribution and the mass flux driven by the hydrodynamic wind. The actual energy need (black dotted line in Figure 4) is very slightly higher compared to the steady-state approximation (Equation (25) and dotted gray line) since the gravitational contraction of the bottleneck layer produces extra thermal energy to fuel the wind.

Therefore, in the presence of inadequate intrinsic luminosity deposited in the deep radiative atmosphere, the mass-loss rate of such an intrinsic-limited boil-off is estimated by

$$\begin{eqnarray}&&{\dot{M}}_{e-{\rm{lim}},{\tau }_{h}}=\displaystyle \frac{{L}_{{\rm{int}}}}{{{\rm{GM}}}_{p}/{R}_{{\rm{rcb}}}-{{\rm{GM}}}_{p}/{R}_{{\tau }_{h}}},\end{eqnarray} \tag{ 27 }$$

according to the energy needed by the bottleneck region in Equation (25). We suggest that a more comprehensive model could employ a mass-loss rate

$$\begin{eqnarray}&&\dot{M}=\min ({\dot{M}}_{{\rm{parker}}},{\dot{M}}_{e-{\rm{lim}},{\rm{bol}}},{\dot{M}}_{e-{\rm{lim}},{\tau }_{h}})\end{eqnarray} \tag{ 28 }$$

that incorporates both mass-loss regimes we discussed.

We have demonstrated the energetics of the radiative atmosphere and H/He envelope. We think that bolometric-driven escape might be a better name for the long-lived mass-loss phase, rather than core-powered mass loss, as the outflow is not powered by the core energy, nor is the envelope re-equilibration limited by the interior cooling energy. The thermal energy from the interior is still an important energy source owing to the bottleneck effect, but only in the deep part of the radiative atmosphere and at a young age. We emphasize that this conclusion does not imply that core luminosity is irrelevant—when L_core efficiently heats up the envelope, the planet remains inflated longer, which can contribute to increased mass loss. We merely conclude that mass loss is not limited by core luminosity.

3.1.4. Consequence of Different Wind Assumptions

3.1.5. Role of Core Luminosity

Though we have shown that the core luminosity does not limit the mass-loss rate by restricting the amount of energy available to supply the H/He mass to the outflow at the RCB, we agree with previous work (S. Ginzburg et al. 2016; W. Misener & H. E. Schlichting 2021; J. G. Rogers et al. 2024) demonstrating that core luminosity plays a role as an energy reservoir to delay the thermal evolution and enhance the mass-loss rate by slowing the thermal contraction of the planet (Figure 2). This effect especially becomes prominent in the later long-lived mass-loss stage when thermal evolution is coupled with the mass-loss process. We note that this effect is not unique to this particular mass-loss mechanism. Any other hydrodynamic wind, like XUV-driven escape, can also be enhanced by the large radius owing to the existence of the core (E. D. Lopez et al. 2012), given that the core makes up a large fraction of a planet's energy reservoir. As such, together with our previous discussion on the envelope re-equilibration in Sections 3.1.1 and 3.1.2, we prefer to not call the long-lived mass loss "core-powered mass loss," but instead we use the term "core-enhanced effect" to evaluate the role of core luminosity in amplifying mass loss.

We now reevaluate the role of core luminosity for our boil-off, which we argue is the mass-loss phase that is decoupled from thermal evolution, as discussed in Section 2.3.2. W. Misener & H. E. Schlichting (2021) argue that the gravitational binding energy of the envelope accounts for the majority of the envelope cooling of boil-off, with a negligible contribution from the core. As a comparison, in the bottom-right panel of Figure 8 we find that it is the core luminosity (dotted red line) that constitutes the main component of the envelope cooling (dashed and solid lines) during a vigorous boil-off phase. This is because the increased temperature contrast between the core and the base of the envelope enhances the flux of energy out of the core, when the temperature at the base of the envelope declines (note that this is not a result of thermal cooling but rather due to mass loss, which results in re-equilibration to an adiabat with a lower base temperature as described in Section 3.1.2). Indeed, the core luminosity can be comparable to the advective envelope cooling (solid line) before 10⁵ yr, inhibiting any envelope thermal contraction (solid black line in the bottom-left panel) until the radiative cooling (dashed line) starts to dominate the envelope energy budget close to 1 Myr, after which boil-off ceases. This makes the core luminosity a major interior energy source in boil-off (bottom-right panel).

This core luminosity, however, only has a limited role of enhancing the boil-off mass loss owing to the decoupling in timescale (bottom-left panel) between the mass loss (dashed line) and thermal evolution (solid line). Boil-off is essentially not core enhanced before a transition phase (dotted black line), which is defined as when the core luminosity potentially warms up the envelope faster than the mass-loss timescale ${t}_{\dot{M}}$ (Equation (7)):

$$\begin{eqnarray}&&{t}_{\mathrm{core}}=\displaystyle \frac{{{GM}}_{c}{M}_{\mathrm{env}}}{\alpha {R}_{\mathrm{rcb}}{L}_{\mathrm{core}}}\leqslant {t}_{\dot{M}}.\end{eqnarray} \tag{ 29 }$$

Once a planet evolves into this transition phase that happens close to the end of boil-off at ∼1 Myr, the core luminosity becomes effective in enhancing the mass loss. As the thermal evolution and mass-loss timescales would be comparable if no core luminosity were present, the presence of the core luminosity and hence the transition phase prolongs boil-off by delaying the thermal evolution, resulting in potentially greater mass loss. However, we find that this transition phase is typically shorter-lived than the early boil-off, with a substantially lower mass-loss rate, leading to a minor time-integrated mass-loss enhancement, as shown in Figure 8, in which we display the planetary evolution with (solid black line) and without the core luminosity (solid gray line). The decrease in the final mass fraction due to the core luminosity is typically less than 30%–40%.

Zoom In Zoom Out Reset image size

After the transition, a hydrodynamic mass-loss process is coupled with thermal evolution when ${t}_{\dot{M}}\gtrsim {t}_{\mathrm{cool}}$. However, in this stage the long-lived Parker mass loss has an exponentially increasing mass-loss timescale. Although the core luminosity can still constitute a large fraction of the planet's energy budget and enhances the mass-loss rate by a significant fraction, we find that it does not impact the absolute time-integrated mass loss of the later evolution owing to the orders-of-magnitude-slower mass-loss rate compared to the thermal contraction rate. Consequently, Parker mass loss decouples (meaning that the mass-loss rate becomes negligible, rather than being insensitive to core luminosity) from thermal evolution again (with ${t}_{\dot{M}}\gg {t}_{\mathrm{cool}}$) after a brief coupling phase lasting only a few megayears post-transition. Therefore, from the timescale aspect we did not find an efficient long-lived core-enhanced wind, as opposed to GSS18.

To improve the numerical stability for our evolution model, as a simplification, we do not allow a planet to thermally inflate over time, which forces the core luminosity L_core to be no greater than the envelope cooling luminosity. It is important to note that in this work "thermal inflation" specifically refers to an increase in envelope specific entropy. Thermal inflation does not necessarily result in radius inflation, as mass loss reduces the envelope mass and consequently decreases the planetary radius. Similarly, J. G. Rogers et al. (2024) did not find that the core luminosity thermally inflates their model planets with their parameterized core luminosity calculation. We find that with the simplification relaxed the core luminosity can only thermally inflate planets for a brief period of time, lowering the final mass fraction (leading to an overestimation of the final mass fraction in our model) by at most 50% rather than changing it exponentially. Therefore, we argue that our simplification does not impact the general behavior we investigate here. The validity of our core treatment and future improvements are discussed in Section 6.2.

A large initial H/He mass means a large radius planet, leading to a vigorous mass loss and more advective cooling (J. E. Owen & Y. Wu 2016). With our initial entropy defined in Section 2.4, although we find that the advective cooling rate is orders of magnitude higher than the radiative cooling rate, it is still negligible compared to the large energy reservoir. Due to the decoupling effect discussed in Section 3.1.4, our model is not sensitive to the assumed advective cooling when choosing different initial mass fractions, as these models will converge to the same final mass fraction. In this case, the strength of boil-off is only controlled by the initial entropy.

However, in contrast to our finding, J. E. Owen & Y. Wu (2016) reported a more significant advective cooling effect that leads to a cooler interior at the end of boil-off and therefore more substantial final mass fraction when choosing a higher initial mass fraction. This behavior must correspond to a thermal contraction timescale comparable to the mass-loss timescale. To explain the different results, we attribute the difference to initial conditions. Overall, the physical parameter that controls the effectiveness of the advective cooling is the initial RCB radius, which largely enhances the advective cooling rate and hence the post-boil-off mass fractions if the initial RCB radius is a large fraction of the sonic point radius. This corresponds to models with a hot interior state (with entropy >13k_b /atom for a 10% H/He envelope). However, it is unlikely to happen for a realistic planet, as it would correspond to a contraction rate, due to either the mass loss or advective cooling, much shorter than the disk dispersal timescale.

For comparison, starting with what we believe are more appropriate initial conditions, a typical low-mass planet has an initial RCB radius of only about 10%−20% of its sonic point radius, and the influence of the advective cooling is then limited. We show that the ratio between the timescales of the advective cooling and mass loss ${t}_{\mathrm{cool},\mathrm{ad}}/{t}_{\dot{M}}$ represents the relative importance of the advective cooling, or how strongly the advective cooling is coupled with the mass loss, which scales with a planet's initial RCB radius R_rcb,0 over the sonic point radius R_s :

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{{\rm{cool}},{\rm{ad}}}}{{t}_{\dot{M}}}=\displaystyle \frac{{{\rm{GM}}}_{p}{M}_{{\rm{env}}}}{\alpha {R}_{{\rm{rcb}},0}{L}_{{\rm{adv}}}}\displaystyle \frac{\dot{M}}{{M}_{{\rm{env}}}}=\displaystyle \frac{2(\gamma -1)}{\alpha \gamma }\displaystyle \frac{{R}_{s}}{{R}_{{\rm{rcb}},0}}.\end{eqnarray} \tag{ 30 }$$

For the planets with our initial radius, the timescale ratio is typically around 10, indicating a decoupling, whereas J. E. Owen & Y. Wu (2016) initial conditions with R_rcb,0 = R_s correspond to a strongly coupled boil-off mass loss. Additionally, we suggest that the ignored core luminosity is important in the J. E. Owen & Y. Wu (2016) setup owing to the coupling effect (see Section 3.1.5). Including the core luminosity would greatly increase the interior cooling timescale (Equation (9)) and weaken the coupling between the thermal evolution and mass loss, leading to a less significant advective cooling effect. The quantitative details and a direct comparison to J. E. Owen & Y. Wu (2016) are included in Appendix B.

Due to this sensitive behavior of boil-off, we always initialize boil-off using our self-consistent initial entropy in the later discussion unless otherwise specified. An example of the evolution track is shown in Figure 9, where the greater role of the advective cooling L_ad for the 30% planet is seen at 10² yr, which rapidly cools the planet, slowing mass loss and leading to a larger post-boil-off mass fraction compared to the 10% model. These enhancements are considered to be small. Note that the initial specific entropy of the 30% model is set to be the same as the self-consistent entropy of the 10% model to eliminate the initial entropy effect. Moreover, this effect is largely reduced at a lower initial mass fraction (≤20%) that is predicted by accretion models. In the presence of core heating, the net cooling effect is further diminished. These yield a relative final mass fraction difference smaller than 10% between the models with and without the advective cooling. Therefore, we can safely ignore the advective cooling when choosing a different initial mass fraction ≤20% with our determined initial entropy, as the final H/He mass fraction is insensitive to the initial mass fraction.

Zoom In Zoom Out Reset image size

Our finding above with the assumption that those planets start with the same initial entropy (as shown in Figure 9) still holds true with a more sophisticated initial condition. The specific entropy required to start with a Kelvin–Helmholz cooling timescale of 5 Myr is slightly higher for envelopes with greater initial mass fractions. For planets with reasonable initial mass fractions in the range of 10%−20%, we find that the variation in the self-consistent initial entropy results in only a marginal change of only 10% in the post-boil-off envelope mass compared to the results using the fixed initial entropy. Therefore, this effect can be neglected.

Moreover, we find that the effectiveness of advective cooling depends on not only the timescale ratio ${t}_{\mathrm{cool},\mathrm{ad}}/{t}_{\dot{M}}$ but also the ratio between the cooling timescale t_cool and the current planetary age. If the cooling timescale is longer than the current age, the entropy is unable to efficiently decline through cooling until the planet becomes older. For instance, in the bottom-right panel of Figure 9 we show the evolution of timescales as a function of time. Even though the advective cooling is always 10 times slower than mass loss for the 30% model (solid black) before 1 Myr, the entropy can still significantly drop when the cooling timescale becomes comparable to a certain fraction of the current age between 1 and 100 yr. We find that a 1/10 ratio between the timescales and age (dotted red) gives a good match for the behavior. On the contrary, for the 10% model the advective cooling timescale is always much longer than the age to drive significant cooling, which is the reason that the ≤20% models are not susceptible to advective cooling. Similarly, around 1 Myr, when the cooling timescale becomes shorter than the age, radiative cooling kicks in and thermal contraction becomes considerable again, after which boil-off shuts off. This behavior also applies to mass loss. In the top-right panel, the H/He mass fraction for the 30% model (black) barely declines with time until the mass-loss timescale (dashed black line in the bottom-right panel) becomes comparable to the current age after 1 yr.

To summarize, our treatment of the advective cooling effect is based on the treatment (Equation (2)) from J. E. Owen & Y. Wu (2016). We find that it is negligible for boil-off with our self-consistent initial conditions. Additionally, under our model assumptions and following our theoretical discussion (see Section 3.1.2), we do not find that the energy advection ${L}_{\mathrm{ad}}=h\dot{M}$ has a cooling effect, as the energy advection is a consequence of mass loss from the envelope, which does not affect the thermal evolution. Therefore, our results in the later sections are independent of the treatment of advective cooling.

Here we set up our numerical model with self-consistent initial conditions to assess the conditions a planet will have after boil-off. As we showed that the post-boil-off mass fraction is independent of initial mass fractions ≤20%, as a simplification we start our boil-off models with a constant initial H/He mass fraction of 10% to constrain the final mass fraction after boil-off, unless specified. With our initial entropy choice (see Section 2.4), the final mass fraction after boil-off depends solely on the scale height of the radiative atmosphere, which is a function of core mass and incident stellar flux.

We shut off boil-off either at the time t_end that sets the transition into the long-lived mass loss in Equation (8) or at the time t_XUV when the XUV heating starts to dominate the upper radiative atmosphere (Section 2.3.3). Practically, we find that t_XUV happens slightly earlier than t_end, both around 1 Myr under self-consistent initial conditions, with the time-integrated mass loss negligible between t_XUV and t_end, as the mass-loss rate has greatly diminished with time toward the end of boil-off. Therefore, we suggest that photoevaporation should start right after boil-off around 1 Myr, when the isothermal Parker wind starts to transition into the long-lived mass-loss phase. We find that the post-boil-off Parker wind is briefly as efficient as photoevaporation right after boil-off for a few megayears and then exponentially decays, making it a secondary mass-loss mechanism over long timescales. The role of the long-lived—what we call bolometric-driven—escape in the absence of the photoevaporation is discussed in Section 3.5.

With the initial conditions we chose, we find that the duration of the boil-off phase is typically in the range of 0.5–1 Myr. Right after boil-off, the planet has a contraction timescale of 10–20 Myr owing to the radiative cooling that happens around the end of boil-off, compared to the initial contraction timescale of 5 Myr shown in Table 1. For comparison, J. E. Owen & Y. Wu (2016) find a longer 50 Myr timescale starting with a substantially larger radius. We attribute this difference to a different criterion for terminating boil-off. If we let our boil-off model evolve for another 2 Myr as they do, we find results consistent with J. E. Owen & Y. Wu (2016). Therefore, the cooler appearance of a planet after boil-off is independent of the initial conditions and the types of cooling mechanisms. In fact, we find in our model that it is the radiative cooling that prolongs the contraction timescale in our model, compared to the advective cooling as argued in J. E. Owen & Y. Wu (2016). We note that in our model most cooling happens near the end of boil-off, when the radiative thermal contraction starts to be efficient, shown in Figure 9.

Table 2 shows a grid of the final mass fraction that a planet can have after boil-off and before photoevaporation or the post-boil-off bolometric-driven escape, as a function of core mass and bolometric flux. Low-mass planets and highly irradiated planets tend to lose more mass owing to their large scale height of the radiative atmosphere. Note that most of the models start with an envelope mass fraction of 10% except for heavy cores and cold planets. Since they are less susceptible to boil-off, which allows them to hold a final mass fraction greater than 10% if born with a sufficient amount of H/He mass, we assess them using a higher initial mass fraction. The final mass fraction also represents a critical mass fraction, in that we find that planets with higher initial fractions than the critical mass fractions shown in the table always converge, "boiling off" until such final mass fractions are reached. On the other hand, planets formed with initial mass fractions lower than the critical mass fractions will not undergo boil-off. Heavy cores and less irradiated planets are also invulnerable to boil-off, potentially holding any substantial final mass fraction. In this case, photoevaporation and/or post-boil-off bolometric-driven escape starts immediately after the disk dispersal.

Remarkably, we find that boil-off is so powerful that it is capable of removing the entire convective envelope, leaving behind either a bare rock/iron core or a super-Earth with only a layer of radiative atmosphere. This behavior typically happens to planets with core mass <4 M_⊕ and stellar insolation >100 F_⊕, due to the nature of the large scale height, resulting in a less bound and unstable atmosphere. It indicates that boil-off can possibly contribute significantly to the transformation of sub-Neptunes into super-Earths and can potentially have an impact on the formation of the radius valley. On the other hand, beyond 1000 F_⊕ and 20 M_⊕, our model is consistent with the existence of warm Neptunes, such as GJ 436b, which are planets that survived the boil-off phase with heavy envelopes and are rapidly losing mass through photoevaporation.

Our caveat for our treatment and Table 2 is that our boil-off model implicitly assumes the collisional regime for hydrodynamic winds. When the atmosphere is tenuous enough, the gas starts to become collisionless, known as Jeans escape. The lowest mass fractions (lower left corner) are so small that they would be subject to a Jeans escape, which would greatly slow down the mass-loss rate. Therefore, the quantitative values for them may not be accurate. Despite that, low-mass and highly irradiated planets can never hold a convective envelope based on our model and will be eventually turned into super-Earths in a short time in the presence of photoevaporation and/or post-boil-off bolometric-driven escape, once the envelope is thin enough.

To assess the impact on planet population from boil-off, we evolve planets in the absence of subsequent mass loss by continuing thermal evolution for another 10 Gyr with the initial conditions from the end of boil-off (Table 2). This constrains the largest possible radius allowed as a function of stellar bolometric flux and core mass at each given age. Intuitively, low-mass sub-Neptunes tend to possess a large radius and incredibly low bulk density even with a small amount of H/He gas, due to their large scale heights of the radiative atmosphere compared to their physical sizes. For these planets, the boil-off phase results in very thin envelopes or, in most cases, a complete loss of envelopes. On the other hand, heavier planets that have experienced a boil-off hold more massive envelopes, but their greater surface gravity allows them to have a high bulk density and remain relatively compact. Therefore, boil-off plays a role in setting an upper boundary in the radius distribution of the sub-Neptune population, which may contribute to carving the radius cliff, the decrease in planets with radii above ∼4 R_⊕ from observation.

We show the radius evolution in Table 3. Most of the planets shrink to radii smaller than 5 R_⊕ by 10 Gyr, roughly consistent with the radius cliff but with a larger radius. As a comparison, without a boil-off phase, the radii of low-mass and highly irradiated planets, which dominate the sub-Neptune population, are unbounded, potentially exceeding 10 R_⊕ even at a late age. Note that the massive planets that are invulnerable to boil-off are rare in the observed population and may eventually evolve into the warm Neptune population. Our result implies that, rather than photoevaporation or core-powered escape (the post-boil-off Parker wind), boil-off is likely the main cause of the radius cliff, which significantly removes the H/He envelope from the planets with large physical sizes. However, other mass-loss mechanisms are needed to further strip planets' envelopes and shrink their radii. The possibility that the post-boil-off bolometric-driven escape sculpts the radius cliff is discussed in the next section.

After the boil-off ends, planetary thermal evolution becomes more significant than mass loss in controlling the planetary radius. Consequently, the scale height of the radiative atmosphere rapidly diminishes, leading to an exponentially decaying mass-loss rate (J. E. Owen & Y. Wu 2016; see also our Figure 7, which shows a low mass-loss rate after 10⁶ yr), which almost evolves into a hydrostatic equilibrium. However, in principle, a planet can still lose a substantial amount of mass on a longer timescale. Here, to evaluate the importance of the post-boil-off Parker wind as a long-term mass-loss mechanism, we continue to evolve the planets after boil-off, but with the Parker wind mass loss turned on for another 3 Gyr. For all of the modeled planets in Table 2, we find that the time-integrated mass loss in the presence of the core luminosity out to a few gigayears due to the post-boil-off Parker wind is a negligible amount (≤0.1%) compared to that from boil-off (∼10%; see also Figure 1 for a similar evolutionary track), leading to a nearly constant mass fraction at late ages as demonstrated in Figures 10 and 11.

Zoom In Zoom Out Reset image size

In these figures, we assess the planet's initial mass fraction for the post-boil-off mass loss based on the final mass fraction at the end of boil-off (triangle) to separate the outcomes from the two mass-loss mechanisms. This treatment is essential, as different initial mass fractions lead to divergent H/He loss histories after a few gigayears of evolution owing to the mass-loss process being coupled to significant thermal evolution, which greatly depends on the conditions it starts with at the end of boil-off. As a comparison, in GSS18 the post-boil-off (their initial) mass fraction was derived analytically with particular assumptions, such as that the envelope radius R_rcb is twice the core radius R_c , and without directly modeling a boil-off phase. These authors found a scaling function for the H/He mass fraction, f, where $f=0.05{({M}_{c}/{M}_{\oplus })}^{1/2}$, with M_c the core mass. As a direct comparison to the work developed in GSS18, we initialize planets using their analytic H/He mass fraction scaling function (black curve) to examine the validity of their initial conditions to start the subsequent long-term evolution. This group of planets has the Parker wind turned on. They are compared to the models initialized with our 10% default initial mass fraction (gray curve), with the Parker wind turned off. The choice of the initial conditions does not affect the evolutionary results from our numerical model, as with self-consistent initial entropies those models from different initial mass fractions converge to common final mass fraction values at the end of boil-off as seen in Figure 10 (see also Section 3.2). The initial entropies are chosen to be the same for both models.

A single representative evolution curve for the H/He mass fraction between the two models is shown in Figure 10. We find a very tiny difference between the final mass fractions after 3 Gyr (crosses; based on the GSS18 initial mass fraction with the bolometric-driven escape turned on) and those from the end of boil-off around 1 Myr (triangles; based on our 10% boil-off models; same data as in Table 2). More examples with a wide range of core masses and incident stellar bolometric fluxes are shown in Figure 11, showing that there is very little Parker wind mass loss during this long time. Therefore, we argue that the post-boil-off Parker wind is negligible on a timescale of gigayears. This is because the long-lived mass loss does not have enough time to significantly remove envelope mass with ${t}_{\dot{M}}\gg t$, where t is the planetary age, due to the orders-of-magnitude difference between the mass-loss timescale and thermal contraction timescale (which is comparable to the age t) at most times of the post-boil-off evolution.

Zoom In Zoom Out Reset image size

In rare cases, especially for planets with low core masses and high equilibrium temperatures, our model finds that an envelope can only be totally lost through the post-boil-off Parker wind if the post-boil-off H/He mass fraction is small enough, ∼a few hundredths of a percent of H/He mass. This happens swiftly within a few megayears right after boil-off, rather than the few gigayears that is argued for the timescale for core-powered mass loss, if the physical choices are properly made.

Additionally, a self-consistent initial mass fraction should start planetary evolution right out of the boil-off phase, with thermal evolution comparable to or faster than mass loss. However, the GSS18 post-boil-off mass fraction scaling function has a rather flat profile as a function of core mass, in contrast to the steep post-boil-off mass fraction distribution from our model (and W. Misener & H. E. Schlichting 2021) that quickly declines to zero toward the low-mass end (Table 2). This underestimates the post-boil-off mass fraction on the high-mass end and overestimates it on the low-mass end. In Figure 11, the GSS18 initial mass fraction function is shown by a dotted line, while the self-consistent initial fractions for the long-term evolution that our model finds correspond to triangles, which indicates that a boil-off would happen, instead of a long-lived core-powered mass loss, for nearly all low-mass planets if they started with the GSS18 initial fractions and were assessed with what we believe are more appropriate decoupled Parker winds (Section 3.1.4) and self-consistent initial entropies. We let these planets continue to boil off until these planets enter the "thin regime." According to GSS18 and W. Misener & H. E. Schlichting (2021), this is defined as when the envelope radius shrinks to twice the core radius, R_rcb = 2R_c . We still find that some planets are boiling off with a short mass-loss timescale <1 Myr. This indicates an inconsistency with the GSS18 assumption that the planets have already exited the boil-off phase and entered the long-lived mass-loss phase.

The above initial condition problem is quantitatively identified based on our numerical model. With the GSS18 initial conditions (radius and envelope mass fraction), we calculate the ratio between the mass-loss timescale of the decoupled wind and the Kelvin–Helmholtz thermal contraction timescale at the beginning of their thin regime evolution, which we find is many orders of magnitude smaller than 1 for highly irradiated planets and low-mass planets, shown in Table 4. From the previous discussion, a similar behavior is found at their transition time R_rcb = 2R_c , if the initial envelope mass fraction is instead self-consistently calculated from boil-off as in W. Misener & H. E. Schlichting (2021), indicating that R_rcb = 2R_c may not be a good criterion for the transition, as it can correspond to boil-off physical conditions.

In summary, in this section, since we show that the long-lived post-boil-off Parker wind is inefficient with a wide range of core mass and incident flux given that a boil-off phase is included, we suggest that atmospheric photoevaporation, whose importance in impacting the demographics has been widely studied (E. D. Lopez & J. J. Fortney 2013; J. E. Owen & Y. Wu 2013, 2017), should be the dominant mass-loss mechanism that takes place right after boil-off, which further sculpts the small-planet population and in part forms the observational features such as the radius valley and radius cliff. A modeling work shows a good fit of the combined evolution with both boil-off and post-boil-off photoevaporation to the observed radius gap (L. Affolter et al. 2023). Additionally, our numerical model indicates that core-powered escape is not a long-lived atmospheric escape but part of our decoupled Parker wind boil-off under the GSS18 initial conditions if modeled with the decoupled Parker wind, which yields initial planets with too large radii for low-mass planets. A direct analysis of the analytical model will be further discussed in Section 4.

We showed in Section 3.1 that the energy needed for adiabatic expansion is only the amount needed to expand the radiative atmosphere from the RCB (see Equation (19)), compared to the GSS18 assumption that energy is required to lift material from the deep envelope right above the core (Equation (11)). In Figure 12, we show an experiment based on the implementation for the analytical core-powered mass-loss model from GSS18, with different combinations of physical choices. Here we compare the GSS18 wind (red) to our decoupled (by our default non-energy-limited) isothermal Parker wind (dashed black), with the initial conditions the same as in the original setup. The delayed and prolonged boil-off seen with their analytical model is similar to that in Figure 10 with our numerical model. We show that by 30–40 Myr the planet with their coupled wind is still quickly boiling off, instead of starting the evolution within the long-lived mass-loss phase. This is induced by the overestimated initial H/He mass fractions and/or initial radius, as demonstrated in Section 3.5. Once the extended but less intense boil-off is over, when the wind transitions into the Bondi-limited regime and thermal evolution starts to dominate, the planet holds a larger quantity of H/He mass and thermally contracts faster than the mass loss, leading to an insignificant mass loss in the later evolution. As a comparison, the planet under the decoupled isothermal Parker wind holds a substantially less massive envelope at the end of boil-off <1 Myr, with the mass loss nearly completely quenched in the later gigayear. Planets under their coupled mass loss never converge to the final mass fraction value found by any of the decoupled Parker wind models, due to the coupling with the thermal evolution.

Zoom In Zoom Out Reset image size

As a planet loses H/He envelope mass, the gravitational energy and the internal thermal energy of the escaping matter should be correspondingly lost from the system, leading to a decline of the energy reservoir, which is an energy-loss term in addition to the thermal cooling that is not implemented in the original work. From GSS18, the energy reservoir is given by

$$\begin{eqnarray}&&{E}_{\mathrm{cool}}=g{\rm{\Delta }}R\left(\displaystyle \frac{\gamma }{2\gamma -1}{M}_{E}+\displaystyle \frac{1}{\gamma }\displaystyle \frac{\gamma -1}{{\gamma }_{c}-1}\displaystyle \frac{\mu }{{\mu }_{c}}{M}_{c}\right),\end{eqnarray} \tag{ 31 }$$

where ΔR = R_rcb − R_c is the thickness of the envelope, μ and μ_c are the mean molecular weights of the envelope and the core, respectively, and γ and γ_c are the adiabatic indices. The energy-loss rate as a result of mass loss can then be estimated by the derivative of the energy reservoir equation with respect to the envelope mass M_env:

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{adv}}={\dot{E}}_{\mathrm{cool}}=g{\rm{\Delta }}R\displaystyle \frac{\gamma }{2\gamma -1}{\dot{M}}_{E}.\end{eqnarray} \tag{ 32 }$$

${\dot{E}}_{\mathrm{adv}}$ represents an energy advection, similar to Equation (2), and ${L}_{\mathrm{adv}}=h\dot{M}$ in Section 3.1.2. The thermal evolution within a time step is therefore implemented in our version of the analytical model as

$$\begin{eqnarray}&&{E}_{{\rm{cool}}}\to {E}_{{\rm{cool}}}-L{\rm{\Delta }}t-{\dot{E}}_{{\rm{adv}}}{\rm{\Delta }}t.\end{eqnarray} \tag{ 33 }$$

In Figure 12, we demonstrate the importance of the energy-loss term ${\dot{E}}_{\mathrm{adv}}$ by comparing two models with (solid black line) and without it (dashed black line), both implemented with decoupled isothermal Parker winds, with all other conditions the same as the original GSS18 setup. It is seen in the dashed black line that without ${\dot{E}}_{\mathrm{adv}}$ a planet cannot efficiently contract and can even greatly inflate its radius in a short time interval in the boil-off period, because the total energy E_cool stays nearly constant in boil-off phase, due to both the inefficient thermal cooling and the lack of energy-loss term ${\dot{E}}_{\mathrm{adv}}$, and consequently a decreasing envelope mass M_env leads to an expanding envelope ΔR in Equation (31). This unrealistic radius inflation results in an exaggerated time-integrated mass loss. This behavior becomes especially pronounced, which generates a spike around 10⁴ yr as seen by the dashed black line, as we switch to a decoupled isothermal Parker wind from the GSS18 coupled mass loss (red). With the mass-loss-driven energy-loss term ${\dot{E}}_{\mathrm{adv}}$ implemented shown in black, we find a higher final mass fraction by 3 Gyr compared to the planet without it (dashed black line), with most of the envelope mass lost in 1 Myr, which is consistent with the behavior we find for boil-off based on our numerical model. The unrealistic inflation behavior (dashed black line) is also eliminated by ${\dot{E}}_{\mathrm{adv}}$, implying its significance. The black model with both physical improvements is what we argue is a more appropriate setup for modeling boil-off and the subsequent long-lived mass loss.

The inclusion of the ${\dot{E}}_{\mathrm{adv}}$ term leads to an enhanced final mass fraction, while a decoupled isothermal Parker wind lowers the final mass fraction. Interestingly, with the physics of both thermal evolution and mass loss improved (black), these two effects largely cancel each other out, producing very similar final mass fractions and final radii to those from the original setup without both improvements (red) if GSS18 initial conditions are used. Therefore, the impact from the improved physics is hard to evaluate on a population level (see Section 5) at gigayear ages, but from an evolutionary perspective, the net effect with both improvements leads to a mass-loss timescale shift from gigayear to megayear, indicating a dramatic mass-loss timescale difference between our decoupled Parker wind and the GSS18 coupled Parker wind if assessed without the ${\dot{E}}_{\mathrm{adv}}$ improvement. The timescale difference can be tested against the observed population in future studies.

With the improved physics for thermal contraction and mass loss applied to the analytical model, we compared the evolution between different models. Although the analytical model is consistent with our numerical model in terms of mass-loss timescale ∼1 Myr, we find that it generally predicts a higher amount of mass loss than our numerical model. As shown with the gray line model in Figure 12, this behavior results from the fact that a planet from the GSS18 formulation of the energy budget in Equation (31) does not shrink efficiently until its envelope mass fraction falls below 0.1% as the H/He mass fraction is gradually lost (see radius versus H/He mass fraction in the bottom-left panel). As a comparison, all model planets from our numerical model continuously shrink in radius during boil-off, as shown in Figure 1, due to the mass loss (rather than thermal contraction). This holds true even if we account for thermal inflation.

The constant RCB radius observed in the analytical model must be attributed to a continuous thermal inflation that increases planetary radius more efficiently than the mass loss decreases it. As a comparison, a planet with its mass dominated by the core losing mass at a fixed entropy is expected to undergo a decrease in radius. This inflation behavior corresponds to a tremendous heating from the core to the envelope, which exceeds the intrinsic luminosity with ${L}_{\mathrm{core}}\gg {L}_{\mathrm{int}}$, over the whole mass-loss phase. Note that the core luminosity and specific entropy are not explicitly calculated in the analytical model. The increased entropy therefore leads to more mass loss, due to the initial entropy effect discussed in Sections 2.4 and 3.2. This continuous and significant inflation is not found in our model if we relax our assumption and allow the planet to thermally inflate.

Additionally, the inflation in the analytical model also reduces the RCB pressure by a factor exceeding 1000 (bottom right panel of Figure 12). This suggests a dramatic increase in specific entropy in their model, since a shallow RCB corresponds to a hotter interior. In comparison, our numerical model suggests that the RCB pressure varies by at most a factor of 2 during boil-off with the thermal inflation included. This indicates that our numerical model predicts a much more modest thermal inflation and a significantly less pronounced role for core luminosity compared to the GSS18 model.

In this section, we reexamine the GSS18 coupled wind and, with improved physics, obtain results that diverge from those of GSS18. In W. Misener & H. E. Schlichting (2021) (Figure 2), these authors show that the GSS18-type energy-limited wind exhibits a strong coupling between the mass loss and thermal evolution, with comparable timescales ${t}_{\dot{M}}\sim {t}_{\mathrm{cool}}$. From Equations (9) and (12), we show that the wind treatment leads to a slowly (linearly) changing ratio between the mass-loss timescale ${t}_{\dot{M}}$ and radiative thermal contraction timescale t_cool in the boil-off phase:

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\dot{M}}}{{t}_{\mathrm{cool}}}=\displaystyle \frac{{M}_{\mathrm{env}}}{{L}_{\mathrm{int}}/{{gR}}_{c}}\displaystyle \frac{{L}_{\mathrm{int}}-{L}_{\mathrm{core}}}{{{GM}}_{c}{M}_{\mathrm{env}}/(\alpha {R}_{\mathrm{rcb}})}=\displaystyle \frac{\alpha (1-\beta ){R}_{\mathrm{rcb}}}{{R}_{c}},\end{eqnarray} \tag{ 34 }$$

where the luminosity ratio $\beta ={L}_{\mathrm{core}}/{L}_{\mathrm{int}}$ dictates the strength of the core-enhanced effect (this is true only when the thermal evolution is coupled with the mass loss with ${t}_{\dot{M}}/{t}_{\mathrm{cool}}\sim 1$).

In the top panel of Figure 13, we show the difference in the timescale ratio ${t}_{\dot{M}}/{t}_{\mathrm{cool}}$ between the linear behavior of the GSS18 coupled wind (red) and the exponential behavior (black) of the decoupled isothermal Parker winds. Note that for the decoupled Parker wind the ratio only becomes unity (coupled) briefly in the correspondingly core-enhanced transition phase (close to the end of boil-off; see Section 3.1.5). It is otherwise not core enhanced despite the greater β during most times of the boil-off, as a result of the decoupling between the mass loss and thermal evolution. As a comparison, the GSS18 coupled Parker wind can be much more core enhanced for a longer period of time, leading to a higher mass-loss enhancement. This is because a planet is imposed by the wind assumption to not have a vigorous boil-off phase but instead experience a mild and significantly extended transition phase according to Equations (7) and (29), which would be much shorter and more abrupt if assessed with what we believe is more appropriate decoupled Parker wind.

Zoom In Zoom Out Reset image size

Since the energy-limited regime happens early in the boil-off phase owing to the high energy demand, during which the intrinsic energy L_int is argued to be mostly from the gravitational binding energy of the envelope rather than the core thermal energy reservoir (S. Ginzburg et al. 2016; W. Misener & H. E. Schlichting 2021), β was considered negligible in the previous work in this mass-loss phase. In contrast, in our numerical model using the same coupled Parker wind assumption, we find that the core luminosity is always a significant fraction of the total cooling. For most planets β is typically no greater than 0.9 before transitioning into the non-energy-limited (what they call Bondi-limited) regime, except for those planets nearing complete atmospheric stripping. Consequently, we find that the timescale ratio is typically of order unity during the energy-limited phase. Once the coupled wind transitions into the non-energy-limited phase, the timescale ratio starts to increase exponentially, driving a limited amount of time-integrated mass loss after then, which we find with our numerical model is at most a few tenths of percent of H/He mass.

Compared to GSS18 (red in Figure 12) and W. Misener & H. E. Schlichting (2021), which show a significant number of completely stripped planets, our numerical model (red curve in Figure 10) finds a much lower susceptibility for these planets to the GSS18 coupled wind. A similar pattern emerges when using the GSS18 model with the inclusion of the energy-loss term ${\dot{E}}_{{\rm{adv}}}$, as shown in Figure 12, where we find that the GSS18 coupled wind removes the envelope mass significantly slower (gray), compared to the original approach without ${\dot{E}}_{{\rm{adv}}}$ (red). This underscores the importance of the ${\dot{E}}_{\mathrm{adv}}$ term in the analytical model. We find that a sub-Neptune is more frequently the final evolutionary product after a few gigayears. The divergent results from the original work are because their lack of ${\dot{E}}_{\mathrm{adv}}$ term magnifies the importance of the core luminosity by thermally inflating planets, leading to an exaggerated demographic feature.

The above behavior is expected from Equation (34) because the mass-loss process barely has enough time to entirely remove the envelope ${t}_{\dot{M}}\sim {t}_{\mathrm{cool}}\sim t$, where t is the planetary age. We find that a complete loss of the envelope rarely happens and only happens when the core luminosity is very comparable to or greater than the intrinsic luminosity with β ≳ 1 (the constant radius found with the analytical model as a planet loses mass in Section 4.3 corresponds to β > 1), which leads to a significant decrease of the timescale ratio. These stripped planets have very low core masses, <2 M_⊕, with low initial mass fractions <5% and high insolation >100 F_⊕. This suggests that the GSS18 coupled wind, when combined with improved thermal evolution physics, is not effective in forming a pronounced radius gap.

In Figure 14, we vary the mean molecular weight of the core μ_c in the analytical model, where μ_c defines the core heat capacity, from 34m_H (solid, which corresponds to a large core luminosity) to 340m_H (dashed, negligible core luminosity) based on the GSS18 initial conditions with (black) and without (gray) both the physical improvements, i.e., the decoupled Parker wind and the ${\dot{E}}_{\mathrm{adv}}$ term. The planetary evolution with the improvements is only marginally affected by the thermal behavior of the core, which we find is not sensitive to different combinations of initial conditions. Similar to the results from our numerical model in Section 3.1.5, the core-enhanced mass-loss effect only moderately impacts the transition phase around 10⁵ yr rather than the boil-off and the long-lived mass-loss phases before and after then, respectively (black curves). Consequently, this effect can only change the final mass fraction by a factor of order unity rather than altering the evolutionary fate of a planet. In contrast, without both improvements, the evolution with the GSS18 physics is very sensitive to core heat capacity (gray curves).

Zoom In Zoom Out Reset image size

In GSS16, these authors emphasized that a planet switches to the long-lived core-powered escape from the boil-off phase once it reaches the thin regime with R_rcb = 2R_c . After that point, the thermal contraction is argued to be greatly hindered by the core luminosity, which constitutes the majority of the envelope cooling, leading to a core-enhanced mass loss. Their criterion for the transition, however, is not self-consistent and is not consistent with our numerical results in Section 3.5. We find that the radius ratio at the end of boil-off, which is assumed to be a fixed number by their criterion, can vary greatly with different planetary parameters, i.e., core mass and stellar insolation, if assessed with what we believe is more appropriate decoupled Parker wind, tabulated in Table 5.

W. Misener & H. E. Schlichting (2021) used a different criterion, such that boil-off coincides with the energy-limited (coupled) regime whereas the "core-powered" mass loss exclusively refers to the Bondi-limited (decoupled) regime. Here we reevaluate these criteria with our numerical model. As their coupled wind assumption (see Section 4.4) essentially corresponds to a prolonged ∼1 Gyr and gentle transition phase, before and after both the GSS18 transition R_rcb = 2R_c (crosses in Figure 13) and the W. Misener & H. E. Schlichting (2021) transition (triangles), we find that the luminosity ratio β does not change significantly, shown in the bottom panel, and therefore has a similar impact on the mass-loss enhancement at both stages (see Section 3.1.5) shown in the middle panel. This is in contrast with their argument that the core luminosity is only important after their transitions. We suggest that a good criterion for transitioning into a significantly core-enhanced phase is β ∼ 1 (or with 1 − β ≪ 1 if a thermal inflation is not considered), which decreases the timescale ratio in Equation (34) by orders of magnitude, leading to a significant change of the mass loss. Since our findings show that β hardly exceeds 1 with improved physics, we do not consider a need to distinguish between the core-powered phase and boil-off phase for their coupled wind.

Consequently, both the GSS18 and W. Misener & H. E. Schlichting (2021) criteria for the transition give rise to a large uncertainty for the definition of their long-lived core-powered mass loss (shown after the crosses/triangles and before the circles on the red curve), which we argue should both correspond to our decoupled Parker wind boil-off if the mass-loss physics is improved, as shown in Section 3.5. As a comparison, our criterion for the transition in Equation (8) is independent of different numerical model setup and planetary parameters and reflects the characteristic mass-loss timescales. According to our transition criterion, the flat sections of the H/He evolution in Figures 12 and 13 account for the actual long-lived Parker wind phase.

Due to the coupling effect of the GSS18 wind, any initial conditions would lead to an initial timescale ratio ${t}_{\dot{M}}/{t}_{\mathrm{cool}}$ very close to 1 (Equation (34)). In this case, the importance of the initial conditions could be easily overlooked. Moreover, as discussed in Section 4.5, the uncertainty from the criteria for the transition used in the previous work further complicates the initial condition problem for the post-boil-off mass loss. Consequently, with these treatments, the role of boil-off and the long-lived escape are largely entangled. In Sections 4.6.1 and 4.6.2, we show the consequence of the GSS18 initial mass fractions and initial radii. The improvement of the initial conditions is shown in Section 4.6.3.

4.6.1. Initial RCB-to-core-radius Ratio Is Fixed

4.6.2. Initial H/He Mass Fraction Scaling Function Is Flat

Since the intensity of an isothermal Parker wind is controlled not only by the planet radius but also by the RCB density, the initial condition problem is complicated by the initial H/He mass fraction. In GSS16 and GSS18, these authors analytically derived an envelope mass fraction scaling function $f=0.05{M}_{c}^{1/2}$, which they found is roughly 30%–50% of the initial mass fractions from planet formation. As a comparison, in Section 3.5 (Figure 11), with our numerical model, we demonstrated a dramatic variation in the final mass fraction at the end of boil-off. This exponential behavior contradicts with their flat power-law envelope mass fraction scaling function. This greatly overestimates the initial mass fraction for low-mass and highly irradiated planets and underestimates the fraction for massive and cold planets in GSS18.

As a result, the combined effect from both the overestimated initial mass fraction and initial radius would correspond to boil-off initial conditions for low-mass (M_p ≤ 4 M_⊕) and/or highly irradiated (F ≥ 100 F_⊕) planets, leading to an exaggerated amount of mass loss, which is presented as a long-lived core-powered escape in GSS18 (see also Section 3.5 for a similar finding based on our numerical model and Section 4.5). This is evaluated directly based on the analytical model here. In Figure 15 we show evolution tracks for planets with the physics improved. In the left panels, the model with the original GSS18 initial conditions (red) overestimates both the post-boil-off radius and H/He mass fraction compared to the model with initial conditions that are self-consistently calculated from a boil-off phase (gray) for a 3 M_⊕ planet. On the other hand, the mass loss for a planet that is more massive and/or colder tends to be underestimated owing to the underestimated initial radius and mass fraction, which, however, has a less significant consequence, as these planets are less susceptible to mass loss. However, this yields different final mass fractions as shown in the right panels.

Zoom In Zoom Out Reset image size

4.6.3. Self-consistent Evolution

In the analytical work, the planetary radius is defined at the RCB, with the additional radiative atmosphere above ignored in the structure calculation. However, for a low-gravity planet, its scale height is often a large fraction of its total radius, especially with the variable surface gravity effect taken into account, which further increases the scale height at a higher altitude. Moreover, a sub-Neptune can readily possess a deep RCB at old ages, requiring a greater number of scale heights to reach its optical radius above the RCB. As a result, we suggest that its radiative atmosphere contributes up to 30%–40% to its total optical transit radius, far different from a giant planet, where the thickness of the radiative atmosphere can be typically ignored in modeling work. The importance of the radiative atmosphere becomes more marked for highly irradiated and low-mass planets. In Figure 16, we directly compare the RCB radius to the radius at 20 mbar, a typical optical transit radius. The difference in radius is significant. Without the radiative atmosphere, the location of the sub-Neptune population on an occurrence diagram shifts to the smaller radius end.

Zoom In Zoom Out Reset image size

The isothermal Parker winds from both our numerical model and the analytical model are sensitive to the choice of the mean molecular weight of the H/He atmosphere. For a less metal-rich atmosphere, yielding a smaller μ, the sonic point R_s = GM_p μ/2k_b T_eq penetrates deeper into higher-density regions of the radiative atmosphere. With the smaller R_s , the density ρ_s increases exponentially, which actually leads to a more vigorous mass loss $\dot{M}=4\pi {R}_{s}^{2}{\rho }_{s}{c}_{s}$ for the more metal-poor planet. This effect is especially marked when the planet is still young and inflated with a large scale height. In the GSS18 work, a pure H₂ atmosphere was assumed (μ = 2m_H), which we find is much more susceptible to mass loss compared to our default H/He atmosphere with μ = 2.35m_H. This is shown in Figure 17, in terms of the mass-loss timescale. Perhaps surprisingly, the pure H₂ atmosphere is completely stripped at least two orders of magnitudes faster, which could result in a more pronounced radius gap.

Zoom In Zoom Out Reset image size

Based on the GSS18 analytical model, a range of population-level work (S. Ginzburg et al. 2018; A. Gupta & H. E. Schlichting 2019, 2020, 2021; A. Gupta et al. 2022) has been done with similar physical treatments and initial conditions. With certain choices for core mass and equilibrium temperature distributions, they reproduce a radius valley and radius cliff very similar to that from observations. However, as we discussed above, with a number of physics and initial conditions improved, the final population from the analytical model can be potentially impacted. In this section, the planet occurrence rate distributions from these improved versions of the analytical evolution model are produced and directly compared to our numerical model on a population level, to reassess the importance of the post-boil-off bolometric-driven escape or "core-powered escape."

To evaluate how different approaches impact the demographics of small planets, we run a population of 10,000 randomly generated planets following the same equilibrium temperature distribution used in GSS18 that has a constant probability density ${dN}/d\mathrm{log}{T}_{\mathrm{eq}}$ for T_eq in the range of 500–1000 K and a power law of ${T}_{\mathrm{eq}}^{-6}$ within 1000–2000 K. For the core mass distribution, a variety of distribution functions are evaluated. The simplest one is a Rayleigh distribution that peaks at 3 M_⊕:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dM}}_{c}}\propto {M}_{c}\exp \left(-\displaystyle \frac{{M}_{c}^{2}}{2{\sigma }_{M}^{2}}\right),\end{eqnarray} \tag{ 35 }$$

with σ_M = 3 M_⊕, which truncates with a negligible probability density beyond 10 M_⊕. In GSS18, the authors also find that substituting a flat core mass function below 5 M_⊕ also yields a similar final distribution. To match observations, these authors find that a high-mass tail of core masses that extends out to 100 M_⊕ in addition to the Rayleigh distribution is needed to replicate the "radius cliff" feature. With this core mass and equilibrium temperature function, we compute the final mass fractions and planetary radii at 3 Gyr by different approaches, namely our numerical model-based approach, the fully analytical approach, and the semianalytical approach. The consequence of a variety of physics choices will also be evaluated in this section.

First of all, we implemented the GSS18 analytical evolution model with the same physics, initial conditions, and planet distribution functions. To validate our implementation, a single evolution curve was directly compared to the A. Gupta & H. E. Schlichting (2019) version of the analytical model, which has very slight differences between this version and that in GSS18 (A. Gupta 2024, personal communication). Our implementation yields identical evolution tracks. To further gain confidence, we ran a simple test distribution with 10,000 planet models with the same setup as in GSS18 for 3 M_⊕ planets with T_eq = 1000 K, producing a final distribution that is a clear match to Figure 1 in GSS18. In this test run, the initial envelope mass fraction follows an even distribution ranging between 10⁻⁵ and 100%, the same as in GSS18.

For a more comprehensive population-level comparison, using the same physical choices and planet distribution functions as documented in Section 5.1 and GSS18, our version of the analytical model reproduces a highly similar (but not identical) bimodal distribution to that from GSS18. The planet occurrence rate as a function of planet radii and H/He mass fractions with and without the mass loss is shown in Figure 18. Compared to their Figure 3 (the gray bars), our distribution shows a slightly higher sub-Neptune population occurrence in the range of 2.0−3.5 M_⊕, while the GSS18 distribution exhibits two nearly even peaks. Since the original model is not publicly available, the observed differences may stem from specific details in the realizations of the core mass distribution and equilibrium temperature distribution that we do not have access to.

Zoom In Zoom Out Reset image size

After having gained confidence in our implementation, we then step through a number of physical effects that we believe could be improved, as discussed in Section 4. With the same planet distribution functions employed as in GSS18, the computed final distributions of planetary radii and H/He mass fractions with different combinations of physics are shown in Figure 19. For all of the population-level comparison runs, we always use a mean molecular weight of the atmosphere for a H/He mixture with μ = 2.35m_H unless specified otherwise. The planetary radius is defined at the RCB, following GSS18.

Zoom In Zoom Out Reset image size

Our version of the original GSS18 setup is in yellow. Next, we add in the ${\dot{E}}_{\mathrm{adv}}$ energy-loss term that is missing from the GSS18 model, leading to a single peak around 2–2.5 M_⊕ without a radius gap feature (shown in red). This demonstrated that the GSS18 energy-limited mass loss ("core-powered escape") cannot form a radius gap with what we believe is more appropriate thermal evolution, which was suggested earlier in Equation (34), Section 4.4, and Figure 12. Then, with the ${\dot{E}}_{{\rm{adv}}}$ term, we switch to a decoupled isothermal Parker wind (the non-energy-limited wind unless otherwise stated) instead of the GSS18 coupled wind, which interestingly yields a distribution (black) nearly identical to that from the original setup in yellow. However, a significant timescale difference is seen on an evolutionary level in Figure 12, with the mass-loss timescale shifted to very young ages (∼1 Myr), corresponding to the boil-off timescale, with negligible mass loss thereafter.

Lastly, with the revised physics, we address the initial condition problem for the post-boil-off mass-loss phase by self-consistently including a boil-off phase in the analytical evolution model, which starts with a larger initial planetary radius and faster initial Kelvin–Helmholtz contraction timescale of 5 Myr. This blue histogram corresponds to our fully analytical approach, which we find leads to fewer planets with massive envelopes >10%, due to boil-off. The corresponding mass fractions are shown in the bottom panel. Although a pileup of blue planets in the range of 1.2–1.8 R_⊕ is caused by boil-off, this radius grouping has more planets than in the original yellow model, due to the larger sample of sub-Neptunes in the 10⁻⁴ to 10⁻² envelope range that never lose their envelopes, rather than making a larger super-Earth population (envelope-free) compared to the yellow model.

As discussed in Section 4.3, we find that the GSS18 analytical model yields an RCB radius that barely decreases throughout the whole boil-off phase, overestimating the mass-loss rate for low-mass planets, resulting in a deeper radius gap. The nearly constant RCB radius, however, is not seen in our numerical model. Moreover, since the total energy available for cooling estimated in Equation (31) is only valid when the thickness of the envelope is comparable to or smaller than the physical size of the core, extrapolating beyond this regime may lead to an evolutionary uncertainty. Therefore, the fully analytical model may not appropriately model the population.

To overcome the downside of the fully analytical approach, we investigated a semianalytical approach as described in Section 4.6.3. We set the initial mass fractions by interpolating in the post-boil-off results from our numerical boil-off model (Table 2). We then start the analytical evolution at the beginning of the long-lived bolometric-driven escape phase with boil-off properly incorporated. The initial radius is set by imposing an initial Kelvin–Helmholtz thermal contraction timescale of 50 Myr, which is the post-boil-off condition found in Section 3.3.

Additionally, we ran a grid of evolution calculations based on our numerical model with self-consistent physical conditions incorporating both boil-off and the long-lived escape, and we record the final mass fractions and 20 mbar optical radii at 3 Gyr as a function of planet core masses and bolometric fluxes. We then interpolate within the grid of evolution models for those 10,000 randomly generated planets to compute the final distributions, which represents our fully numerical approach.

In Figure 20, we compare the distributions of the final planet radii and envelope mass fractions at 3 Gyr from our fully numerical model (blue), the fully analytical approach (red; this was the blue distribution in Figure 19), and the semianalytical approach (black). Note that our numerical model employs a different core mass function, which does not have the high-mass tail of GSS18 but a simpler Rayleigh distribution. As a comparison, we also show the final distributions without any long-lived escape for both the semianalytical approach (yellow) and the fully analytical approach (gray), finding no significant changes compared to the black and red distribution with the long-lived wind, suggesting its insignificance.

Zoom In Zoom Out Reset image size

Importantly, we find no clear radius gap between 1.5 and 2 R_⊕ for the numerical distribution but a flat distribution down to 1 R_⊕, which we argue mostly results from boil-off. The planets that are on the left of the peak, in the range of 1.0–2.0 M_⊕, however, are a mixture of low-mass planets with very little H/He mass and bare cores (super-Earths) with a relatively smaller super-Earth population compared to the analytical models. Furthermore, the analytical approaches have a much greater population, with high initial mass fractions >10% largely contributing to the sub-Neptune population and constituting the radius cliff, which is a result of too many heavy planets >8 M_⊕ concentrated in the high-mass tail of the core mass function, as most of them remain intact after boil-off.

Additionally, the numerical model without the high-mass tail for the core mass distribution exhibits a cutoff at 5 R_⊕, shifted to a higher radius compared to the radius cliff feature in B. J. Fulton & E. A. Petigura (2018). This requires other mass-loss mechanisms to shrink planetary radii further, which suggests a role for XUV-driven escape. The radius cutoff we found is not a result of the massive planets >10 M_⊕, as is needed in the analytical models, but a consequence of our numerical boil-off. Lastly, we find that most sub-Neptunes have an intermediate H/He mass fraction in the range of 10⁻³ to 10⁻¹ after boil-off and by 3 Gyr.

Also in Figure 20, we find that the semianalytical approach (black) cannot form a radius gap but forms a peak around 1.5–2.0 R_⊕. It produces a H/He mass fraction distribution more similar to the numerical one, which has a larger population of planets that have intermediate H/He mass fraction 10⁻³ to 10⁻¹ compared to that from the fully analytical model. This population of planets corresponds to the peak at 1.5–2.0 R_⊕. As a comparison, the semianalytical approach still has the feature from the massive planet >10 M_⊕ tail, which forms another minor peak around 2.5 M_⊕. However, we argue that the peak around 1.5–2.0 R_⊕ accounts for the sub-Neptune peak around 2.0–2.5 R_⊕ from observations by B. J. Fulton & E. A. Petigura (2018) but is shifted to a smaller radius by about 30% owing to the lack of radiative atmospheres in the analytical model. Lastly, below 1.5 M_⊕ the super-Earths form not a clear peak but a gentle slope.

In summary, each of our models suggests that there is no significant long-lived bolometric-driven escape ("core-powered escape"). Boil-off is a powerful mass-loss process that greatly sculpts the planet population (blue) somewhat into the shape that we see today, but it alone is not able to reproduce the exact shapes of the observational features, including the deep radius gap and the more contracted radius cliff. Other mass-loss mechanisms such as XUV-driven escape are required to further shrink planet radii, deepen the radius gap, and shift the radius cliff to smaller sizes. Since the planets that are right beyond the observed radius cliff in the range 4–5 R_⊕ are sub-Saturn planets with more massive core masses >5 M_⊕, their mass loss is only subject to strong XUV heating (T. Hallatt & E. J. Lee 2022). We expect the hot (with an orbital period ∼a few tens of days) and less massive (<10 M_⊕) planets to eventually shift below the radius cliff, with the cold and more massive objects constituting the observed sub-Saturn population, which can be readily tested in future studies.

In Sections 4.7 and 5.3, we show the importance of the radiative atmosphere and the role of the core mass distribution. In this section, we evaluate the impact on the population distribution from these two choices.

In order to replicate the radius cliff, the cutoff in the frequency of planets with radii above 4.0 R_⊕, GSS18 used a power-law high-mass tail for the core mass distribution with a relatively large population of core masses >5 M_⊕, concatenated to the right of a Rayleigh core mass distribution. However, in our numerical model, we find that a 20 M_⊕ core by 10 Gyr has a large radius far beyond the "radius cliff." Moreover, these massive planets that are modeled remain unaffected after boil-off with a massive H/He envelope >15%. These larger-radius planets should be regarded as Neptune-class or larger and therefore excluded from a small-planet population study.

In Figure 21, we show the distributions from our numerical model in blue and the original GSS18 version of the analytical model in yellow. Our numerical model does not need the high-mass tail of core masses to present a radius cliff feature. On the contrary, including the high-mass tail in our numerical setup (shown in black) fails to replicate a radius cliff and unrealistically extends a flat radius distribution beyond 5 R_⊕, which suggests that the massive planets with H/He mass fraction >10% are unnecessary in our model, while it is important to the analytical models. Additionally, we demonstrate that the radiative atmosphere is a key to properly determining the shape of the planet distribution, due to its significant contribution of the physical size of a low-mass planet (see Section 4.7). As a comparison, we show a distribution of the RCB radii without the radiative atmosphere calculation based on the blue numerical model, which is shown in red. The radius peak shifts significantly to a smaller radius at 1.5–2.0 M_⊕. This peak location is consistent with that from the semianalytical model that does not a have radiative atmosphere (see black in Figure 20). It exhibits a similar cutoff at 4 R_⊕ as also seen in the analytical models (yellow), but without the need for the high-mass tail of core masses in the numerical work.

Zoom In Zoom Out Reset image size

Lastly, in Figure 22, we include the radiative atmosphere in the semianalytical model (black) and compare the distributions to the one without it (blue; same as black in Figure 20). We find that with the radiative atmosphere the peak shifts to a reasonable location at a larger radius, with a radius gap feature shown in the range of 1.3–1.7 R_⊕, although the envelope-free population <1.3 M_⊕ seems small compared to the observation. This is different from our numerical results that show a flat distribution (blue in Figure 20). The discrepancy is likely caused by the fact that the semianalytical model lacks the population with intermediate envelope mass 10⁻⁴ to 10⁻², which in comparison smooths out our numerical distribution. Therefore, the improved semianalytical model indicates that boil-off alone cannot reproduce the observed radius gap, consistent with our numerical results. Another mass-loss mechanism is required to further deepen the gap and shape the distribution into the one from observations.

Zoom In Zoom Out Reset image size

In our model, the boil-off atmospheric escape is assumed to be a non-energy-limited isothermal Parker wind at the planet's equilibrium temperature throughout the entire radiative atmosphere (see Equation (4)). As we discuss in Section 3.1, in the case of an energy deficiency, the internal energy of the wind that is quickly dissipated by the adiabatic expansion would lead to a cooler and lower-density deep atmosphere (below τ_h ) and would therefore hinder the hydrodynamic outflow. To tackle this problem in our model, the energy-limited prescriptions in Equations (22) and (27) are employed to assess the mass-loss rate, where the ratio of energy being turned into usable work is assumed to be 100%. Although it is shown to not significantly impact the post-boil-off mass fraction and therefore the later evolution after boil-off, this might introduce uncertainty into the early-stage evolution. This calls for a real hydrodynamic calculation to further study this issue.

Moreover, in reality, the wind is not completely isothermal as we assumed, such that the modestly warm deeper atmosphere has more bolometric energy deposited, whereas the optically thin region is a cooler isotherm known as the "skin temperature" region. Especially at a late age when the mass-loss rate has been greatly reduced, the wind base is pushed to a higher altitude, mostly located in the colder "skin temperature" region, resulting in a less efficient wind. This requires a detailed energy transfer model integrated into a hydrodynamic code to capture the physics discussed.

A wide range of sub-Neptune evolution models have assumed a simplified isothermal core that is set to have the same temperature as the bottom of the adiabatic envelope (N. Nettelmann et al. 2011; E. D. Lopez et al. 2012; S. Ginzburg et al. 2016) over the entire evolution, which we term a coupled core evolution, which corresponds to a perfectly conducting core. This simplification, however, can lead to a large uncertainty for the core evolution, as it may lead to different core cooling timescales (A. Vazan et al. 2018a), which can potentially impact our model results. This is another assumption that needs more attention, and we discuss three points below.

First, we discuss the consequence of our coupled core evolution. In reality, the core energy transport is limited by convection in the silicate outer layer (the geophysical "mantle") instead of conduction. The efficiency of energy transport depends on the physical conditions at the thermal boundary layer at the top of the mantle, e.g., the temperature contrast between the core temperature T_c and the temperature at the bottom of the envelope T_b , as well as the mantle viscosity. Therefore, our current coupled core evolution can potentially overestimate the core luminosity by instantly releasing all of the core thermal energy for cooling available, especially when the temperature at the bottom of the envelope T_b rapidly declines as a result of the vigorous boil-off mass loss. Note that such an overestimation typically does not affect our boil-off evolution, as the core luminosity can hardly impact the mass-loss process, due to the decoupling between the thermal evolution and the mass loss as demonstrated in Section 3.1.5.

Second, in the numerical model, we do not allow the core luminosity L_core to exceed the intrinsic luminosity L_int (see Sections 2.1 and 3.1.5), which may underestimate its role of thermally inflating a planet, potentially impacting the boil-off mass loss. This impact would typically occur in the transition phase, when the core luminosity can thermally inflate a planet faster than the mass loss decreases planetary radius in Equation (29) and with ${L}_{\mathrm{core}}\gt {L}_{\mathrm{int}}$ (i.e., β > 1). This would couple the mass loss with the thermal inflation, leading to more mass loss. We find that such thermal inflation can potentially affect a low-mass (≤4 M_⊕) and highly irradiated (≥100 F_⊕) planet, which are the ones more vulnerable to boil-off. To capture this behavior requires more sophisticated core evolution physics in future work, with a proper choice of the decoupled Parker wind as we suggested in Section 3.1 because of the coupling behavior. Additionally, we argue that the thermal inflation can only last for a brief period of time during boil-off, as vigorous mantle convection rapidly brings down the temperature contrast, terminating the thermal inflation. With the mantle convection considered, we expect a factor of order unity change of the final mass fraction after boil-off, rather than meaningfully affecting the evolutionary fate of planets (e.g., super-Earths or sub-Neptunes) that would otherwise correspond to an exponential change of the final mass fraction. Our treatment in this work therefore does not qualitatively impact our numerical results.

Third, if the initial core temperature is higher than that at the bottom of the envelope at the beginning of boil-off (∼8000 K), compared to what we assumed, this leads to a large thermal energy reservoir to cool and potentially more mass loss during the core-enhanced transition phase (see Section 3.1.5). We assess the consequence of using such a potentially underestimated initial core temperature (dictated by the initial temperature at the base of the envelope), by starting the evolution with a higher core temperature assessed with a decoupled core evolution such that T_c ≥ T_b . We evaluate the maximum possible mass-loss enhancement with the highest core temperature ${T}_{c,\max }$ that is theoretically allowed, such that the thermal energy ${{kT}}_{c,\max }$ is comparable to the gravitational binding energy of the core GM_c /R_c :

$$\begin{eqnarray}&&{{kT}}_{c,\max }=\displaystyle \frac{{{GM}}_{c}}{{R}_{c}},\end{eqnarray} \tag{ 36 }$$

which yields ${T}_{c,max}\sim \rm{18,000}\,{\rm{K}}$ (such a scenario assumes no core cooling or energy loss during planet formation and corresponds to an initial Bondi radius as small as the core size).

In Figure 23, we examine a decoupled core evolution (black) with a ${T}_{c,\max }$ start and show the evolution of temperature, H/He mass fraction, and RCB radius. The results are compared to the coupled core evolution (red) that starts with the lower core temperature assessed at the base of the envelope. To decouple the core evolution without implementing a convective/conduction transport scheme within the core (which is beyond the scope of this work), and to maximize the core-enhanced mass-loss effect without thermally inflating a planet and therefore the total mass loss, we assume that the core heating keeps the envelope at constant specific entropy and what would be a constant planet radius, meaning that any radius evolution is due only to boil-off mass loss. Under this assumption, the envelope is unable to thermally contract until the core temperature equals that at the bottom of the envelope (after which the core evolution transitions to being coupled). In this case, the core luminosity simply equals the envelope cooling luminosity, ${L}_{\mathrm{core}}={L}_{\mathrm{env}}$.

Zoom In Zoom Out Reset image size

In the presence of such an extreme initial core temperature, we find that the transition phase (which starts after the dotted black line, assessed with Equations (7) and (29)) between boil-off and the long-lived Parker wind is prolonged to at most 10 Myr. In the top panel, we show that before the transition phase the core temperature remains nearly constant, essentially meaning a negligible core-enhanced effect. The core temperature under the decoupled core evolution quickly declines and converges back to the temperature at the bottom of the envelope (gray) within the prolonged transition phase. This leads to a coupled core evolution thereafter, with only a small amount (∼0.1%) of mass-loss enhancement shown in the middle panel. As a comparison, our fully coupled evolution for the same model planet has a negligible mass-loss enhancement within the transition phase (middle) and a faster contraction rate (bottom) within 1–10 Myr as it exits the boil-off phase (the transition times are denoted with circle shapes; see Equations (7) and (9)), due to the thermal contraction. We note that such a mass-loss enhancement under the decoupled evolution and high initial core temperature is even smaller for a more massive planet and for a planet with a lower core mass, as boil-off would already remove its entire envelope. A similar behavior applies to a colder or hotter planet, respectively. Therefore, we argue that such decoupled core evolution with a much higher core temperature does not significantly affect our numerical model results.

In future work, a more detailed treatment of energy transport and core cooling is needed to better understand the interior and evolution of sub-Neptunes (A. Vazan et al. 2018b). For instance, an inner iron core could be partially liquid or solid and could crystallize over time. Similarly, parts of the rocky mantle could be liquid or solid, necessitating convective and/or conductive energy transport. Furthermore, there could be heat production from crystallization, as well as decompression as a result of boil-off. These processes can be better modeled in the future.

In a wide range of work, with mutually similar initial conditions and physical assumptions (S. Ginzburg et al. 2018; A. Gupta & H. E. Schlichting 2019, 2020, 2021; A. Gupta et al. 2022), a significant role for core-powered mass loss was found. We find that the results are in contrast with our model because their initial conditions would correspond to a boil-off condition, which was taken as core-powered mass loss. With an improvement of initial conditions, including a boil-off phase to self-consistently calculate the initial conditions for core-powered mass loss, and interior physics, W. Misener & H. E. Schlichting (2021) still found an efficient long-lived core-powered mass-loss phase, as a result of both the wind assumption, which is similar to that from GSS18 (Equation (12)), and the lack of the ${\dot{E}}_{{\rm{adv}}}$ term. Those treatments make their mass loss sensitive to the core luminosity and lead to an overestimation of the mass loss, discussed in Section 4.4.

In J. G. Rogers et al. (2024), this key assumption for the core-powered mass loss, however, was not assumed. But they instead used a similar non-energy-limited (decoupled) isothermal Parker wind assumption for both boil-off and core-powered mass loss, with similar initial entropies to those from our work. The energy source for the Parker wind was argued, without a detailed validation, to be from the interior cooling, which ultimately releases the gravitational binding energy when a planet contracts. However, this must correspond to a cooling timescale comparable to the mass-loss timescale, which we have shown to be unlikely to happen with these initial conditions (Section 3.2). In Figures 5 and 6 of J. G. Rogers et al. (2024), they still found a "core-powered" mass-loss phase. However, this phase was shown to be very short-lived, ≪1 Myr.

According to these authors, the "core-powered" phase results in a mass-loss enhancement by a factor of two (as seen in the top panel of their Figure 5), due to the core luminosity. This effect is more significant than what our model predicts but less impactful compared to other studies on core-powered mass loss that used the GSS18 energy-limited wind model.

We argue, however, that the reported difference is not solely due to core luminosity. The reasons are as follows. (1) The core-enhanced mass-loss effect should occur when the core luminosity maintains a larger planetary radius compared to the model without it. However, the planetary radii for the two model planets, with and without core luminosity, remain identical throughout the mass-loss phase, as shown in their bottom panel. (2) The second panel from the bottom indicates that the cooling timescale is much longer than the mass-loss timescale, suggesting that thermal evolution is inefficient, and thus the core luminosity is decoupled from the mass-loss process. (3) The initial mass fractions between the two models differ by a few percent despite having the same initial radii. This suggests that the two planets have different initial envelope entropies, with the model including core luminosity being initially hotter. We suggest that the difference in final mass fraction is due to the initial entropy effect, as demonstrated in our work, rather than a core-enhanced effect. (4) Similarly, at the end of boil-off, these two planets should have the same entropy owing to the decoupling effect if they start with the same initial entropy. Our model predicts that, with the same planetary parameters, their 2% model (without core luminosity) should be at least 20% larger in radius than the 1% model (with core luminosity). However, they found identical radii for both planets, implying different entropies. With the initial entropy properly set to the same value, we would expect a much smaller mass-loss enhancement in their model predictions.

Therefore, we consider this mass-loss phase as part of boil-off, which we termed the transition phase, due to this insensitive behavior to the core luminosity. More importantly, the previous work argued that core-powered mass loss can be separated from boil-off by two other characteristics, i.e., the energy to lift H/He mass and the relation to the thermal evolution. From these aspects, as demonstrated and reevaluated in Section 3.1, we do not think that such a mass-loss phase as found by J. G. Rogers et al. (2024) should be interpreted as core-powered mass loss, given their non-energy-limited wind assumption.

In spite of differences in the treatment of physics (i.e., core cooling, disk dispersal, and accretion) used in J. G. Rogers et al. (2024), they found qualitatively similar results to our model that are modestly different in the quantitative details. First, they used a different core luminosity formulation from ours, with a constant core cooling timescale over time to parameterize the core luminosity. Following the discussion in Section 6.2, a different core luminosity treatment did not qualitatively affect our evolution in boil-off, given the decoupling between the thermal evolution and mass loss.

Their model has a real disk dispersal phase and accretion phase, compared to our model, which we argue only slightly impacts our model results. In fact, a different disk dispersal time in their Figure 4 only changes the final mass fraction linearly, compared to the exponential behavior of the boil-off mass loss. The mass loss is controlled more by the initial entropy if the disk dispersal timescale is shorter than 1 Myr. Our initial entropy treatment yields a very comparable boil-off duration of 1 Myr and initial radius (entropy) to those from their work. In Section 3.2, we showed that our boil-off model is not sensitive to the initial mass fraction, and therefore the compositional (but not thermal) evolution history of the accretion phase is largely erased for a planet that has experienced a boil-off phase. Therefore, we find that the uncertainty from these different treatments is small.

Another physical process that needs additional attention is that of the internal energy advection. In Section 3.1.2, we suggested that the energy transport through the advective bulk motion is a consequence of the mass loss when the interior material is transported from a deep part of the envelope adiabat to a shallower region, which is eventually converted into the gravitational energy for the lifted material, instead of leading to an interior cooling. Even if such a cooling is assumed, we have demonstrated that it does not impact boil-off in Section 3.2. Here we discuss its potential role as an energy source to power the outflow in the radiative atmosphere.

This advective energy transport is completely inefficient in the assumed isothermal outflow in the radiative atmosphere. However, in reality, the deep radiative atmosphere is not completely isothermal, which gradually transitions from the adiabatic envelope T–P profile, rather than the sharp transition between the adiabatic profile and the isothermal radiative atmospheric profile that is assumed in our model. The deepest part of the escaping radiative atmosphere below the critical optical depth to visible photons τ_h , where it requires an additional energy source rather than the stellar heating, can be possibly located on a temperature gradient. In this case, the energy advection by the outflow can potentially play a role in draining the local internal energy to power the wind in this bottleneck region. Note that the energy advection in this case does not cool the interior as well. The efficiency of the advective energy transport depends on the temperature gradient between the RCB radius, where the H/He material is taken from, and the τ_h radius. Note that the discussion here is relevant only if τ_h occurs above the RCB in radius, which we find commonly happens to super-Earth-mass (>1 M_⊕) planets.

This effect would assist in the atmospheric escape in the bottleneck region. Therefore, given that the total amount of the internal energy advection and the intrinsic cooling energy is sufficient in the deep atmosphere to fuel the wind, the bottleneck of the boil-off mass loss then is the total stellar energy available to the radiative atmosphere (what we call bolometric limited; Equation (22)), instead of being intrinsic luminosity limited (Equation (27)). Understanding the relative importance of the internal energy advection and the radiative intrinsic cooling in the deep radiative atmosphere below τ_h requires a coupling between a radiative transfer model and a hydrodynamic wind model.

In our model, our treatments of the temperature profile of the radiative atmosphere and RCB always correspond to an intrinsic-limited boil-off that converges to common final physical conditions that are found at the end of boil-off. The subsequent long-term evolution is therefore unaffected if we switch between the different decoupled wind assumptions as discussed in Section 3.1.4. Consequently, our discussion above does not quantitatively affect our results.

We have shown that the decoupled isothermal Parker wind starts to transition into the XUV-driven phase right after boil-off, when the planet starts to efficiently thermally contract. At the same time, the post-boil-off bolometric-driven escape comes to be inefficient and be completely quenched shortly later. As discussed above (Section 5.3), boil-off alone generates a planet occurrence distribution that is somewhat similar to, but cannot match, the observed radius valley. A better match requires XUV-driven escape to further deepen the gap and shrink planetary sizes.

Moreover, in the next few megayears after the transition time, there should exist a period when the post-boil-off Parker wind and XUV-driven escape are comparably efficient. In this phase the hydrodynamic wind is launched by the absorption of intrinsic cooling energy from the deep atmosphere and then stellar bolometric energy above τ_h and the wind is accelerated further, and therefore the mass loss is enhanced by the heat deposited in the form of photoionizing high-energy flux at nbar pressure. Remarkably, such an enhanced isothermal Parker wind is coupled to the thermal evolution. We propose that these physical processes including boil-off, photoevaporation, and the enhanced escape should be further studied together, with a real hydrodynamic code and more complex energy transfer models in the core. This should also be integrated into the planetary thermal evolution framework in future generations of sub-Neptune evolution models.

With all of the improvements and physics discussed above included, one could generate a population of sub-Neptune models. The modeled demographics of the formed sub-Neptunes and super-Earths can be compared to the observed population and shed light on main physical processes sculpting their evolution.