Under the light of a new star: evolution of planetary atmospheres through protoplanetary disc dispersal and boil-off

The atmospheres of small, close-in exoplanets are vulnerable to rapid mass-loss during protoplanetary disc dispersal via a process referred to as `boil-off', in which confining pressure from the local gas disc reduces, inducing atmospheric loss and contraction. We construct self-consistent models of planet evolution during gaseous core accretion and boil-off. As the surrounding disc gas dissipates, we find that planets lose mass via subsonic breeze outflows which allow causal contact to exist between disc and planet. Planets initially accrete of order $\sim 10\%$ in atmospheric mass, however, boil-off can remove $\gtrsim 90\%$ of this mass during disc dispersal. We show that a planet's final atmospheric mass fraction is strongly dictated by the ratio of cooling timescale to disc dispersal timescale, as well as the planet's core mass and equilibrium temperature. With contributions from core cooling and radioactivity, we show that core luminosity eventually leads to the transition from boil-off to core-powered mass-loss. We find that smaller mass planets closest to their host star may have their atmospheres completely stripped through a combination of boil-off and core-powered mass-loss during disc dispersal, implying the existence of a population-level radius gap emerging as the disc disperses. We additionally consider the transition from boil-off/core-powered mass-loss to X-ray/EUV (XUV) photoevaporation by considering the penetration of stellar XUV photons below the planet's sonic surface. Finally, we show that planets may open gaps in their protoplanetary discs during the late stages of boil-off, which may enhance mass-loss rates.

To correctly model the consequences of photoevaporation and core-powered mass-loss on the small exoplanet population, one must carefully consider their initial conditions.It is likely that these processes are not the first time that small close-in exoplanets are at risk of losing atmospheric mass.Whilst protoplanetary discs survive for ★ E-mail: jamesrogers@ucla.edu∼ 3 − 10 Myr, the inner regions (≲ 1 AU) have been observed to disperse on timescales of ∼ 10 5 yrs (e.g.Kenyon & Hartmann 1995;Ercolano et al. 2011;Koepferl et al. 2013).This rapid disc dispersal has strong implications for the atmospheric evolution of close-in exoplanets since the sudden reduction in the disc gas' confining pressure results in dramatic atmospheric escape (Ikoma & Hori 2012).This process, in which a planet is illuminated by its young star for the first time, is referred to as 'boil-off' (Owen & Wu 2016) or 'spontaneous mass-loss' (Ginzburg et al. 2016;Misener & Schlichting 2021).
Population-level evidence for boil-off was hinted at in Rogers et al. (2023a), in which Kepler, K2 and TESS planets with well-constrained masses and radii were used to statistically infer the relation between planetary core masses and initial atmospheric masses at the end of the disc lifetime.Their results were found to be consistent with the predictions from Ginzburg et al. (2016), in which analytic arguments had been used to understand the atmospheric evolution of small planets prior and post-disc dispersal.Ginzburg et al. (2016) showed that boil-off will likely set the initial conditions for late-time escape mechanisms such as photoevaporation and core-powered mass-loss, which are thought to control small exoplanet demographics.It is, therefore, imperative to better understand the governing physics at play during boil-off.
In this paper, we investigate the effects of boil-off on the exoplanet population by explicitly modelling the transition from gas accretion to atmospheric escape, for a planet embedded in a dissipating gas disc.
Previously, boil-off was modelled in Ikoma & Hori (2012) under the assumption that planet atmospheres were in hydrostatic equilibrium with the surrounding disc.This approach typically resulted in planets losing their entire atmosphere during disc dispersal since the only way to remain in hydrostatic equilibrium with a vanishing disc is to host a vanishing atmosphere.At the point of disc dispersal  disp ∼ 3 − 10 Myrs, the cooling timescale for a planet's atmosphere is likely to have an upper limit of  cool ≲  disp since cooling timescales approximately track a planet's age and planets can form at any point during the disc's lifetime.Recall, however, that typical inner disc dispersal timescales are  disp ∼ 10 5 yrs; thus, a planet's cooling timescale  cool is likely greater than the disc dispersal timescale  disp .If this is the case, a planet may be unable to thermally respond to any changes in external pressure due to disc dispersal.It will undergo adiabatic expansion on timescales of order  disp .External heating (from the star and/or disc) will heat the expanding material, causing atmospheric escape via a rapid hydrodynamic outflow.
In Owen & Wu (2016), the case of non-hydrostatic equilibrium was explored with numerical evolutionary models.Based on the prior timescale argument, they assumed the disc dispersed instantaneously and material was removed from the planet via an isothermal Parker wind (Parker 1958).They found that up to ∼ 90% of atmospheric mass was removed on timescales of ∼ 10 5 yrs; however, this was based on initial conditions in which the planet's photospheric radius  ph was equal to its Bondi radius  B : where  p is the planet mass,  is the gravitational constant and  s is the isothermal sound speed, defined as: where  B is the Boltzmann constant,  eq is the planet's equilibrium temperature,  H is the mass of a hydrogen atom, and we use  = 2.35 (e.g.Anders & Grevesse 1989).The aforementioned assumption of  ph =  B meant that the planet was initially extremely inflated and fully convective, i.e. in a maximum entropy state, thus inducing considerable initial mass-loss rates.In reality, a planet will cool and accrete gaseous material before disc dispersal, thus forming a radiative outer zone that limits cooling (e.g.Rafikov 2006).Furthermore, whilst an instantaneous disc dispersal is warranted given that  disp <  cool , the introduction of a finite dispersal timescale may allow the planet to remain in hydrostatic equilibrium with the disc for at least some of its evolution.This possibility is backed up by the fact that the mass-loss timescales found in Owen & Wu (2016) were comparable to disc dispersal timescales.This paper addresses many of the assumptions adopted in the aforementioned works.We use a self-consistent numerical evolution models model for gaseous accretion and boil-off that accounts for a noninstantaneous disc dispersal process.Beyond this, we then develop our models to include core luminosities and study the subsequent transition from boil-off to core-powered mass-loss (e.g Ginzburg et al. 2016;Gupta & Schlichting 2019;Misener & Schlichting 2021).Finally, we consider stellar X-ray/EUV heating and the transition from boil-off or core-powered mass-loss to photoevaporation.We layout our numerical methods in Section 2, with results in Section 3, discussion in Section 4 and conclusions in Section 5.

METHOD
In order to model boil-off in a self-consistent manner, one must also consider the gas accretion phase that occurs before disc dispersal.Gas accretion is controlled by a planet's ability to cool (e.g. Lee & Chiang 2015).Hence, the initial thermodynamic state of a planet at the onset of boil-off is affected by accretion.In Owen & Wu (2016), boil-off evolution models were performed using Modules for Experiments in Stellar Astrophysics, MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018) ) and assumed planets undergo mass-loss via a hydrodynamic isothermal Parker wind (Parker 1958).These models did not, however, self-consistently model the gas accretion phase before disc dispersal and instead assumed planets began from an initial inflated state in which  ph ∼  B .This can be thought of as a young, fully convective, maximum entropy planet that is inflated due to an accretion luminosity from the core's formation.
As discussed, another assumption of Owen & Wu (2016) was to assume an instantaneous disc dispersal process.Whilst the disc dispersal timescale is typically much shorter than the disc's lifetime, the disc will require a finite amount of time to reduce its gas surface density such that the planet is exposed to a vacuum. 1 A gradual reduction in confining pressure will make the boil-off process less extreme since the planet has more time to respond thermodynamically.In this case, the planet may be able to remain in a quasi-hydrostatic equilibrium with the disc (implying causal contact), as was assumed in Ikoma & Hori (2012), for at least part of the planet's evolution.
The combination of these two assumptions: no gaseous accretion and hence maximally entropic initial planet, and instantaneous disc dispersal, imply that the mass-loss rates in Owen & Wu (2016) were likely upper limits.The presence of an upper radiative zone and a non-instantaneous disc dispersal will act to reduce mass-loss.

Atmospheric Structure Model
Similar to Owen & Wu (2016), we perform evolution models with MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018)).MESA is an open-source one-dimensional stellar evolution code that can handle certain hydrodynamic processes.More recently, it has been used to understand planet evolution (e.g.Owen & Wu 2013, 2016;Chen & Rogers 2016;Kubyshkina et al. 2020;Kubyshkina & Fossati 2022;Malsky & Rogers 2020;Owen 2020).For a small mass exoplanet with H/He dominated atmosphere, it uses the SCVH equation of state from Saumon et al. (1995) and dust-free opacity tables from Freedman et al. (2008).For each timestep in MESA, the total mass in the simulation domain is adjusted following any mass-loss/accretion (as described in Section 2.3).In addition to the timestep constraints within MESA, we add the requirement that the change in the atmospheric mass fraction is ≤ 0.1%.The system then relaxes by re-solving the stellar structure equations, which include additional hydrodynamic terms to account for advection.
An important limitation of MESA is that it performs radiative transfer using the diffusion approximation, requiring the atmosphere to be locally optically thick in the entire domain.Specifically, this requires the photon mean-free path to be much shorter than the pressure scale height for the entire domain.In addition, the code can only solve the stellar structure equations in sub-sonic flow since the implicit numerical scheme can not accurately resolve the necessary sound waves to model shocks and supersonic flow.As discussed in Section 2.3, we adopt a transonic outflow model that spans optically thin and thick regimes.In other words, the flow is launched in an optically thick region but extends beyond the photosphere.This limitation poses a severe challenge to modelling mass-loss.A solution exists, however, if one carefully considers where to place the outer boundary condition of the MESA model.

Outer Boundary Condition
When modelling boil-off and hydrodynamic escape in general, one must consider supersonic flow.As discussed, MESA sets two constraints determining where the simulation's outer boundary can be placed.Namely, the atmosphere below this point should be locally optically thick and subsonic.As described in the next Section, massloss is modelled via an isothermal Parker wind beyond this boundary, adding a third constraint: the atmosphere above the boundary should be approximately isothermal.To satisfy these constraints, we position the outer boundary condition at a fixed pressure  BC and temperature  BC .The choice in pressure is arbitrary, although it is chosen to ensure the gas is locally optically thick and in the upper radiative (and therefore approximately isothermal) layer.For all models presented in this paper, we adopt a value of  BC = 0.1 bar, but confirm that this choice does not affect our results.To ensure that this pressure satisfies the required conditions, we show that the corresponding boundary condition radius  BC at ( BC ,  BC ) is beneath the Bondi radius (subsonic), in the radiative layer (isothermal) and that the photon mean-free path is much shorter than the pressure scale height for the entire domain (locally optically thick) in Figure 1.
The outer boundary temperature, which is in the outer radiative layer of the atmosphere, is set to a constant equilibrium temperature  0 =  eq since it is in radiative equilibrium with the local protoplanetary disc and/or host star.In addition to a fixed temperature and pressure, MESA also requires the optical depth to out-going radiation at this boundary.This relies on knowing the density structure above  BC , which is determined by mass-loss/accretion, wich we discuss in the next Section.

Accretion and Mass-Loss
In order to model gas accretion and boil-off for a planet immersed in a protoplanetary disc, we require the ability to model both accretion and mass-loss self-consistently.As discussed, mass-loss in MESA is treated by adding/removing mass to the domain and then allowing the system to relax.2Given a set of boundary conditions (Section 2.2), one calculates and prescribes a mass-flux due to accretion/mass-loss.The question is, how to determine the mass-flux?
The distinction between mass-loss and accretion will depend on the physics of the disc.If the planet is pressure-supported by the local disc gas and can cool, then the atmosphere contracts, leading to an under-pressured system that cannot support the disc material above it, causing mass to flow onto the planet.On the other hand, if the disc's pressure support is insufficient, the upper atmosphere will expand and undergo mass-loss.We model the accreting/escaping atmosphere above the boundary condition with an isothermal layer, also set to  eq , allowing one to solve the appropriate isothermal fluid equations to determine the flow structure and hence mass-flux into/out of the MESA domain.The problem is set up as follows.As in In the upper panel, a temperature-pressure profile is shown for an atmospheric model with 5 ⊕ core and atmospheric mass fraction of ∼ 10% at ∼ 2 Myr.The boundary condition (BC) is set with  BC = 0.1 bar and  BC = 1000K, as shown with grey dashed lines and blue points.One can see that the temperature in the inner atmosphere is approximately adiabatic but flattens out in the outer radiative zone.Beyond this, the atmosphere is assumed to be isothermal.In addition, the location of the Bondi sphere is shown as a grey dot-dashed line.These properties imply that the MESA boundary condition is in an isothermal (radiative) and subsonic (interior to Bondi sphere) region.The lower panel shows the photon mean-free path,  ph = 1/( ) (shown in black-solid), where  is gas opacity, which is much shorter than the pressure scale height (black-dashed) for all radii, meaning that the entire MESA domain is also locally optically thick.
the original Parker wind problem (Parker 1958), one must solve the steady-state mass and momentum conservation equations: (3) where  is the planet mass-loss rate.We adopt an ideal gas equation of state for constant temperature and combine to find: which is substituted into Equation 4 to find the equation of motion for a steady state flow: This equation is trivially satisfied with  =  s .In this case, the lefthand side becomes zero, allowing us to solve for the location of this critical point.One finds that the flow velocity is equal to the sound speed when: which represents the transonic Parker solution (Parker 1958).Equation 6 is subject to an inner boundary condition provided by the outer boundary condition of the MESA domain at  BC , and an outer boundary condition which is controlled by the local gas density of the disc,  disc .Strictly speaking, the location of this boundary condition is min( B ,  Hill ), where the Hill radius is defined as  Hill = ( p /3 * ) 1/3 for semi-major axis  and stellar mass  * , since these are the points where gas dynamics are controlled by the disc (and therefore the host star) and not the planet.We simplify the mathematics of the problem by always adopting the Bondi radius  B as our outer boundary location, which is justified for the vast majority of our parameter space range if one relates a planet's equilibrium temperature  eq to semi-major axis  for a young, inflated pre-main sequence solar-mass star (as is done in Section 3.3). 3As such, the density at the outer boundary for this isothermal layer is set to  disc at  B .One can then relate  disc to the disc gas surface density Σ g : where  =  s /Ω K is the disc pressure scale height for a Keplerian angular frequency Ω K .As shown in Cranmer (2004), Equation 6 can be solved with Lambert W functions: where  is the function-branch and  () is: Here,  = 4 ln  B + 2 ln  crit − ( crit / s ) 2 + 4 is a constant and controlled by the fluid velocity,  crit , at the critical point,  B .If | crit | <  s then the solution is a subsonic breeze, whereas if | crit | =  s , then the solution corresponds to in-flowing transonic Bondi accretion (Bondi 1952) or an out-flowing transonic wind (Parker 1958).Examples of such solutions with varying critical velocity,  crit , are shown in Figure 2. Breeze solutions remain subsonic for all radii, whilst transonic solutions exceed the sound speed beyond/interior to  B .Note that for the transonic accretion solution,  6) are shown with radius.Here, six breeze solutions are shown, which remain subsonic for all locations ( / s < 1 for all / B ). Positive and negative velocities correspond to mass-loss and accretion, respectively.The two transonic solutions exceed the sound speed above/below the Bondi radius  B .Planets can accrete/lose mass in our adopted evolution models via any of these solutions.
the velocity increases in magnitude as one approaches the planet and would eventually become supersonic and shock inside the MESA domain.As discussed, MESA cannot treat a supersonic fluid.However, since we are not modelling planets that enter a runaway accretion phase, the accretion occurs via breeze solutions that are sufficiently subsonic to not be of concern.
The difference between requiring a breeze solution or a transonic solution depends on the fluid pressure at the outer boundary (e.g.Velli 1994;Keto 2020).In Owen & Wu (2016), only the transonic solution was considered due to the vacuum boundary condition and the inherent assumption of an instantaneous disc dispersal.In reality, the planetary flow will increase in launch velocity and transition from breeze solutions to the transonic wind as the disc gradually disperses.The faster this dispersal process happens, the quicker the outflow will tend to the transonic solution.In the case of hydrostatic equilibrium (HSE), the density profile in the isothermal layer beyond the MESA domain is given by: Strictly speaking, this is a solution to the fluid equations if  crit = 0. Now, if one manually decreases the density at the Bondi radius such that ( B ) <  HSE ( B ) (which in this problem is akin to reducing the gas density of the disc), then  crit will increase in an attempt to return to equilibrium.In the special case of an isothermal flow, the solution will remain as an out-flowing subsonic breeze if  −0.5 ≤ ( B )/ HSE ( B ) < 1 (Lamers & Cassinelli 1999).Suppose the density is reduced below this minimum value.In this case, the flow is no longer pressure-supported, the planet and disc lose their causal connection, and the only viable solution is a transonic wind.
For this problem, it follows that for a given ratio of  BC / disc there exists a solution to the fluid equations in non-hydrostatic equilibrium with a unique  crit ∈ [− s , + s ].From mass continuity in Equation 3 at the two boundaries, we have: Recall that there is only one unique solution to the isothermal fluid equations for a given  crit .Hence one can calculate  BC for a given  crit .Rearranging: we may now solve this equation to find  BC .Since  BC is a function of  crit , this equation is solved numerically with a brentq algorithm, adopting a numerical fractional tolerance of 10 −12 .Finally, the massflux  P due to Parker-type solutions (including subsonic breezes) at the MESA outer boundary is: Naively, one can adopt this value for each timestep of the MESA model; however, this results in large, unstable oscillations in planet mass.This is because the value calculated for  P in Equation 14 is an instantaneous mass-loss rate that is only appropriate for a set of boundary conditions at a given time.Recall, however, that MESA is an implicit numerical scheme that takes large time steps.The value of  P at the beginning of this timestep may, therefore, be very different to that at the end, particularly when noting that it is very sensitive to  BC / disc .This sensitivity means that the value is a dramatic overestimate of the required value averaged over the timestep, meaning that the planet loses too much mass and then requires accretion to return to pressure support.This unstable oscillation continues, resulting in unphysical mass fluxes.We adapt MESA with the following implicit mass-loss scheme to remedy this.Firstly, MESA takes a step forward in time with  = 0.The state of the system at the end of this timestep is then used to calculate  P , which represents the mass-loss that the Parker-type solutions can provide. 4Then, MESA goes back and repeats the same timestep with iterative values of  until it converges on a value such that the system has returned to hydrostatic equilibrium with the disc i.e. |(( B )/ disc ) − 1| < , where  = 10 −3 is the numerical tolerance.This mass-loss rate is labelled as  HSE .Finally, MESA performs the same timestep one more time, adopting the following mass-loss rate: In other words, if the subsonic Parker-type solutions can provide enough mass-flux to return to hydrostatic equilibrium within each timestep, then MESA adopts  HSE and remains in causal contact with the disc.If this is not the case, then the system cannot return to equilibrium and  P is adopted.Now that the mass-loss rate is determined for each timestep, the density structure beyond the MESA domain can be calculated, allowing one to evaluate the optical depth to out-going radiation as required in the MESA outer boundary condition (see Section 2.2).For each evaluation of the model at a given timestep, the adopted value of  provides  BC .Therefore, one can calculate the density profile in the isothermal layer from mass continuity: where () is calculated using Lambert W function for a given  BC .
The optical depth contribution from this isothermal layer is thus: Since MESA utilises the dust-free opacities from Freedman et al. (2008) in its atmospheric structure calculations, we also adopt the same opacities for the optical depth calculation in Equation 17for consistency, evaluated at  eq and ().This dust-free opacity choice is motivated by the fact that during inner disc dispersal, dust rapidly drains from the inner disc due to gas drag (e.g.Alexander & Armitage 2007;Owen & Kollmeier 2019).We highlight that changing the gas opacities, such as the inclusion of dust, will act to change planetary cooling rates and thus the amount of atmospheric mass that is accreted onto the planet (e.g.Ikoma et al. 2000;Piso et al. 2015).Regardless of this fact, however, boil-off will still operate since the process is triggered by a change in a planet's outer boundary condition, regardless of the exact amount of accreted gas.The contribution to optical depth from the disc is approximated to be  disc ≈  disc Σ g /2 where  disc = 0.01 cm 2 g −1 for a dust-free environment (e.g.D' Alessio et al. 2001).Note that this choice in disc opacity has little effect on the evolution, although we comment on the potential effects of dust in Section 4.7.Finally, the total optical depth to out-going radiation beyond the MESA boundary is: In addition to contributing to optical depth, the isothermal layer contributes a small amount of atmospheric mass to the planet.Recall that the outer boundary  BC for the MESA domain is defined at a fixed temperature and pressure, rendering it meaningless in terms of the physical mass of the planet.The true mass is calculated as the sum of the mass inside the MESA domain and the mass in the isothermal layer out to the Bondi radius.Note that our model differs from the work of Valletta & Helled (2020) due to our outer boundary existing at fixed pressure and temperature.In the latter study, accretion was modelled for growing giant planets by considering an outer boundary at a fixed radius, with mass flux prescribed to ensure the latter condition is satisfied.As a result, there was no direct consideration of the fluid equations that govern accretion.Our approach allows for a self-consistent view of the interaction of a planet's atmosphere with its local protoplanetary disc, including the transition to mass-loss for low-mass planets.

Initial MESA models
Constructing initial models in MESA for young planetary atmospheres is a convoluted process.Constructing a low-mass, low-luminosity and yet highly inflated body with a solid core requires multiple steps, which are outlined below and follow previous works of Owen & Wu (2013); Chen & Rogers (2016); Kubyshkina et al. (2020); Kubyshkina & Fossati (2022); Malsky & Rogers (2020).
(i) Initial model: To begin, the create_initial_model routine is used to produce an initial H/He dominated model of 2 Jup and 6 Jup .We choose a hydrogen mass fraction of 0.74, helium mass fraction of 0.24 and metal mass fraction of 0.02 to reflect typical nebular conditions.This model is then evolved for 10 12 yrs such that the planet has a low enough entropy in order for atmospheric mass to be removed stably in Step (iii).This allows one to produce a planet with the desired initial atmospheric mass fraction.
(ii) Add planetary core: Next, the relax_core routine is used to add an inert core of mass  c .To determine the core radius  c , we adopt the mass-radius relations of Fortney et al. (2007) for a composition of 1/3 iron, 2/3 silicate.Note that MESA treats the core as inert and does not model its interior.In section 3.2, we introduce core luminosity due to cooling and radioactive decay.
(iii) Choose initial atmospheric mass fraction: To choose an atmospheric mass fraction , the relax_mass routine is used to reduce the initial 2 Jup model to a desired value.Note, however, that this atmospheric mass fraction is arbitrary at this point since we allow the system to come into equilibrium with the disc and then accrete material prior to disc dispersal.The initial atmospheric mass fraction at the end of Step (v) is of physical interest to the problem.
(iv) Impose outer boundary condition: At this point, we switch from the standard photospheric boundary condition implemented in MESA (see Paxton et al. 2011) to one with fixed pressure and temperature as described in Section 2.2.This boundary condition is prescribed with the other_surface_PT routine.At this point, we also perform tests to ensure this boundary is in an optically thick, subsonic and isothermal region (see Figure 1).
(v) Inflate and equilibrate: To simulate a young planet in equilibrium with its nascent disc, one must increase the model entropy from the previous steps.We use the other_energy routine to introduce a core luminosity that acts to inflate the planet, which can be interpreted as the accretion luminosity,  acc , that would have been present during core formation.We iterate the core luminosity until the density at the Bondi radius is equal to that of the local disc density |(( B )/ disc ) − 1| <  where  = 10 −3 is the numerical tolerance.
(vi) Turn off core luminosity and accrete: The final preparation step is to switch off the core luminosity from Step (v), reset the planet's age to a specified value  pl within the disc's lifetime, and allow gas accretion to commence, as described in Section 2.3, which is implemented in the extras_check_model.We additionally switch on the hydrodynamic fluid terms with change_v_flag to allow for advection.The model can now self-consistently switch from accretion to mass-loss as the disc begins to disperse.At this point, which marks the onset of boil-off, we record the atmospheric mass fraction as  init .
The resulting model is representative of a young planet that has formed and accreted a H/He dominated atmosphere while immersed in a protoplanetary disc.The planet parameters at this point are the core mass,  c , initial atmospheric mass fraction,  init , cooling timescale  cool , equilibrium temperature,  eq , semi-major axis, , and planet age  pl .Note that this latter variable is, in fact, a crude estimate of the actual planet age, but instead represents the amount of time for which we have allowed accretion to take place.In essence, it is a variable that controls the initial atmospheric mass fraction of a planet before dispersal begins.

RESULTS
Our results are split into three parts: in Section 3.1 we aim to isolate the physics of boil-off and thus consider minimal disc models with planets at constant equilibrium temperature and no luminosity contribution from the core.In Section 3.2 we introduce core luminosity and consider the transition from boil-off to core-powered mass-loss and/or XUV-driven photoevaporation.Finally, in Section 3.3, we employ more realistic disc models to investigate the effect of a changing equilibrium temperature throughout disc evolution.

Pure boil-off: the race between planet cooling timescale and disc dispersal timescale
To begin, we evolve the local disc gas surface density Σ g with the following parameterisation: where Σ g,0 is the initial gas surface density,  disp is the time at which the disc starts to disperse and  disp is the dispersal timescale.To investigate the effect of disc dispersal rate, Figure 3 shows the evolution of a 5 ⊕ core at 0.1AU, with a constant equilibrium temperature of 900K.The disc starts dispersing at  disp = 3 × 10 6 yrs, but with three separate dispersal timescales:  disp = 10 4 , 10 5 and 10 6 yrs.In addition, we begin planet evolution at different times: 1 × 10 6 , 1.4 × 10 6 and 2×10 6 yrs.As a result, these planets have a range in atmospheric mass fractions at the point of disc dispersal since they have had varying amounts of time to cool and accrete.Figure 3 reveals that, despite having initial atmospheric mass fractions between 30 − 50% at the onset of disc dispersal, the final atmospheric mass fraction is significantly reduced by greater than ∼ 90%.Planets that undergo dispersal timescales of  disp = 10 4 , 10 5 and 10 6 yrs are left with ∼ 2%, ∼ 2.3% and ∼ 4% atmospheric mass fractions respectively.Thus, for a given core mass, equilibrium temperature and dispersal timescale, planets all host approximately the same atmospheric mass fraction, independently of how much atmosphere they initially accreted.The key quantities controlling mass-loss are the disc dispersal timescale  disp , and the atmosphere's cooling timescale  cool (which is a proxy for planet entropy).This latter variable is important since it controls the phase of accretion that occurred prior to boil-off and ultimately sets the ability of the planet to thermally respond to the changing disc boundary condition.In the case of no internal energy sources (such as a cooling planetary core, as considered in Section 3.2.1), the cooling timescale is the time it would take to radiate the total energy  tot of the atmosphere away, with contributions from internal thermal energy and gravitational binding energy: where  i is the internal thermal energy, (< ) is the mass interior to radius , and  is the radiative luminosity of the planet, evaluated at the radiative-convective boundary.If the ratio of  cool / disp is large, the planet cannot thermally respond to the disc dispersing, promoting adiabatic expansion and increased mass-loss rates.In the bottom panel of Figure 3 we also show the cooling timescale of the planet, which displays a sudden drop at the point of boil-off due to the sharp decrease in atmospheric mass.In Figure 3 we also note the ratio of cooling timescale to dispersal timescale, recorded at the time of disc dispersal onset at 3 × 10 6 yrs.The larger this ratio, the more mass is lost.We note that in the case of very rapid disc dispersal  disp = 10 4 yrs, planets that accreted more atmosphere retain a slightly more massive atmosphere by the end of boil-off.
A different regime exists, however, if the ratio of cooling timescale to dispersal timescale becomes less than unity.In this case, the planet can remain in thermodynamic equilibrium with the disc, and would not undergo the significant mass-loss and associated cooling processes involved in boil-off.Whilst a small amount of accretion would continue after the disc begins to disperse, the planet would then release this atmosphere to stay in equilibrium with a vanishing disc.This is indeed the case considered in Ikoma & Hori (2012), which resulted in very small atmospheres after dispersal concluded.One can see that in the case of the longest dispersal timescale in Figure 3, a small amount of accretion occurs after the onset of disc dispersal at 3 Myrs.
The results in Figure 3 suggest that in the case of a long cooling timescale compared to dispersal timescale (as is suggested by the disc observations), boil-off sets the initial conditions for the post-disc thermal evolution of close-in exoplanets.To understand this further, Figure 4 shows the evolution of Mach number at the Bondi radius and entropy at the base of the atmosphere for the case of planets beginning evolution at 2×10 6 yrs (green lines) from Figure 3.The Mach number evolution demonstrates that the disc dispersal timescale controls the boil-off process.Initially, boil-off presents as a low Mach number outflow while the planet and disc can maintain causal contact via an instantaneous pressure-support mediated through mass-loss (see Equation 15).However, as the disc density drops further, the Mach number must increase to maintain pressure support at a rate that is controlled by the dispersal timescale.Eventually, the Mach number must rise so high that a causal connection cannot be maintained and a transonic wind is launched with unity Mach number.
Intuitively, the disc dispersal rate must control the outflow speed since the ratio of planet density and disc density controls the Mach number and therefore mass-loss rate (see Equation 13).If a planet were to have too much mass removed at a given time, then massloss would slow until the disc disperses to a sufficient level.Massloss would then resume at the rate determined by disc dispersal.Therefore, this highlights a significant result from this work: breeze solutions allow a planet to remain in causal contact with the local disc, allowing the disc dispersal process to dictate the final atmospheric mass fraction.In essence, this then builds upon the separate effects from Ikoma & Hori (2012) and Owen & Wu (2016), which we discuss in Section 4.2.disc dispersal show an increase in  rcb at the onset of dispersal.This effect is due to a brief period of d expansion, which sets up a steep temperature gradient, thus enhancing radiative cooling.As previously highlighted, the magnitude of the d expansion and, thus, radiative cooling is based on the balance of cooling timescales and dispersal timescales.For planets with a larger ratio of  cool / disp , the greater the magnitude of d expansion since the initial evolution is better approximated as adiabatic.Nevertheless, and as initially found in Owen & Wu (2016), we confirm that advective cooling is the dominant source of energy transport deep in the atmosphere for rapid disc dispersal.

The interplay between boil-off, core-powered mass-loss and XUV photoevaporation
Whereas in Section 3.1, we did not consider any luminosity emanating from the planetary core; in reality, the base of the atmosphere will interact with a young molten surface and cool as a coupled system.To address this, we follow Lopez et al. (2012); Chen & Rogers (2016); Vazan et al. (2018) and assume the core is isothermal with temperature   , which is initially set equal to the base temperature of the atmosphere (which is determined through MESA's structure model).We quantify the core cooling luminosity  cool as: where we assume a heat capacity appropriate for a silicate mantle  V = 1.0 J K −1 g −1 (e.g.Alfè et al. 2002;Valencia et al. 2010).
Numerically determining d c /d is not straightforward during rapid changes in an atmospheric structure such as those involved in boil-off.
The core acts as a significant reservoir of thermal energy due to its large heat capacity, meaning that only a small change in core temperature  c over the timescales that boil-off acts can deposit significant energy into the atmosphere, which poses a complex numerical challenge (particularly for an implicit integration scheme).To address this, we exploit the notion that the magma surface is very unlikely to remain in thermal equilibrium with the base of the atmosphere on the very short timescales that boil-off operates ≲ 10 5 yrs.In fact, Markham et al. (2022); Misener & Schlichting (2022); Misener et al. (2023) have shown that atmospheric convection may be inhibited near the core due to high mean-molecular weight gradients, leading to an interior radiative zone.One can speculate that the inefficient cooling of such a radiative zone may provide a mechanism for a planet's core and atmosphere to cool out of equilibrium.Nevertheless, in the absence of constraints on the equilibration timescales of magma-atmosphere boundaries, we parameterise d c /d with a Newtonian cooling term: where  atm,base is the temperature at the base of the atmosphere and  c is the thermal equilibration timescale of the core-atmosphere boundary.During periods of dramatic atmospheric cooling due to massloss, we adopt core temperature evolution according to Equation 22.
Once the core and atmospheric base have returned to equilibrium, we set  c =  atm,base such that d c /d = d atm,base /d for the remaining evolution.We find that unphysical numerical solutions exist for  c ≲ 10 4 yrs, although we do not imply that physical solutions do not exist below this value, merely that our numerical scheme cannot find them.In the range 10 4 ≲  c ≲ 10 7 yrs, we find that atmospheric evolution is weakly sensitive to  c .We therefore set the equilibration timescale to  c = 10 5 yrs for all models.We note, however, that if  c ≳ 10 7 yrs, then the core-atmosphere boundary remains out of equilibrium until after boil-off has concluded, leading to a different evolution.In this case, we reproduce the results of Vazan et al. (2018) with delayed planet contraction over timescales of  c .
In addition to core cooling, heating also acts as a result of radioactive decay in the core  radio , for which we consider the contributions of 238 U, 235 U, 232 Th, 40 K and 26 Al, which have half-lives of 4.5 Gyrs, 0.7 Gyrs, 14.1 Gyrs, 1.3 Gyrs and 0.7 Myrs respectively.We adopt meteoritic abundances for all nuclides from Anders & Grevesse (1989), apart from 26 Al, which is taken from Lee et al. (1976); Lee et al. (1977).This latter isotope dominates radiogenic heating during boil-off due to its short half-life.Nevertheless, the cooling luminosity is typically orders of magnitude larger than that of radioactivity.The total core luminosity is given by: The effect of core heating is to increase mass-loss rates since the additional luminosity supplies more material to the base of the outflow (or, in other words, increases the atmosphere density at the base of the outflow).

Transition from boil-off to core-powered mass-loss
The introduction of core luminosity allows us to consider the process of core-powered mass-loss in detail (Ginzburg et al. 2018;Gupta & Schlichting 2019, 2020, 2021;Gupta et al. 2022;Misener & Schlichting 2021).It is crucial to note that the physics of boil-off and corepowered mass-loss are fundamentally linked.Both processes are driven by bolometric heating of a planet's upper atmosphere, with mass-loss throttled by the rate at which material can be removed from the planet via a thermally driven Parker wind/breeze.However, energy is required to transport material from deep in the planet's gravitational well to the base of the outflow.During boil-off, transport energy is supplied via advection, as well as rapid cooling and contraction of the planet's atmosphere, which thus releases binding energy.On the other hand, core-powered mass-loss relies on core luminosity to transport material to the base of the outflow.In other words, the removal of mass is not 'core-powered'; instead, the core luminosity allows material to be supplied to the top of the atmosphere, wherein it is heated and physically removed by the host star.This effect is precisely why core-powered mass-loss predicts a radius  24) with solid lines.Note that  atm is initially positive and then changes sign to be negative ∼ 9 × 10 5 yrs after disc dispersal began.The core luminosity, including contributions from core cooling  cool and radiogenic heating  radio , is shown in dot-dashed and dotted, respectively.In the third panel, the cooling timescale is shown, as defined in Equations 20 and 25 for the case of with and without core luminosity, respectively.Finally, planet radius evolution is shown in the bottom panel, with the photospheric radius and radiative-convective boundary shown in solid and dot-dashed, respectively.The transition between boil-off and core-powered mass-loss is defined as the point at which the core luminosity becomes greater than the rate of change of atmospheric energy  core ≥  atm .This point is shown for each model with coloured squares.
Figure 5 demonstrates the effects of adding core luminosity.We run evolution models with a 5 ⊕ core with constant equilibrium temperature of 900K, disc dispersal onset at  disp = 3 × 10 6 yrs and dispersal timescale of  disp = 10 5 yrs.We compare the scenario with and without core luminosity in black and red, respectively.In the top panel, we show the evolution of atmospheric mass fraction, highlighting that allowing the core and atmosphere to cool as a coupled system enhances mass-loss.
In the second panel, we show the rate of change in atmospheric energy as solid lines, with core luminosities from cooling and radiogenic heating in dot-dashed and dotted, respectively.We demonstrate that core cooling dominates over radiogenic heating, as found in similar works (e.g.Lopez et al. 2012;Chen & Rogers 2016).Note that the definition of the rate of change in total atmospheric energy is as follows: In other words,  atm is the rate of change in total atmospheric energy, which includes change in internal energy, where  i is the internal specific energy; change in binding energy, where (< ) is the mass interior to radius ; and change in mass, where  pl is the total mass of the planet (which, crucially, changes with time) interior to the Bondi radius.This definition differs from the planet luminosity shown in Figure 4, representing the radiative luminosity released at the radiative-convective boundary.In the case of no mass-loss, the two definitions are identical since energy is only lost due to radiation.During mass-loss, the planet loses negative gravitational potential energy, hence  atm > 0. Once mass-loss ceases, the planet loses energy due to cooling and gravitational contraction.As a result, the  atm < 0 after ∼ 9 × 10 5 yrs, which can be seen as a change in sign in Figure 5.
With the planetary energy budgets in mind, we choose to define the transition between boil-off and core-powered mass-loss as the point at which an atmosphere's rate of change in energy (as defined in Equation 24) becomes less than the luminosity released from the core.In order words, core-powered mass-loss begins when the core dominates the planetary energy budget.Mathematically, this transition occurs when  core ≥  atm , which we highlight in Figure 5 with a square marker.Note, however, that the bottleneck for both boil-off and core-powered mass-loss is the same, being the Parkerwind mass-loss rates from Equation 14.Under some circumstances, the energy liberated from the core can provide the limit for mass-loss (e.g.Ginzburg et al. 2018); however, we do not find any examples of such cases.As a result, we do not expect the behaviour of atmospheric escape to change dramatically during the transition from boil-off to core-powered mass-loss.
The third panel of Figure 5 shows the cooling timescale for the planet's atmosphere with and without the influence of core luminosity.In the case of planetary evolution with core luminosity, the definition of cooling timescale changes from that of Equation 20, to: where we have additionally accounted for the thermal energy reservoir stored in the planetary core  c .Note that here we approximate the core as incompressible (as is assumed in our MESA model), mean-ing that the change in binding energy of the core does not contribute significantly to the thermal evolution of the planet.As with Equation 20,  atm is the atmosphere's total energy and  is the radiative luminosity at the radiative-convective boundary.Initially, the cooling timescale is larger in the case of a cooling core (red line) since the planet has a larger energy reservoir.As mass-loss continues, however, the planet hosts a less massive atmosphere and has a higher radiative luminosity due to the contributions from the core, meaning that the cooling timescale drops to lower values during boil-off than in the case of no core-cooling.Once mass-loss has concluded, we still find that planets have prematurely cooled and display a longer cooling timescale than their age, as identified in Owen & Wu (2016).Finally, the bottom panel of Figure 5 shows the planet's photospheric radius and radiative-convective boundary evolution in solid and dot-dashed lines, respectively.Note, however, that prior to the disc completely dispersing (coinciding with ∼ 10 6 yrs after disc dispersal began in Figure 5), the photospheric radius of a planet is physically meaningless.

Transition from boil-off or core-powered mass-loss to XUV photoevaporation
Whereas boil-off and core-powered mass-loss rely on bolometric heating of escaping atmosphere, photoevaporation relies on Xray/EUV (XUV) photons to drive hydrodynamic outflow (e.g.Erkaev et al. 2007;Murray-Clay et al. 2009;Owen & Jackson 2012).The transition between these heating regimes is characterised by the depth of XUV photon penetration into the planet's atmosphere.As shown in Owen & Schlichting (2023), photoevaporation dominates over a bolometric-heated wind (such as boil-off or core-powered mass-loss) when XUV photons penetrate below its sonic surface of the bolometrically heated region at  B .This allows the gas to be heated to typical temperatures of ∼ 10 4 K in the subsonic (and thus causally connected) portion of the outflow.To quantify this effect, we follow Owen & Schlichting (2023) in determining the XUV penetration radius  XUV , which occurs at typical column (number) densities of  ≈ 10 22 cm −2 (e.g.Ercolano et al. 2009) for soft X-rays, which dominate the mass-loss at early ages (e.g.Owen & Jackson 2012;Owen & Wu 2013).Hence, we state that the following is true at  XUV : where () is the gas number density between the planet and star, including contributions from the atmosphere and protoplanetary disc.
Beyond  XUV , we assume the outflow is isothermal with a temperature of  XUV = 10 4 K and accelerated to the sound speed  s,XUV = ( B  XUV / H ) 1/2 .In order to determine the density at the boundary between bolometric and XUV-heated regions, we apply conservation of momentum at  XUV , yielding: where subscripts 'XUV' and 'bol' refer to variables in the XUV and bolometric heated region, respectively.Here, we have neglected the ram pressure term in the bolometric heated region since the flow is sufficiently subsonic ( bol / s ≪ 1).We also assume the flow in the XUV heated region is transonic at  XUV , hence  XUV ( XUV ) ≈  s,XUV .The density profile in the XUV heated region beyond  XUV is calculated with: Time Since Disc Dispersal Began (yrs) Figure 6.Atmospheric mass fraction evolution is shown for planets with core masses of 2 ⊕ , 4 ⊕ , 6 ⊕ and 8 ⊕ (moving from left to right) at constant equilibrium temperatures of 600K, 900K, 1200K, 1500K in blue, green, orange and red respectively.Time is measured from the point at which disc dispersal began, in this case at  disp = 3 × 10 6 yrs.We assume a disc dispersal timescale of  disp = 10 5 yrs.Lines demonstrate the transition from boil-off driven mass-loss (solid lines), to core-powered mass-loss (CPML; dashed lines, with transition denoted with squares) and XUV-driven photoevaporation (grey dotted, with transition denoted with circles).Note that although we mark the transition to photoevaporative mass-loss as circles in the models, we do not include the relevant photoevaporative mass-loss rates after this point in our models.
where, in order to calculate  XUV (), we set  =  s,XUV at  XUV in the general solutions to isothermal outflow from Equations 9 and 10.Equations 26, 27 and 28 are solved simultaneously with a brentq solver at each timestep of our evolution calculations to find  XUV , adopting a numerical fractional tolerance of 10 −12 .As in Owen & Schlichting (2023), we assert that the planetary outflow is dominated by photoevaporation over boil-off or core-powered mass-loss if  XUV <  B .Note that explicitly including the photoevaporativedriven mass-loss rates beyond this transition is beyond the scope of this work and are therefore not included in our evolutionary calculations.Instead, we record planet properties at the point of XUV penetration and consider these as the initial conditions for photoevaporation models.
Figure 6 shows a small parameter study for planets with core masses of 2 ⊕ , 4 ⊕ , 6 ⊕ and 8 ⊕ at temperatures of 600K, 900K, 1200K, 1500K.Under our assumptions of dust-free gas, planets further out in the disc will be at lower equilibrium temperatures and can, therefore, cool more efficiently, allowing accretion to occur at a quicker rate.In the same vein, a boil-off process will be less dramatic at lower equilibrium temperatures since it is the local irradiation that actually removes the material from the planet's influence (note the dependence of sound speed and thus the temperature in outflow velocity from Equation 10).However, accretion and boil-off are also strongly controlled by the core mass, since this sets the gravitational potential and, thus, the ease with which mass can enter/leave the planet.A clear trend can be seen in Figure 6 in which planets with cooler equilibrium temperatures accrete a more massive atmosphere.Similarly, planets with larger mass cores follow the same trend.Indeed, the 8 ⊕ model at 600 K undergoes runaway accretion before disc dispersal onset, terminating the evolutionary model, and is not shown in Figure 6.
In Figure 6, we demonstrate the transition from boil-off driven mass-loss (solid lines) to core-powered mass-loss (dashed lines, with transition denoted with squares) and XUV-driven photoevaporation (grey dotted, with transition denoted with circles).Figure 6 demonstrates that boil-off can perform mass-loss at early ages for all considered masses and equilibrium temperatures.As discussed in Section 3.2.1, the energy/work done in driving atmospheric mass-loss dur-ing this boil-off phase (solid lines) is dominated by the release of binding energy due to atmospheric contraction and not contributions from the core.
Recalling that the vast majority of mass-loss under boil-off/corepowered mass-loss occurs while the protoplanetary disc is present and optically thick, these results suggest that some small-mass, closein planets may have been stripped before the disc has completely dispersed, independently of XUV photoevaporation.We discuss this point in Section 4.5.
We present the properties of planet models from Figure 6 at the onset of XUV photoevaporation in Figure 7. Since these models adopt a disc dispersal timescale of  disp = 3 Myrs and XUV photons penetrate beneath the Bondi sphere ∼ Myr later, these conditions represent the planets at ∼ 4 Myrs.Note that the exact time of XUV photoevaporative dominance would change for different disc ages and dispersal timescales.In the left-hand panel, we show atmospheric mass fractions compared to the results of Owen & Wu (2016) in green squares.We discuss the differences in Section 4.2; however, we highlight that the models of Owen & Wu (2016) adopt different initial atmospheric mass fractions as a result of their different assumptions, many of which we have individually addressed in this work.
In the right-hand panel of Figure 7, we show the photospheric radius of each planet at the onset of XUV photoevaporation.For planets stripped of their hydrogen-dominated atmosphere, their photospheric radius equals their core radius, as highlighted with the mass-radius relation for Earth-composition cores from Fortney et al. (2007) in grey dashed.Since contributions from the core dominate planet mass, this panel can be interpreted as a mass-radius distribution for young planets at the end of disc dispersal and the onset of XUV photoevaporation.Planets at lower equilibrium temperatures are larger in size since they host more atmosphere post disc-dispersal.Note, however, that the exact size of a planet at a given time will depend on when the disc disperses and the rate at which this occurs.

Physical Disc Model
The final development of our evolution model is to consider a more realistic disc evolution sequence with a varying planet equilibrium Earth composition T eq = 600K T eq = 900K T eq = 1200K T eq = 1500K Initial Conditions for XUV Photoevaporation ( ∼ 4 Myrs) Figure 7.A summary of planet properties at the onset of XUV photoevaporation from models in Figure 6 with constant equilibrium temperatures of 600K, 900K, 1200K, and 1500K in blue, green, orange and red respectively.We assume the disc disperses at  disp = 3 × 10 6 yrs with a disc dispersal timescale of  disp = 10 5 yrs, all of which lead to XUV photoevaporation dominating at ∼ 4 Myrs, which coincides with the local protoplanetary disc having completely dispersed in these models.The left-hand panel shows the atmospheric mass fraction.Green lines with squares represent results from Owen & Wu (2016), for which boil-off was performed at 900K with initial atmospheric mass fractions of 0.3 and 0.1 in dot-dashed and dotted lines, respectively (note that initial atmospheric mass fractions in our models are calculated self-consistently through gaseous accretion).The right-hand panel shows planet sizes, including those that have been stripped of their atmosphere (atmospheric mass fraction < 10 −4 ) and sit on an Earth-like composition line (grey dashed).This panel can be interpreted as a primordial mass-radius distribution.
temperature.Assuming that our newly formed planet sits at the midplane of the disc, there are three primary heat sources to consider: viscous heating due to stellar accretion, irradiation from dust in the upper disc atmosphere and direct stellar irradiation once the disc is optically thin through its midplane.Firstly, the flux released at the disc midplane due to active accretion (e.g.Cassen 1994;Pringle 1981) is: where the factor of disc vertical optical depth  = Σ g /2 accounts for the flux at the midplane.The stellar accretion rate  * is related to the gas surface density Σ g via: 5 where  is the kinematic viscosity, parameterised by the standard dimensionless -viscosity parameter (Shakura & Sunyaev 1973): where we assume  = 10 −2 , which is appropriate for the inner disc where the viscosity is likely driven by the magneto-rotational 5 Typically, one sees an additional factor of 1 − √︃  *  in this expression and in Equation 29, which would arise from the assumption of a zero-torque inner boundary condition.Whilst appropriate for accretion discs orbiting black holes, the inner disc edge is truncated by the star's magnetic field in protoplanetary discs.
instability (e.g.Jankovic et al. 2021;Delage et al. 2023).Secondly, we follow the 2-layer model of Chiang et al. (2001) for heating at the disc midplane due to passive irradiation from dust in the upper disc.A radiative balance is found in which: where  * and  * are the stellar effective temperature and radius respectively,  s is the opacity of the disc interior to irradiation from the disc surface,  i is the opacity of the disc surface to irradiation from the disc interior,  i is the disc interior temperature,  = 0.5 is the fraction of the stellar hemisphere that is seen by a dust grain and  is the grazing angle due to irradiation incident on a flaring protoplanetary disc: where  irr is the vertical height at which irradiation is absorbed in the disc.For simplicity, we assume  irr = 4 (Chiang & Goldreich 1997), where  is the disc vertical scale height.In hydrostatic equilibrium: where  g ≡   *  H / B  * .We follow the numerical prescription from Chiang et al. (2001) to solve Equations 32, 33 and 34 for  irr and  i simultaneously on a grid of semi-major axis.For consistency with the rest of our dust-free gas calculations in this work, we adopt values of  s = 0.1 cm 2 g −1 and  i = 0.001 cm 2 g −1 , although we confirm that our results are insensitive to our choices in disc opacities.Instead, temperatures are predominantly controlled by the host stellar radius and temperature.Although minor, the flux due to this source of dust-driven irradiation is given by: The final source of heating is direct stellar irradiation once the disc is optically thin through its midplane, with flux given by: where  * is the stellar luminosity and  mid is the optical depth through the midplane of the disc: where we assume the disc's inner edge is truncated at  in = 10 * and the density at the disc midplane  mid () is related to the disc surface density via Equation 8. To account for the pre-main-sequence (PMS) evolution of stellar luminosity  * , effective temperature  * and radius  * , we incorporate the MIST evolution tracks for a 1 ⊙ star from Dotter (2016); Choi et al. (2016).The equilibrium temperature  eq as a result of all three heating terms is thus given by: where  visc ,  dust ,  * are given by Equations 29, 35 and 36 respectively.To determine the evolution of the disc surface density and planet equilibrium temperature, we take the models of XUV-driven disc dispersal from Owen et al. (2011).This dispersal mechanism relies on the trade-off between gas accretion onto the central star and gas being driven away vertically via photoevaporative winds generated by the high-energy flux from the star.As disc mass and accretion rate drop towards the end of the disc lifetime, there will come a point where gas is more likely to be removed by the disc wind than accreted.As a result, a gap is opened at this radius, typically ∼ 1AU, which disconnects an inner and outer disc -causing the two to evolve separately.Since accretion discs evolve on a viscous timescale which scales linearly with disc radius, 6 the inner disc rapidly drains onto the host star on timescales ∼ 10 5 yrs, which triggers boil-off.
In order to evolve the disc surface density, we assume an initial stellar accretion rate of 5 × 10 −8  ⊙ yr −1 and use this to calculate an initial gas surface density profile Σ( = 0, ) from Equation 30.We then normalise the models of disc evolution and dispersal from Owen et al. (2011) to have initial values of Σ( = 0, ), which then provide Σ(, ) for all times and stellar separations.The resulting disc surface densities are used to calculate a planet's equilibrium temperature from Equation 38.As a demonstration of this disc model, we calculate evolution for a planet at an orbital separation of 0.1 AU and set the MESA outer boundary temperature  BC (see Section 2.2) to these values as a function of time.
In Figure 8, we plot the disc surface density, atmospheric mass fraction, and equilibrium temperature with contributions from viscous, dust and stellar heating in dotted, dashed and dot-dashed lines respectively.We show this as a function of time, as opposed to time since disc dispersal as in Figures 4 -6 since the disc dispersal starts gradually and without a distinct onset in these more realistic disc evolution models.As in our previous simpler disc models from Sections 3.1 and 3.2, the planet accretes gas and then undergoes a rapid massloss phase as the disc disperses.One can see that the equilibrium 6 For a constant , the kinematic viscosity  ∼  and hence the viscous timescale  visc ∼  2  scales linearly with radius. .Evolution of disc surface density, atmospheric mass fraction, equilibrium temperature and thermal mass for a 3 ⊕ planet at 0.1 AU.In the third panel, we show the contributions of planet equilibrium temperature from viscous heating (dotted), passive heating from dust in the upper protoplanetary disc (dashed), and direct stellar heating which dominates once the disc is optically thin through its midplane (dot-dashed).In the bottom panel, the thermal mass represents the ratio of the planet's Hill radius to the disc's vertical scale height.If this ratio rises above unity, it indicates that a gap may open in the protoplanetary disc.
temperature initially falls at ∼ 2 Myr since the disc is optically thick and the surface density, accretion rate and therefore viscous heating are all dropping.In addition, PMS stellar evolution at this stage also contributes to a decreasing level of stellar irradiation, which reduces heating due to dust.The disc becomes optically thin at ∼ 4 Myr, revealing the planet to direct stellar irradiation for the first time.The equilibrium temperature rapidly rises and is dominated by the stellar luminosity.
In the bottom panel of Figure 8, we show the evolution of the planet's dimensionless thermal mass  th (e.g.Lin & Papaloizou 1985), defined as the ratio between a planet's Hill radius and its nascent disc's vertical scale height: The importance of thermal mass relates to the ability of a planet to open a gap in the protoplanetary disc.If  th > 1, then the Hill sphere extends beyond the disc scale height and implies sufficient local gravitational influence that gap formation might occur.Figure 8 implies that gap opening could occur, even for a small 3 ⊕ planet.
The reason is due to the strong dependence of thermal mass on disc temperature, which drops considerably before the host star directly irradiates the planet.Rapid gap opening may induce enhanced boiloff, as atmospheric mass is lost to the void opened in the planet's vicinity.We discuss this in Section 4.6.

DISCUSSION
We have presented new evolutionary models of small, close-in exoplanets undergoing gaseous accretion, boil-off and core-powered mass-loss.They highlight the importance of disc dispersal in setting the initial conditions for planetary evolution in the post-disc epoch.
They also demonstrate the need for self-consistent modelling, from gas accretion to direct irradiation from the host star.We have shown that planets will typically accrete a few to tens of per cent in atmospheric mass before boil-off begins, although we highlight that the exact values depend on parameters such as core mass, orbital separation, disc evolution, gas opacities and the time at which gaseous accretion began.Furthermore, once the disc begins to disperse, boiloff is triggered and removes ≳ 90% of atmospheric mass, again with exact values depending on core mass, orbital separation and gas opacities.Crucially, the ratio of cooling timescale to disc dispersal timescale is imprinted on the final atmospheric mass fraction of the planet, with quicker dispersal times inducing a more significant change in atmospheric mass for a given planet's cooling timescale.
In addition, we have considered the role of core luminosity and the transition from boil-off, to core-powered mass-loss and XUV-driven photoevaporation.We now discuss the implications of these findings.

A Consistent Evolution Sequence
In Jankovic et al. (2019), planet formation in the inner protoplanetary disc was studied for small-mass exoplanets.They modelled dust evolution, core accretion and then XUV-driven photoevaporation of planet atmospheres following disc dispersal, but excluding boil-off.
The final mass and radii of these planets were shown to be inconsistent with current observational data, with the models producing a higher atmospheric mass fraction than is observed.Then, in Rogers & Owen (2021) and Rogers et al. (2023a), XUV photoevaporation was exploited to rewind the clock for Kepler planet evolution.To be consistent with the present-day demographics, they showed that the initial atmospheric mass fraction of planets, once their discs had dispersed, must scale positively with core mass, consistent with the analytic findings of Ginzburg et al. (2016).When considering the entire exoplanet population with a peaked core mass distribution, this scaling results in typical sub-Neptune initial atmospheric mass fractions of ∼ 2% (see also Wu 2019) Ginzburg et al. 2016;Lee et al. 2018;Fung & Lee 2018), which is in tension with the value of ∼ 2% found above.Note that this discrepancy is likely not an indication that XUV photoevaporation underestimates mass-loss rates since planets beyond orbital periods of ∼ 30 days (which are not strongly affected by photoevaporation) have mass and radius measurements consistent with hosting a few per cent in H/He mass, instead of ≳ 10% (e.g.Weiss 2014;Jontof-Hutter et al. 2016;Wolfgang et al. 2016).As we have shown in this work, boil-off can provide a potential solution to this discrepancy via significant atmospheric escape during disc dispersal.We also highlight that, in contrast to some of the photoevaporation models discussed above, core-powered mass-loss models are typically initialised with the atmospheric mass fraction scaling of Ginzburg et al. (2016), thus implicitly assuming boil-off in their initial conditions (Ginzburg et al. 2018;Gupta & Schlichting 2019, 2020).As such, the aforementioned discrepancy was not seen in their work.
Whilst we have shown boil-off is an important process and may provide a solution to the aforementioned discrepancy, other mechanisms do exist that may also resolve this tension.For example, one important consideration that has not been included in this work is the effects of gas opacity on accretion (e.g.Ikoma et al. 2000;Piso et al. 2015).As additionally considered in multiple studies (Ikoma & Hori 2012;Lee et al. 2014;Lee & Chiang 2015, 2016;Lee et al. 2018;Fung & Lee 2018), dust acts to increase gas opacity in the planet's atmosphere and surrounding protoplanetary disc, which thus slows cooling and the rate of gas accretion.In this case, typical accreted atmospheres are ≲ 10%, thus reducing the significance of the discrepancy between accretion and mass loss.Additional ways to resolve tensions include: forming the planet at the very end of the disc lifetime (Lee & Chiang 2016); accretion-induced shocks (Lee 2019), increased planetary luminosity as a result of giant mergers (Liu et al. 2015;Inamdar & Schlichting 2016); or recycling of high-entropy material from the disc, although we highlight that the efficacy of this mechanism is currently under debate (Ormel et al. 2015;Fung et al. 2015;Cimerman et al. 2017;Kurokawa & Tanigawa 2018;Ali-Dib et al. 2020;Chen et al. 2020;Bailey & Zhu 2023;Savignac & Lee 2023).Indeed, we find that, in the lack of such mechanisms, many planets undergo runaway accretion before disc dispersal, particularly in the case of more realistic disc evolution models.Since we observe many more sub-Neptunes with atmospheric mass fractions ≲ 5% compared to Jovian planets, there must exist mechanisms such as those described that slow or delay the rate of accretion prior to disc dispersal.

A Gentle Breeze
Previous studies that considered mass-loss during disc dispersal have taken different approaches, yielding different results.In Ikoma & Hori (2012), accretion and mass-loss were considered in a hydrostatic equilibrium regime.They found that some planets will continue to lose their entire atmospheric masses as the disc disperses to maintain this equilibrium.Owen & Wu (2016), on the other hand, argued that planets will not be able to stay in hydrostatic equilibrium with the disc since dispersal occurs on timescales of 10 5 yrs (Kenyon & Hartmann 1995;Ercolano et al. 2011;Koepferl et al. 2013).In contrast, the cooling timescale of a planet has upper limits of  cool ∼ 10 6 − 10 7 yrs.In the case that hydrostatic equilibrium cannot be maintained, mass-loss will ensue via a hydrodynamic wind.Owen & Wu (2016) assumed the planet was initially fully convective out to the Bondi radius, with an escaping isothermal Parker wind (Parker 1958) and an instantaneous disc dispersal.Unlike Ikoma & Hori (2012), they found that planets do not lose all their mass, since mass-loss shuts off once the planet shrinks to small sizes ∼ 0.1 B .However, their assumptions implied that their calculated mass-loss rates were likely upper limits.
In this study, we have built upon the models of Ikoma & Hori (2012) and Owen & Wu (2016), by self-consistently modelling the planet through gas accretion and disc dispersal.Notably, the planet is allowed to remain in causal contact with the disc via an isothermal outflow and non-hydrostatic equilibrium if this is not physically possible.When introducing this self-consistent approach, we find that the planets require physics from both Ikoma & Hori (2012) and Owen & Wu (2016).Planets begin by cooling and accreting, before losing mass via subsonic breeze solutions, allowing the planet to sustain causal contact with the disc as it disperses.As this continues, the Mach number of the outflow must continually increase such that the escaping gas eventually becomes transonic (see Figure 4).We find that once this transition occurs, the planet leaves causal contact with the disc and, similar to Owen & Wu (2016), the mass-loss then quickly shuts off once the planet has shrunk to ∼ 0.1 B .We highlight that our models suggest that the vast majority of mass-loss occurs while the disc survives.Therefore, boil-off/core-powered mass-loss has all but concluded by the time the disc has dispersed in most cases.As a result, boil-off is likely unaffected by the dynamical instabilities highlighted in Wang & Lin (2023) due to mass-loss in multi-planet systems.During periods of high mass-loss rate, the disc is present and likely dissipates any orbital instabilities that arise due to changes in angular momentum.Although, we highlight that more work should be done to investigate this.
We find that the majority of boil-off cooling occurs due to advection (as initially shown in Owen & Wu 2016) and enhanced radiative cooling as atmospheric mass is rapidly lost, inducing steep temperature gradients due to d expansion.We additionally corroborate the findings of Owen & Wu (2016) in that the typical cooling timescales for planets that have undergone a boil-off are  cool ∼ 10 8 yrs at 10 Myrs (see Figures 3 and 5).Thermodynamically speaking, planets appear older than their age would suggest due to significant premature cooling caused by mass-loss during boil-off.As shown in Owen (2020), this means that observations of young planets with measured masses, radii and ages may allow us to constrain their initial entropies and, therefore, determine whether they did or did not undergo a boil-off phase.

Boil-off sets the initial conditions for planetary evolution
One of the key findings from this study has been the impact of disc dispersal timescale on atmospheric evolution.For a given planet mass and equilibrium temperature, the final state of a planet after disc dispersal is predominantly controlled by the ratio of a planet's cooling timescale to disc dispersal timescale.As we have shown in Figure 4, larger ratios of cooling timescale to disc dispersal timescale  cool / disc induce greater levels of advection and radiative cooling due to d expansion.This then causes planet evolution to approximately converge, independently of initial conditions.Nevertheless, very large ratios lead to a small dispersion in the final atmospheric mass (see the middle-left panel of Figure 3).This is indeed the conclusion that Owen & Wu (2016) came to under their assumption of an instantaneous disc dispersal, suggesting that our models are consistent with theirs in the limit of  cool / disc → ∞.
In the more realistic case of a non-instantaneous disc dispersal based on observed timescales of ∼ 10 5 yrs (Kenyon & Hartmann 1995;Ercolano et al. 2011;Koepferl et al. 2013), our results suggest that the entropy, and therefore thermal structure of the planet, is set by boil-off.Information about the ratio of a planet's cooling timescale to disc dispersal timescale is imprinted on the atmospheric structure of the planet, which will affect how a planet evolves after the disc has dispersed.This opens up the possibility of using demographic surveys of exoplanets to perform physics-based inference with boiloff, core-powered mass-loss and XUV photoevaporation models to place constraints on disc dispersal timescales, which is left for future work.

Comparison with Ginzburg et al. (2016)
In Ginzburg et al. (2016), a purely analytic approach, based on radiative-convective equilibrium models of H/He atmospheres, was used to predict the atmospheric mass  ≡  atm / core at the end of the boil-off phase (referred to as 'spontaneous mass-loss' in their study).Unlike Owen & Wu (2016), their initial conditions were analytically calculated self-consistently with gaseous accretion prior to disc dispersal.Under the assumption of constant gas opacity, ideal gas equations of state and an adiabatic index of  = 7/5, they predict: As shown in Figure 7, we find a steeper dependence of atmospheric mass fraction with core mass after boil-off and an inverse dependence on equilibrium temperature.In other words, we find that hotter planets lose more atmosphere.This latter difference arises from our more realistic modelling of opacities from Freedman et al. (2008), which vary with temperature.Indeed, if one performs the analysis of Ginzburg et al. (2016) with temperature-dependant opacities, an inverse dependence of atmospheric mass fraction with equilibrium temperature can be recovered.We speculate that the steeper massscaling found in this work compared to Ginzburg et al. (2016) and inferred in Rogers et al. (2023a) is likely due to the fact that we are unable to model the gas-accretion phase on a population level over enough cooling timescales before the onset of disc dispersal.We leave the details of understanding these different mass scalings for future work.

Evolution beyond boil-off
Once the initially rapid mass-loss, cooling and contraction involved in the boil-off phase have ceased, the atmospheric evolution of small, close-in planets will be dictated by atmospheric processes such as core-powered mass-loss (e.g.Ginzburg et al. 2018;Gupta & Schlichting 2019) and/or XUV-photoevaporation (e.g.Lopez & Fortney 2013;Owen & Wu 2013).In Section 3.2, we considered the effect of core luminosity (arising from core cooling and radiogenic heating) and its influence on atmospheric evolution.As boil-off progresses and the core cooling luminosity increases, it will eventually dominate the planet's energy budget.As such, the core provides the energy required to lift gas to the base of the bolometrically heated outflow.This changeover offers a convenient definition for the transition between boil-off and core-powered mass-loss, as highlighted in Figure 5.We have demonstrated in Figure 6 that larger atmospheres around larger planet cores release more energy and transition to corepowered mass-loss later.Nevertheless, the bottleneck to boil-off or core-powered mass-loss is the speed at which bolometric heating can remove material via a Parker outflow.As such, the final atmospheric mass fractions with and without core luminosities are very similar, as shown in Figure 5.As described in Owen & Schlichting (2023), XUV photons from the host star will eventually penetrate beneath the sonic surface of the planet at  B .In this case, XUV photons can heat and ionise the hydrogen-dominated gas before it reaches transonic velocities and becomes causally disconnected from the planet's interior.As such, a strong XUV-driven photoevaporative wind is launched from beneath the Bondi surface.This defines the transition from boil-off or corepowered mass-loss to XUV photoevaporation, which we highlight in Figure 6.One can see that some smaller-mass planets are affected by all three evolution processes: boil-off, core-powered mass-loss and photoevaporation.This agrees with the speculated evolutionary scenario proposed by Owen & Schlichting (2023).The smallest, highly irradiated planets are stripped of their atmospheres during disc dispersal prior to XUV heating.We note that throughout this work, we have defined a stripped planet as one with a negligible atmospheric mass fraction  < 10 −4 .As shown in Misener & Schlichting (2021), however, small-mass primordial atmospheres ( ≲ 10 −5 ) can be retained in the scenario of core-powered mass-loss.Nevertheless, such small atmospheres are unlikely to be retained in the XUV photoevaporation regime unless one gets significant mass-dependent fractionation.
Our models suggest that many small-mass planets may have been stripped of their primordial atmospheres through a combination of boil-off and core-powered mass-loss.This creates a radius valley prior to the effects of photoevaporation, which typically acts on ∼ 100 Myr to Gyr timescales.Note, however, that this is distinct from the 'primordial' radius valley as predicted Lee & Connors (2021); Lee et al. (2022) which arises purely due to gas accretion and would therefore be present even before boil-off / core-powered mass-loss.Regardless, the combination of these proposed mechanisms suggests that the radius valley should be observable in the demographics of very young planetary systems.

A Cold Sunrise
In Section 3.3, we considered a more realistic disc evolution model that included an evolving disc midplane temperature.As shown in Figure 8, the equilibrium temperature experienced by a close-in planet drops, often considerably, below its future value once the disc is optically thin and the planet is directly irradiated by its star for the first time.The effect of this is two-fold: firstly, cooler temperatures in a dust-free environment increase gaseous accretion rates since planets can cool and contract more efficiently to accrete more material within their Bondi sphere.Secondly, cooler temperatures reduce the sound speed of the escaping gas, which reduces the mass-loss rate as the disc disperses.Nevertheless, we still find that planets accrete ≳ 10% in H/He dominated material and then boil-off during disc dispersal to atmospheric mass fractions of order ∼ 1%.As discussed in Section 4.1, we also find that many planets with larger core masses undergo runaway gaseous accretion due to the lower temperatures.However, since we observe many planets of mass ≳ 5 ⊕ with H/He atmospheric mass fractions ≲ 10%, there must exist mechanisms by which gaseous accretion is reduced or delayed to be consistent with results (e.g. Lee 2019).In Figure 8, we also evaluated the planet's thermal mass to show that even a small 3 ⊕ planet may open a gap in the protoplanetary disc towards the end of disc dispersal.Under such circumstances, accretion can only occur through the gap, possibly via a circumplanetary disc, yielding lower mass fluxes.This latter point raises the question of whether gap-opening, akin to the local gas disc dispersing on very quick timescales, can induce some form of mass-loss?Whilst we leave this for future work, one can consider larger-mass planets that have already opened gaps prior to boil-off.The act of opening a gap may induce mass-loss as the local gas density rapidly drops.On the other hand, when the disc does disperse, as considered in this work, the planet will not experience any further boil-off since the planet has already been exposed to a vacuum in its local vicinity during the gap-opening process.Indeed, Rogers et al. (2023b) found tentative evidence for the lack of boil-off in evolved planets with measured mass and radii around M-dwarfs.An inferred increase in atmospheric mass fraction dispersion for planets ≳ 10 ⊕ suggested that they did not follow the physics of boil-off, leading to far greater initial atmospheric masses ∼ 10% after disc dispersal.This could be caused by slow disc dispersal rates, or gap-opening as discussed above.

Model Improvements
We have shown that self-consistent modelling of gas accretion and mass-loss can shed light on the atmospheric evolution of exoplanets during disc dispersal.However, there is one major assumption that we have made that limits this analysis: an isothermal outflow.As shown in Figure 1, the assumption of isothermality is likely valid in the radiative zone of a planet.However, this may not be the case, particularly in optically thin regions far from the planet.In assuming an isothermal outflow, we inherently assume that energy is inserted into the flow to keep it heated to  eq and thus remove it from the planet's Bondi/Hill sphere.In reality, there may be circumstances where the flow cannot be heated to sufficient temperatures to leave the planet's influence.As discussed in Section 4.2, mass-loss in the current model shuts off due to the planet cooling and contracting.However, this effect may be compounded when considering nonisothermal outflows.One could use hydrodynamic models that consider the radiative transfer and subsequent mass-loss in the upper atmosphere to understand this effect.These models can then be used to construct a grid of mass-loss rates, which can be implemented in MESA.This was indeed the approach of Owen & Wu (2013), applied to XUV photoevaporation, in using the mass-loss rates from hydrodynamic models from Owen & Jackson (2012).This is left for future work.
Finally, a significant effect that has been neglected is that of dust.Thus far, only dust-free gas has been considered by way of the opacity tables from Freedman et al. (2008) utilised in atmospheric structure, and our choices in disc opacities.In reality, dust is also present in the protoplanetary disc and, subsequently, the planet's atmosphere.Since opacity changes the rate of cooling and therefore the rate of accretion and atmospheric escape, it will likely alter the exact amount of mass that is accreted onto, or lost from, sub-Neptune atmospheres during disc dispersal (e.g.Owen & Wu 2016).Investigating this effect is also left for future work.

CONCLUSION
We have considered the atmospheric evolution of small, close-in exoplanets hosting H/He-dominated atmospheres through protoplanetary disc dispersal.As the surrounding disc gas dissipates, a powerful hydrodynamic outflow is triggered from the planet, inducing rapid cooling and contraction of the atmosphere.We have shown that this process, referred to as 'boil-off', is instrumental in setting the initial conditions for other mass-loss mechanisms, such as core-powered mass-loss and XUV photoevaporation.Our main conclusions are as follows: • Planets ≲ 10 ⊕ will typically accrete a few to tens of per cent in atmospheric mass during the disc lifetime, depending on disc conditions.The onset of disc dispersal then triggers boil-off, in which planets typically lose ≳ 90% of their accreted atmospheres.This mechanism provides a possible solution to the documented discrepancy between predicted accreted masses and the initial conditions inferred from Kepler data with XUV photoevaporation, as shown in Jankovic et al. (2019); Rogers & Owen (2021).
• Unlike previous studies of mass-loss during disc dispersal (e.g.Ikoma & Hori 2012;Owen & Wu 2016), we show that boil-off occurs via subsonic breeze outflows, which allow a causal connection to exist between escaping atmosphere and disc.As such, boil-off has all but concluded by the time the disc has fully dispersed.This typically coincides with the planet launching a transonic wind, which yields relatively low mass-loss rates due to the planet's contracted size.
• For a given planet mass and equilibrium temperature, we find that the amount of mass lost during boil-off is directly dictated by the ratio of cooling timescale to disc dispersal timescale.From disc observations, one expects this ratio to be greater than unity; however, the larger the ratio, the more mass is lost.This leads to convergent evolution, largely independent of how much atmosphere was initially accreted.We show that this is due to advective cooling and a brief period of d expansion with subsequent radiative cooling that occurs at the onset of boil-off.
• The effect of core luminosity, with contributions from core cooling and radioactivity, can eventually dominate a planet's energy budget once atmospheric contraction has slowed, marking the transition to core-powered mass-loss (e.g.Ginzburg et al. 2016;Gupta & Schlichting 2019;Misener & Schlichting 2021).We have shown that a combination of boil-off and core-powered mass-loss is capable of stripping small mass (≲ 2 ⊕ ), close-in ( eq ≳ 1000 K) planets of their primordial atmospheres even before the disc has fully dispersed, without the need of late-time mass-loss processes such XUV photoevaporation.
• Eventually, XUV photons can penetrate below a planet's sonic surface and induce a more powerful photoevaporative hydrodynamic escape.We follow the work of Owen & Schlichting (2023) and corroborate that the radii of smaller mass planets are predominantly sculpted by core-powered mass-loss.In contrast, higher-mass planets are sculpted by photoevaporation.
• Finally, we considered a more realistic disc evolution model with evolving midplane temperature to show that planets may open gaps in their protoplanetary disc during boil-off.The gap may enhance massloss for low-mass planets due to rapid reduction in local confining pressure.For high mass planets ≳ 10 ⊕ , this may occur during its accretion phase and limit the amount of gas that can flow onto the planet.More work is needed to understand this process.
We highlight that the assumptions of an isothermal outflow and dust-free opacities have limited this study.In future work, we plan to incorporate such effects to understand the radiative processes that limit/enhance mass-loss during disc dispersal.
Photon Mean-Free Path Pressure Scale Height

Figure 3 .
Figure3.Accretion and boil-off are shown for planets with core masses of 5 ⊕ and constant equilibrium temperature of 900K.The disc begins to disperse at  disp = 3 Myrs, with disc gas surface density shown in the top panels for initial surface density of Σ 0,g = 3 × 10 4 g cm −2 and dispersal timescales of  disp = 0.01, 0.1 and 1.0 Myrs shown in left, middle and right-hand columns respectively.Atmospheric mass fraction is shown in the middle row, with planets forming at 1.0, 1.4 and 2.0 Myrs shown in blue, orange and green respectively.The cooling timescale is shown in the bottom row, with its ratio to the disc dispersal timescale recorded at 3 Myrs.

Figure 4
Figure 4 also shows the entropy at the base of the atmosphere as a function of time after disc dispersal begins.Planets that undergo a more rapid disc dispersal process end with a lower final entropy, implying they have cooled more dramatically.To understand this, we also show the luminosity at the radiative-convective boundary (RCB)  rcb as a function of time.Planets experiencing a quicker

Figure 4 .
Figure 4. Boil-off is shown for planets with core masses of 5 ⊕ and constant equilibrium temperature of 900K.Planets all begin evolution at 2 Myrs, corresponding to green lines in Figure3.Time is measured from the point at which disc dispersal began, in this case at  disp = 3 × 10 6 yrs.Upper-right shows the atmospheric mass fraction evolution, as in Figure3, additionally with the fractional change / init .The lower-right panel shows the Mach number at the Bondi radius of the outflow, whilst the bottom-left shows the entropy at the base of the atmosphere.The upper-left panel shows the luminosity at the radiative-convective boundary.Solid, dashed and dot-dashed lines represent disc dispersal timescales of  disp = 10 4 , 10 5 and 10 6 yrs respectively.

Figure 5 .
PhotosphereRadiative-convective boundary Boil-off → Core-powered mass-loss Figure8.Evolution of disc surface density, atmospheric mass fraction, equilibrium temperature and thermal mass for a 3 ⊕ planet at 0.1 AU.In the third panel, we show the contributions of planet equilibrium temperature from viscous heating (dotted), passive heating from dust in the upper protoplanetary disc (dashed), and direct stellar heating which dominates once the disc is optically thin through its midplane (dot-dashed).In the bottom panel, the thermal mass represents the ratio of the planet's Hill radius to the disc's vertical scale height.If this ratio rises above unity, it indicates that a gap may open in the protoplanetary disc.
. On the other hand, many studies considering accretion under a range of conditions have shown that in the dust-free case, as considered in this work, planets ≳ 5 ⊕ can accrete ≳ 10% in atmospheric mass by the end of the disc lifetime (e.g.Ikoma & Hori 2012; Lee et al. 2014; Lee & Chiang 2015;