The evolution of catastrophically evaporating rocky planets

In this work, we develop a rocky planet interior model and use it to investigate the evolution of catastrophically evaporating rocky exoplanets. These planets, detected through the dust tails produced by evaporative outflows from their molten surfaces, can be entirely destroyed in a fraction of their host star's lifetime. To allow for the major decrease in mass, our interior model can simultaneously calculate the evolution of the pressure and density structure of a planet alongside its thermal evolution, which includes the effects of conduction, convection and partial melting. We first use this model to show that the underlying planets are likely to be almost entirely solid. This means that the dusty tails are made up of material sampled only from a thin dayside lava pool. If one wishes to infer the bulk compositions of rocky exoplanets from their dust tails, it is important to take the localised origin of this material into account. Secondly, by considering how frequently one should be able to detect mass loss from these systems, we investigate the occurrence of sub-Earth mass exoplanets, which is difficult with conventional planet detection surveys. We predict that, depending on model assumptions, the number of progenitors of the catastrophically evaporating planets is either in line with, or higher than, the observed population of close-in (substellar temperatures around 2200 K) terrestrial exoplanets.


INTRODUCTION
Understanding rocky planet interiors is an essential part of exoplanet science.Firstly, their compositions lend insight into planet formation.Secondly, the interior chemistry determines the atmospheric chemistry through outgassing for planets below a few Earth masses that do not possess a H 2 /He envelope.Atmospheric chemistry is influenced both through exchange with the early molten surface (e.g., Abe & Matsui 1986;Gaillard & Scaillet 2014;Lichtenberg et al. 2021), and across the lifetime of the planet through volcanism (e.g., Kite et al. 2009;Noack et al. 2014;Tosi et al. 2017).Furthermore, the thermal evolution of the interior affects the efficiency of outgassing processes and what species get outgassed, as well as providing thermal input into the base of the atmosphere.Thus, the interior is important to explain the range of possible exoplanetary atmospheres, including habitability.Chemical compositions of exoplanet interiors are most directly observed in the atmospheres of polluted white dwarfs (Zuckerman et al. 2010), which show metal lines due to planetary bodies that have recently been accreted and have been shown to have compositions that are consistent with Solar System rock (e.g., Hollands et al. 2017;Harrison et al. 2021).However, as these systems can only track the bulk composition of the pollutant, little information about the compositional structure or influence on the planet's atmosphere can be inferred.Another avenue of interest is the study of plan-ets' atmospheres, particularly those of hot planets where the atmospheric composition is likely dominated by interaction with the surface and interior (e.g., Zilinskas et al. 2022).However, the atmospheric interactions are still not fully understood and are complicated by atmospheric processes such as loss and photochemistry, even if such tenuous atmospheres become observable.The catastrophically evaporating planets, however, are less dependent on such atmospheric processes because they are observed directly through solid material.
Three catastrophically evaporating planets -Kepler 1520b, (formerly KIC1255b, Rappaport et al. 2012), KOI-2700b (Rappaport et al. 2014) and K2-22b (Sanchis-Ojeda et al. 2015) -were discovered by the Kepler/K2 missions (Borucki 2016;Howell et al. 2014).The distinguishing features of the lightcurves of these systems are (i) all three systems have highly variable transit depths indicating that the orbiting body is not a single solid object.(ii) there is an increase in flux directly before transit, which can be explained by the forward scattering of starlight by dust (iii) Kepler 1520b and KOI-2700b show a highly asymmetric averaged transit due to an extended tail (see e.g., van Lieshout et al. 2016, Figure 1).These features are well explained by the transit curves being due to tails of dusty material (van Lieshout & Rappaport 2018).
The short periods of the dust tails imply high surface temperatures for any planet at that orbital distance (see Table 3), sufficient to melt rocky materials.Therefore the origin of the dusty tails is thought to be the evaporation of rocky material from the molten surface of an underlying planet, which then expands and condenses as dust further from the surface.This idea was originally proposed by Rappaport et al. (2012), with further physical modelling by Perez-Becker & Chiang (2013).More recent models of the outflows have focused on photochemistry (Ito & Ikoma 2021), day to nightside flows (Kang et al. 2021, see also Castan & Menou 2011) and variability (Bromley & Chiang 2023).Booth et al. (2023) were notably able to show that dust can form self-consistently during this atmospheric escape under the right conditions.
Since the material in the dusty tails comes directly from the solid portion of the planet, these systems are particularly interesting in the wider context of exoplanets because they may allow the study of interior composition.Thus far, attempts to constrain the composition of tail material have depended on modelling the light curves.For instance, van Lieshout et al. (2014van Lieshout et al. ( , 2016) ) found that the lightcurve of Kepler-1520b is consistent with corundum (Al 3 O 2 ) and KOI-2700b with corundum or fayalite (Fe 2 SiO 4 ), but both are inconsistent with pure iron or carbon, based on the species' sublimation rates.More recent modelling efforts (Campos Estrada et al. 2023) suggest that the tails are best explained by iron-rich silicate dust.In addition, the size of the forward scattering peak can be used to infer grain sizes, which are required to break degeneracies with composition in modelling, and are constrained to a range of 0.1 -1 µm (Budaj 2013).
Further, multi-wavelength observations will allow inferences to become more detailed.JWST infra-red spectra may allow silicate minerals to be identified through their resonant features around 10µm (Bodman et al. 2018).It may also be possible to detect atomic lines from the gas in the tails (van Lieshout & Rappaport 2018), although attempts to find atomic species in the tail of K2-22b have so far been inconclusive (Ridden-Harper et al. 2019).
As has been noted in previous works (e.g., Perez-Becker & Chiang 2013), and as we will also show, the high mass loss rates inferred for these planets means that they can be entirely destroyed within several Gyrs.Models of the mass loss process agree that the mass loss rates should generally increase with decreasing mass, leading to an accelerating mass loss before total destruction.A consequence of this is that the observed systems must be within a region, of temperature, mass and time, where their mass loss rates are high enough to be observed, but the planet has not yet been destroyed.Therefore, the fact that any are observed at all may tell us about the occurrence rate of these planets' progenitors.This is of particular significance because,as we will show (see also Perez-Becker & Chiang 2013), these planets must have low current and initial masses (≲ 0.3M⊕).Therefore, they are part of a distribution not observed by conventional detection techniques (e.g., Christiansen et al. 2014).
In order to fully understand the catastrophically evaporating planets, both in their own right and to make full use of their potential as probes of interior physics, it is necessary to have a model of their interiors.The evolution of the underlying rocky planet may affect the thermal state at the bottom of the mass outflow.Additionally, and more significantly, the interior evolution sets the surface composition, which can be probed observationally through the dusty tails.Different compositions, for instance, through different depths that melting reaches, can vastly affect the chemicals that would be present in the outflows (e.g., Schaefer & Fegley 2009).
In this work, we present our evolutionary model for the interior of catastrophically evaporating planets.In particular, we investigate the evolution of melt, since the molten regions convect over much shorter timescales than solid regions, and so have an important influence on composition by circulating chemicals to the surface.We then use this model to investigate the population of the progenitors of the catastrophically evaporating planets.
In §2, we describe our 1D numerical model for rocky interiors and its specific application to these systems.In §3, we present the results of our thermal evolution calculations, with a discussion of their consequences and limitations in §4.
In §5, we investigate the implications for the the occurrence rate of low-mass planets.We conclude with a summary in §6.

INTERIOR MODEL
We have developed a 1D code for modelling the evolution of a rocky planet with a time-varying mass, which we shall apply to the evolution of catastrophically evaporating planets.The numerical scheme is based on stellar structure codes, as described in Bodenheimer et al. (2006) and Kippenhahn et al. (2012).We solve for the structure of the planet's rocky mantle and iron core.Here we summarise the physics included and its basic operation.

Basic Equations
The essential equations for the internal structure of a spherical body, written with mass, m, as the independent variable, are: (i) Hydrostatic Equilibrium dP dm = − Gm 4πr 4 (1) where P is the pressure, r the radius at the mass point m and G the gravitational constant.
(ii) Mass conservation where ρ is the density.
(iii) Energy conservation whereL is the luminosity, and H is the heat generation rate per unit mass (through radioactive decay in the case of our planets).t is time, T is the temperature and S is the specific entropy, meaning T dS/dt is the rate of heat exchange from a unit mass.The second equality follows from thermodynamic relations.CP is the specific heat capacity at constant pressure, and is a measure of thermal expansivity.
One also needs an equation for the temperature gradient, which will be related to heat transport.Following the conventions in stellar interiors, we define the temperature gradient ∇ ≡ d ln T d ln P (5) and use hydrostatic equilibrium (Equation 1) to write: We will show in §2.3 how we determine ∇ from consideration of the heat transport.Equations 1-3 and Equation 6 form a closed system of differential equations when combined with an equation of state and other material properties, which may be functions of temperature and pressure ( §2.4).The only other information required to solve for the four dependent variables (P , r, L and T ) as functions of m is four boundary conditions.The two inner boundary conditions are simply R = 0 and L = 0 at m = 0.At the outer boundary, we use a fixed P = P0 and a relation between the outer temperature T0 and the other outer properties.We explain these outer boundary conditions for the cases we investigate here in §2.7.
In order to solve this system of equations, we follow a Henyey scheme, the details of which can be found in Bodenheimer et al. (2006), chapter 5, for example.In summary, we set up a mass grid such that the differential equations become a series of simultaneous equations.We solve these through Newton's method, up to a given tolerance threshold.The tolerance criterion we use is that the difference in a dependent variable across a mass grid cell must have a fractional accuracy of 10 −6 or better.The time dependence in the energy equation is solved implicitly.

Melting
Rocks are chemical mixtures and thus can be partially molten even in thermodynamic equilibrium at fixed temperature and pressure.This process is complicated because the melt composition will generally differ from that of the solid rock, and will also depend strongly on both the local conditions and the overall composition of the rock.To simplify the full problem, we introduce a parameter for the mass fraction of melt, ϕ, which is simply a function of P and T .
Following other works (e.g., Abe 1993), we use the simple linear function where T sol is the solidus, the temperature below which the material is completely solid (i.e., ϕ=0 for T <T sol ), and T liq is the liquidus, the temperature above which the material is completely liquid (i.e., ϕ=1 for T > T liq ); both T sol and T liq are functions of P .Following e.g., Abe (1993); Bower et al. (2018), we consider these fixed functions, but again, the reality is far more complex and highly composition-dependent.
For the mantle solidus and liquidus profiles, we use the fit to the Simon and Glatzel equation (Simon & Glatzel 1929) from Andrault et al. (2011) for high pressures, and the curves in Litasov & Ohtani (2002) for low pressures.For use in our code, we tabulate these functions and access values using cubic Hermite interpolation.

Heat flow in the mantle
We include heat flow in our model by finding how the temperature gradient ∇ depends on the energy flux.In practice, we find how the energy flux, F = L/(4πr 2 ), depends on the temperature gradient and invert this function numerically.1

Conduction and Convection
Heat transport occurs in rocky interiors through conduction and convection.We model conduction using Fick's law, so conductive heat flux is given by: assuming a constant conductivity, k.
Convection is extremely important for the evolution of rocky planets, both in the liquid and solid phases.To model convection we use mixing length theory (see e.g.Abe 1995;Kippenhahn et al. 2012) and use the equation: where l is the mixing length, u is the speed of convection and ∇ Ad is the adiabatic, logarithmic temperature-pressure gradient (see Equation 5).The total flux is given by The mixing length prescription requires an estimate of the convective velocity u.To find the velocity, we consider the forces Fi acting on a parcel of the material moving due to convection.If the parcel is moving at a constant speed, the buoyancy force must be balanced by any drag forces.The drag forces are ram pressure, which is most important in the low viscosity limit, and viscous drag, which is more important in the high viscosity limit.These three forces may be given by the following formulae: where g is the gravitational acceleration, ν ≡ η/ρ is the kinematic viscosity, and η is the dynamic viscosity.Here R, A and V are the fluid parcel's radius, cross-sectional area and volume, respectively.Combining these results in a quadratic equation for u, the (positive) solution to which is This has the limits for high and low viscosity respectively.
We then choose geometric factors relating R, A and V to l, l 2 and l 3 in order to produce the velocity limits used by Abe (1995) for high and low viscosity.These are which it can be seen are equivalent to our Equation 13 up to numerical factors.
Our formulation using Equation 12 makes the transition smooth, which aids numerical convergence, as opposed to a switch at a critical value used in other works (e.g., Abe 1995;Bower et al. 2018).

Application in partial melt
A general formula for the adiabatic gradient in a medium is Under the assumption that the melt fraction ϕ is always in equilibrium, the adiabatic gradient is given by the "wet adiabat", where the latent heat alters the values of δ/ρ and CP , and thus changes the value of ∇ Ad (Equation 15).These properties are also adjusted throughout the other equations.
The adjustments to the density ρ, thermal expansivity δ and heat capacity CP may be easily derived by assuming additive volumes and entropies of the melt and solid phases, i.e., V = V l ϕ + Vs(1 − ϕ) and S = S l ϕ + Ss(1 − ϕ), where V and S are specific volume and entropy and subscripts l and s denote liquid and solid properties, respectively.Consequently, the density is given by To find δ and CP under partial melting one only requires the definitions The results are where ∆V ≡ 1 ρ l − 1 ρs is the specific volume change of melting, and ∆S is the specific entropy change of melting.
Consequently, Equation 15becomes For this work, we assume that the timescales of melting are shorter than those of mixing, so melt equilibrium is maintained, and our mixing length formulation fully captures the heat flow in the system.

Physical properties of the mantle
In order to solve the equations in the previous three sections we need to know the values of the relevant physical properties, density, heat capacity etc., for solid and molten rocky materials.In this subsection, we describe how we calculate them.

Equations of state
ρ, CP and δ are all obtained as functions of P and T .For the solid phase, we use the equation of state for enstatite, the Earth's main mantle component, from Stixrude & Lithgow-Bertelloni (2011).For the melt, we use the MgSiO 3 equation of state in Wolf & Bower (2018).
We precalculate these properties on a linear grid of temperature and pressure and retrieve values using bicubic Hermite interpolation.Under partial melting, density, thermal expansivity and heat capacity are given as a combination of the melt and solid properties, as described above in §2.3.2.

Volume and entropy change for melting
In order to describe melting, it is necessary to know the change in specific volume, ∆V , and entropy, ∆S, between the solid and melt (see Equation 17-19).∆V can be calculated from the densities given by the two equations of state.However, calculating ∆S from an equation of state requires an absolute entropy scale, which does not emerge naturally from our prescriptions.Instead, we calculate the entropy change through which is derived in Appendix A. We use the estimate: which simply interpolates between the T − P gradient of the ϕ = 0 line, the solidus, and the ϕ = 1 line, the liquidus.

Viscosity
The effective viscosity varies by many orders of magnitude as the melt fraction changes.A steep change from more liquid to solid-like behaviour is often taken to occur at a critical melt fraction, ϕc (e.g.Solomatov 2007).This melt fraction is essentially the fraction of melt at the point when closely packed spheres of the typical size of crystals can no longer move past each other.We use an experimentally measured value of ϕc = 0.4 from Lejeune & Richet (1995).
Additionally, the viscosity of the solid component is temperature and pressure-dependent, and we model it using an Arrhenius law for diffusion creep (Tackley et al. 2013) where the s subscript once again denotes solid, and the values of parameters are shown in Table 1.
For low melt fractions (0 < ϕ < ϕc), we use an exponential parameterisation for the dynamic viscosity (Kelemen et al. 1997) with ηs(P, T ) as given in Equation 22and evaluated at the solidus, i.e., just before partial melting occurs.We assume that diffusion creep is the dominant mechanism of solid deformation (e.g., Tackley et al. 2013), and so use αη = 26, as found in Mei et al. (2002).
For high melt fractions (ϕ > ϕc), we use the formula from Roscoe (1952) where the l subscript denotes liquid.We take the liquid viscosity η l to be a constant, as its temperature dependence is small relative to the solid phase and to changes in melt fraction (Solomatov 2007).
We combine these formulae into a smooth function, with no singularity at ϕ = ϕc, given by

The iron core
Planets above a few thousand km in radii are likely to have formed an iron core through gravitational settling of the denser iron from the mantle (e.g., Elkins-Tanton 2012).We concentrate on planets with a core mass fraction of 0.3, similar to the Earth.Were planets to have significantly different initial silicate-to-iron ratios to the Earth, or undergo a major collision that might strip the mantle, then it is possible that the core mass fraction would be different for exoplanets.However, we find that differences in evolution are generally not significant enough to alter the conclusions we draw.We assume the core is pure iron and solve for its structure using our Henyey scheme, like the mantle.For iron, we use physical properties for the γ phase of iron from Dorogokupets et al. (2017), as it is the appropriate phase at the pressures in the cores of the low mass planets we will consider (e.g., Tsujino et al. 2013).).The inner part of the planet is considered spherically symmetric, and its structure is solved for by our Henyey scheme.Its edge is at a pressure P 0 , which has a spherically symmetric temperature, T 0 .The planet's surface has a day-to-nightside temperature gradient, Ts(θ), due to the stellar irradiation determined by Equations 26 and 27.For the outer layer (P < P 0 ), we consider radial fluxes F (θ) such that the flux integrated over the surface is equal to the luminosity for the interior, L.
For these calculations, we neglect the fact that at early times the mantle will be molten and set the temperature structure in the core to be adiabatic.The reason for these simplifications is that we are mainly interested in the amount of energy provided to the mantle, which is likely not affected by these assumptions, as will be discussed in §4.3.Modelling the core in full would add extra complexity, especially as its evolution is not well constrained (see e.g., Zhang & Rogers 2022).

Radioactive heating
For Gyr old planets, the energy generated by the decay of long-lived radioactive isotopes is an important contributor to the total luminosity.The significant elements in the Earth are 232 Th, 238 U, 40 K and 235 U (Turcotte et al. 2002).These elements are lithophiles, which means they preferentially dwell in the mantle.Over shorter period(less than a few Myrs), short-lived radioisotopes are more important; in the Solar System, the most significant of these were 26 Al and 60 Fe.Our models use the same concentrations of long-lived radioisotopes as found in the Earth's mantle, including their decay over time (data from Turcotte et al. 2002).We also include 60 Fe in the core and 26 Al in the mantle at the concentration estimated for the Solar System in Lugaro et al. (2018).

Outer boundary conditions
The catastrophically evaporating systems we are investigating are very close to their host stars and are likely tidally locked.Therefore, they have permanent daysides that are highly irradiated.We take this lack of symmetry into account in our evolutionary models by adapting the outer boundary conditions.The full problem of heat transport in the planet is complex due to the mixture of radial and angular heat transport and the effect of heat transport in any magma pool; we, therefore, make a series of assumptions as discussed below.We explain the principles of the approach here, with some additional details supplied in Appendix B.
The first and most important assumption is that there is some depth within the planet below which the planet can be assumed to be spherically symmetric, and thus the 1D structure model described in previous sections can be applied unchanged (see schematic in Figure 1).We designate the pressure this occurs at as P0.For regions at lower pressures, i.e., above this boundary, we allow the symmetry to be broken, and quantities then depend upon the angle θ from the substellar point (the point on the planet's surface closest to the star).This assumption is justified if the inward energy flux, due to external heating from the star, that penetrates deeper than P0 is much less than the outward flux from the planet cooling.We will discuss this further in §4.2, with reference to our results, but the essential reason this is true in most cases is that heat can only be transported inwards by conduction, not convection, and thus the efficiency of inward heat transport is low.
What we require is a relation between the temperature at the boundary of the spherically symmetric interior, T0, and the interior luminosity L (see Figure 1).For this, we must consider the region of the planet close to the surface (P < P0) where there is an angular dependence.In order to model this region, several assumptions must be made.
We begin by considering the surface (P = 0).We assume that there is no angular redistribution of the star's energy here.The first reason for this is that since the planets are tidally locked, there can be no redistribution due to rotation.Secondly, any atmosphere generated through outgassing volatiles will likely be thin (for instance, the atmospheric escape models of Booth et al. 2023 generate atmospheres with maximum pressures of ∼ 10 −5 bar) and so we expect them to be unable to transport heat efficiently.Thirdly, Kite et al. (2016) showed that a surface magma ocean cannot transport enough heat laterally to decrease the temperature gradient imposed by irradiation.Consequently, we treat the surface as a local black body; thus, its temperature is given by where Firr(θ) is the irradiation from the star, and F (θ, T0) is the angle-dependent heat flux from the interior, which depends on T0.)A priori, the value of both F (θ, T0) and Ts(θ) are unknown, and the aim is to find both together.Assuming the planet is far enough from the star, the stellar irradiation is plane-parallel and given by where R * , T * and a are the stellar radius, stellar effective temperature and planet's semi-major axis, respectively.For much of the planet's lifetime F < Firr and so the surface temperature is proportional to Firr(θ) 1 4 , yielding a steep temperature gradient from day to nightside.
Considering now the sub-surface part of the P < P0 region, we assume that the angular flux is much smaller than the radial flux in this region as well.This will be justified in §4.2.This means that rather than a full 2D problem, the heat transport within the P < P0 layer is simply a series of 1D heat transport equations at any angle θ (see Figure 1).In fact, at any given θ the conduction/convection equations (Equations 8-10) become just one ordinary differential equation for T (P ) -Equation B5/B6 -if gravity g and heat flux F are assumed to be vertically constant (see Appendix B1.)This assumption is reasonable if the P < P0 region is thin.Thus one only needs F (θ, T0) to specify the temperature-pressure structure at any angle θ.
The magnitude of any assumed inward heat transport by conduction would depend on the specific choice of P0.We believe that any inward heat flux will be small due to the inefficiency of conduction.So, to avoid the dependence on the arbitrary choice of P0, we set the heat flux at any angle where the surface temperature is higher than T0 to 0. Doing so should not affect the overall evolution (see §4.2.) The desired boundary condition for the Henyey scheme is the temperature at P0 as a function of the luminosity from the interior, i.e., T0(L).One must find the function F (θ, T0) that produces a given spherically symmetric T0.F (θ, T0) is linked to the luminosity, L, by integrating over the surface: where R is the total radius of the planet.
We now have all the information required to find F (θ, T0), under the previous assumptions.If one guesses the interior heat flux at a given angle, F (θ), then Equation 26gives the surface temperature.One can then use this as a boundary condition for integrating dT /dP (Equation B5/B6) to find T0, the temperature at P0.This gives a mapping between T0 and the flux, and so for a specified T0 the corresponding angular function F (θ, T0) can be found.Equation 28 allows us to find the function L(T0), which is the inverse of the relation we seek.
The function F (θ, T0) has to be computed numerically, and we do so by solving Equation B5/B6 to find T0 for a grid of F and θ.We use 5th-order Runge-Kutta integration, with an adaptive step size and a relative tolerance of 10 −6 .This grid is then fit to a spline and F (θ, T0) is found using Brent's method (Brent 1974).
We compute the boundary conditions for various gravitational accelerations, g, and substellar temperatures, Tss.The substellar temperature is defined by: and so is a measure of the stellar flux a planet receives.It is also the temperature of the point on the planet closest to the star given no redistribution of energy or heat flow from the interior.As the assumption of no redistribution is likely well justified for the short-period planets we consider, it is approximately the maximum temperature of the planet at late times when the internal heat flux is small.The resulting relation between luminosity and temperature at the edge of the 1D domain is shown in Figure 2 for a particular Tss. 2 , as detailed in Appendix B2.The basic features of the plot are that at high and low fluxes, the flux is a strong function of temperature, essentially due to black body cooling.Meanwhile, for a large range of intermediate fluxes, the temperature changes little.This is because, when T0 is close to the critical melt fraction, the viscosity at that point, and thus the amount of energy that can be transported, changes by orders of magnitude for very small temperature changes ( §2.4.3).

Grid and timestepping
We solve the equations in §2.1 on a mass grid where the amount of mass enclosed below cell j is mj ∝ j α .α is chosen to give a compromise between cells getting smaller in radius towards the edge, where high resolution is required to deal with the final crystallisation of the mantle, and maintaining a reasonable resolution near the centre of the planet.As density is close to constant, the first condition requires α ≲ 3 (α = 3 would correspond to cells of constant radius if the density is constant).For the second condition, it helps to have α > 1.We find α = 1.5 works well.
The timescales over which processes occur in planetary evolution differ vastly.We must therefore adapt the timestep in order to maintain accuracy and also allow the models to run in a reasonable amount of time.We employ a simple algorithm to do this.When undergoing the convergence steps of the Henyey scheme the timestep is fixed and is chosen at the start.As a first suggestion for the n th timestep, we use the formula: where L is again the luminosity and subscripts denote the timestep number, so ∆tn−1 is the previous timestep.This extends or shrinks the timestep such that the fractional change in the luminosity gets closer to fL.We use fL = 30%.
We then also predict the luminosity for this step using linear interpolation and use this to predict T0 for the current step, which is a function of L ( §2.7).We then impose that T0 does not change by more than a fraction fT , and the new timestep is given by We use fT = 2%.The timestep may be shrunk further to limit the amount of mass lost in this step, as described in the next section.Finally, for the first timestep, since there is no previous timestep to compare to, we use 0.1% of the Kelvin-Helmholz timescale.

Mass Loss
Since our aim is to model the evolution of rocky planets undergoing evaporative mass loss, we incorporate mass loss as follows.
At each timestep, we ascribe a mass loss rate determined by the current mass, radius and surface temperature.We use mass loss rates computed using the method described in Booth et al. (2023), which assumes that the dusty outflows form 1 µm grains with properties similar to forsterite.Example mass loss rates are shown in Figure 3.We pre-tabulate mass loss models for a given substellar temperature (Equation 29) for a grid of planet masses and densities.
When the mass loss rate has been determined, we then calculate the mass change to the next timestep according to our suggested timestep (see §2.8).In order to make the numerical problem easier to solve, rather than removing the mass from only the outermost cell, as would be the closest to the physical reality, we instead shrink a finite number of the outermost cells and keep the mass contained in cells interior to this the same.This prevents a drastic change in the size of the outermost cell, which might be numerically unstable, allowing more mass to be removed per timestep.Before doing so, we check that the mass loss results in the mass contained in these cells being reduced by no more than a certain fraction, and if the mass step is too great, we reduce the timestep so that this condition is satisfied.Our default setup is that the outermost 5% of total cells are shrunk and by no more than 1%.This procedure requires us to merge cells so that we can continue to remove mass after the cells shrink and to split cells, to ensure that individual cells do not become too large a fraction of the total planet's mass.We make these grid changes at the beginning of each timestep, when appropriate thresholds are reached and use piecewise cubic Hermite interpolation (PCHIP, Fritsch & Carlson 1980), where interpolation of values is required.
During mass loss, the gravitational potential of the planet is altered.Since mass loss occurs at the beginning of the timestep, this energy change is not included in Equation 3. We, therefore, introduce an extra source term to that equation, using the prescription used in MESA (see Paxton et al. 2019, §3.3).This approach considers the mass moving through and out of cells and the energy deposited by it, and we also assume that the mass loss timescale is longer than the thermal timescale.We generally find this energy contribution is small.

No mass loss
To demonstrate our model, we first show a fiducial case without mass loss, as the features of the thermal evolution are clearer in this case.We choose a substellar temperature of 2320 K as it is the fiducial value used in Booth et al. (2023) 3 , which is the source of our mass loss model ( §2.9).
Snapshots of the internal structure are shown in Figure 4 and the evolutions of some overall properties are shown in Figure 5.As the planet's temperature decreases (Figure 4, Panel a), it crystallises from the inside out (Panel c), which is essentially due to the shape of the solidus/liquidus relative to the temperature structure, which is close to adiabatic. 4As the melt fraction passes through the critical point of ϕ = ϕc = 0.4, a many order of magnitude change in viscosity is observed (Panel d and see §2.4.3).The increase in melt fraction close to the core-mantle boundary and decrease towards the surface seen at later times is because the mixing length is equal to the distance to the nearest boundary, gen- of a 0.15 M ⊕ planet with a core mass fraction of 0.3 and a substellar temperature of 2320K, but no mass loss.In the bottom panel, regions to the left of the liquidus are entirely liquid, to the left of the ϕ = ϕc line are partially molten, to the left of the solidus behave like a solid but have some partial melting and to the right of the solidus are fully solid.The region above the 1 GPa line is unresolved by the 1D model and will have melting on the dayside surface but be cold and solid on the night side.We also do not plot any partial melting at the core-mantle boundary, which does occur in our models, at certain times, due to the thermal boundary layer there (see Figure 4) because we are more interested in the state towards the surface.erating thermal boundary layers, which can also be seen in Panel a.In addition, the solid viscosity is highly temperature dependent (Equation 22).This results in a positive feedback: the viscosity increases towards the surface, due to lower temperatures there.This makes convection less efficient and the thermal boundary layer stronger, so the temperature at the surface falls, causing the viscosity to further increase.The effect is most noticeable at late times in Panels a and d.
The first two panels of Figure 5 show the evolution of flux and temperature, which are best understood with reference to Figure 2.There is an extended period of time where the temperature at P0 does not change, but the luminosity continues to drop due to the large changes in viscosity around ϕc.Note also that the flux levels off slightly at ∼ 10 8 yrs before decaying over the following 10 10 yrs, which is due to the radioactive decay of long-lived radioisotopes becoming the dominant heat source over primordial energy.
The final panel shows the planet shrinking, which is due to the density increase over time (see Panel b of Figure 4).Both the core and mantle shrink over the whole evolution, but the largest effect on the radius is the change from liquid to solid in the mantle, causing a marked decrease in total radius in the first few thousand years.As can be seen from the horizontal extent of the lines in Panel b of Figure 4, the density increase also results in a pressure increase at the coremantle boundary and centre of the planet.
As well as the planet shrinking, the final panel of Figure 5 also shows the molten state of the planet.One sees that, within a few hundred years, the planet is no longer fully molten (dashed red line) and, within a few thousand, the whole mantle is below the critical melt fraction (dot-dashed line) meaning that its behaviour is close to that of a solid.Partial melting does persist for Gyrs, but, beyond a few Myrs, it is very low (see blue dashed line that marks ϕ = 5%).
In short, despite the high substellar temperature, the planet is almost entirely solid for the majority of its life, because it can cool from the night side.What is not included in this one dimensional plot is that a magma pool will persist on the dayside, since the temperature there, which is approximately the substellar temperature (see §2.7), is sufficient to fully melt the rock (i.e., it is above the liquidus.)It is this magma pool that evaporation occurs from, leading to the observed mass loss.

Mass loss
Figure 6 shows the evolution of the total mass, radius and mass loss rate for a number of models with different initial planet masses and substellar temperatures.Mass loss rates increase with increasing substellar temperature and decreasing initial planet mass (as seen in Figure 3), which is seen in the faster mass and radius decrease for such planets.We mark on the mass loss rate panel an 'observable' cut-off mass loss rate of 10 −1 M⊕ Gyr −1 , motivated by calculations by Perez-Becker & Chiang (2013).We return to considering when the catastrophically evaporating planets are observable in more detail in §5.1.2.However, even with this simple condition, some interesting features can be noted about observed planets.Some planets are observable throughout their lifetimes, but their lifetimes are short (≲ 10 8 yrs).These are planets with sufficiently high substellar temperatures (dotted lines) or low initial masses (blue lines).At the other extreme, much cooler planets (solid lines) or massive planets (red curves) are never observable within the age of the universe.There is also an intermediate regime of planets that become observable late in their lifetimes once they have lost sufficient mass.This was also noted by Perez-Becker & Chiang (2013) from their earlier mass loss models.We discuss the consequence of this on the number of systems we might observe in §5.
In Figure 7, we show the evolution of the molten state of the mantle.In all the models, the mantle still entirely crystallises to a melt fraction below ϕc within ∼ 10 4 yrs, meaning the planet's interiors should be essentially solid by the time they are observed around Gyrs old main sequence stars.Furthermore, in many models, crystallisation occurs fully by around 1 Gyr.The exceptions are the highest temperature and lowest .Evolution of planets' mass and mass loss rates with different initial masses and substellar temperatures.Mass loss rates are calculated using Booth et al. (2023).All planets start with a core mass fraction of 0.3.mass cases, where the planet is evaporated before the melt fraction is zero.
However, even if the melt fraction is not zero, it is very small in all cases, and the planet still essentially behaves like a solid, as shown in the lower panel of Figure 8.At the same time, note that the melt fraction increases slightly at late times in some cases.The reason for this is subtle.The amount the planet's centre can cool is limited by the amount of energy that convection can carry away.The decrease in mass means there is less material above the lower mantle and core, allowing it to cool more easily.Thus the temperature of the core and lower mantle can decrease, causing an increase in luminosity and thus mean flux at the surface, as shown in the top panel of Figure 8.This increased flux induces more melting at 1GPa.This increased partial melting may have some effect on any volcanism on the planet and may alter the detailed composition through solid-melt partitioning, but since the amount of melting is small, it does not affect the overall conclusion that the catastrophically evaporating planets are essentially solid.
To summarise, these highly irradiated planets solidify significantly within a few thousand years and are completely solid within a few Gyrs.This is because they cool easily from their nightsides.Consequently, major melting occurs solely on the dayside, as proposed by e.g., Kite et al. (2016).Mass loss does not have a large effect on this result.

Depth of the magma pool and implications for chemistry
We found in §3 that catastrophically evaporating planets are likely to have almost entirely crystallised by the time they are observed.However, the peak surface temperature is sufficient to melt rock (higher than the liquidus), so there must be a molten region on the dayside.
The most naive way to calculate the depth of this pool is to assume a constant temperature gradient between the surface and a point in the interior where the temperature is known, i.e., conduction into the planet with constant conductivity.This neglects any convective heat transport in the magma pool (which we address below.)The depth of the pool in this case is the point along this temperature gradient where the material becomes solid, best defined as when the melt fraction reaches the critical value (ϕ = ϕc).We use T0, the temperature at a pressure P0, the edge of our 1D model (see §2.7) as the interior temperature, and thus the depth is dependent on what we choose for P0, as shown in Figure 9 (top panel).Even this naive calculation shows that the pool is shallow relative to the planet's size (≳ 1000 km), but, as we shall describe below, including circulation makes the pool even shallower.Kite et al. (2016) argue that lateral circulation in the lava pool greatly reduces the depth of the pool, compared to the above estimate.The reason for this is the thermal structure created by horizontally driven convection.Convection in the pool is driven by the day-to-nightside temperature gradient at the surface.5This configuration is known through experiments and theory to create a thin thermal boundary layer at the surface, as shown in Figure 10 (e.g., Vallis 2006;Hughes & Griffiths 2008).Using scaling laws from Vallis (2006), the depth of this boundary layer is (adapted from Kite et al. 2016, Eq. 8).Here R is the planetary radius, Ω its orbital frequency, θp the angular size of the pool centred on the substellar point, and ∆ρ the density change across the boundary layer.Since there is a large temperature drop across the thermal boundary layer at the surface, the point where ϕ = ϕc will be shallower than in the naive case.
We adopt the assumption in Kite et al. (2016) that the total depth of the pool should not be more than ∼ 10 times this boundary layer.θp may be found, for a given substellar temperature, Tss, by using the fact that the surface temperature varies approximately as cos 1/4 θ (see §2.7), and finding the point that the critical melt fraction is reached on the surface.Tss ∼ 2100 K, representative of the known systems KOI 2700b and K2-22b (Table 3), gives, using our melting curves ( §2.2), a very wide pool with θp ≈ 0.4π.If one takes surface values of the physical quantities, assuming ∆ρ = 10%, a typical value for melt-tosolid density change, at a one-day orbit, then one finds the scaling Since for an individual planet, the other terms in Equation 33do not change, Equation 34 may be used to track the evolution of the pool's depth for that planet.
We plot the evolution of the pool depth, assuming it is equal to 10 δT , for one case in Figure 9 (bottom panel).This demonstrates the scale difference between this case, including circulation, and the simple argument without, as noted by Kite et al. (2016).
It is interesting to note the different evolutions for the two cases.At early times the depth of the pool with both estimates decreases due to cooling and contraction of the planet.The contraction causes the gravity to increase; thus, the pressure at which the surface lava becomes solid is at a shallower depth.At later stages, the naive estimate increases, whereas the estimate with circulation decreases.This is because, at this point, mass loss is the most important factor.For the naive estimate, the corresponding decrease in gravity means that the depth to a certain pressure increases, meaning the pool depth increases.On the other hand, for the estimate with circulation, the pool's depth decreases as the planet gets smaller and denser.This is because the estimate is coupled to the size of the pool's thermal boundary layer which scales with mass and radius in a different way (Equation 34).However, it should be noted that the fractional changes for both cases are relatively small, so the true time evolution may well be altered by effects we have not captured, for instance, compositional changes in the pool.
The important point, however, is that, in either case, it is likely that the planets have only a shallow molten region on their dayside.This has consequences for the composition of the dusty tails, which are the observational signatures of catastrophically evaporating planets.The evaporation that generates the tails must come from the molten region of the planet.Our results show that this is only a shallow region at the surface at any time, which gradually moves into the planet.Were the planet fully molten, relatively volatile el- Simulations that end with a dot had full crystallisation up to 1GPa.
Those that end with a cross reached the maximum bulk planet density that mass loss rates were computed to.
ements such as Na would be lost.However, with a shallow pool, it is possible that they are locked in the solid portion of the planet, which is gradually melted away, and so are still present in the winds.A more detailed consideration of the chemistry is required to determine the compositional evolution, which will come in a future study.

Assumptions in the boundary conditions
As discussed in §2.7, due to the complexity of the radial and angular heat distribution of the star's energy, we made some assumptions about the outer layers of the planets.We argued that the inward heat fluxes should be small compared to the outward ones, which justified having the inner regions of the planet evolve as a spherically symmetric structure, unaffected by the star.At very late times, when all internal energy has dissipated, this assumption is clearly untrue.Therefore, in the following, we shall check whether including any inward or angular fluxes would affect our conclusions.At very late times, for tidally locked planets, the temperature should decrease from the dayside to the nightside, and energy should flow through the planet.Heat flux, in this case, should be of the order where Tn is the nightside temperature, and k is again the conductivity.For a 1000 km planet with Tss = 2320 K, and a negligible nightside temperature, this gives a value of ∼ 5 × 10 −3 Wm −2 .The blackbody temperature required to emit this on the nightside is 17 K, justifying the assumption of the nightside temperature being negligible.One can see in Figure 5 that even after 10 Gyrs, the outward fluxes do not get this low, meaning the assumption of no inward flux is justifiable up to this point.Furthermore, the planet crystallised significantly before this point when the outward flux was at least 20 times higher.Therefore, the additional flux contribution would not be enough to melt the planet.We also assumed that angular fluxes were small, which allowed us to estimate the flux from the nightside simply as that coming from the interior, not the dayside.Angular fluxes in the solid planet should also be of a similar order to the latetime day-to-nightside flux, estimated above; thus those fluxes become important at a similar time.Consequently, any adjustment to the nightside temperature would only become significant once crystallisation has occurred, meaning angular fluxes also cannot affect our conclusion about the planet having solidified.
We have so far only shown examples where the pressure at the edge of the 1D interior model is fixed at P0 = 1 GPa.We chose this pressure as regions deeper than it are unlikely to be affected by stellar heating.A demonstration of this is that assuming conduction inwards from the surface, the flux into P0 will be approximately: with quantities as defined in §2, and evaluated at the surface.Taking a typical value for our planets of g = 4.5 ms −2 and Tss = 2300 K and assuming T0 = 1500 K, which is below the solidus, gives ∼ 0.05 Wm −2 .This is less than the amount of outward flux at the point when the planet has fully crystallised, meaning at least up to that point P0 = 1 GPa is deep enough to assume inward heat flux to deeper in the planet is small.Furthermore, there is heat redistribution in the pool, meaning this is likely an overestimate of any inward flux, as addressed in the previous section.
We have justified that this value of P0 is valid in terms of having a low enough inward heat flux through it.However, the exact choice is still arbitrary, so it would be useful to test the effect that choosing a different valid value would have.If the value was much deeper, the assumptions of constant gravity and flux in the outer region become less accurate since more than 10% of the total mass of smaller planets would be included in the boundary region.6Therefore, to investigate the effect of choosing this arbitrary pressure, we instead compared to a lower pressure case of P0 = 0.5 GPa.We found that, in this case, deep mantle temperatures (close to the core-mantle boundary) change by less than 0.5%, and upper mantle temperatures (around 1GPa) change by less than 5%.However, because the temperature of the upper mantle does not change much between 10,000 years and several Gyrs (Figure 5, middle panel), and the viscosity is strongly dependent on temperature, the time at which crystallisation occurs is nevertheless sensitive to the small change.Crystallisation of the mantle below ϕc takes a few hundred years longer, and full crystallisation (the whole mantle having ϕ = 0) takes a few hundred Myrs longer (differences of ∼ 10%.)This effect is a consequence of some simplifications made in the outer layer with pressure below P0.The outer layer does not include any energy production due to cooling or radioactive decay.Therefore, in the 0.5 GPa case, the planet has more energy available due to the outer layer being a smaller proportion of the planet and so cools slightly slower.
It should be noted, however, that beyond 10 Myrs, while the melt fraction in the mantle is non-zero, it is very low (Figure 5, bottom panel and Figure 8), so the planets are essentially solid.In Figures 5 and 7, we plotted ϕ = 5%, to represent a low melt fraction that could probably be considered solid.This would have to be delayed to 1 Gyr to alter the conclusion that planets are solid when they are observed evaporating around main sequence stars.This point is delayed by a few 10s of Myrs by halving P0, so we deem it unlikely that our model differs sufficiently from the physical reality to change this overall conclusion.

Further caveats
We have used a relatively simple model to investigate the evolution of the catastrophically evaporating planets.Thus there is necessarily some physics missing.For instance, the composition is fixed and uniform in the mantle.In reality, there may be planet-to-planet variation and compositional evolution within the planet.However, these are unlikely to change the density and heat capacity properties enough to alter our conclusions significantly.The composition would also affect the solidus and liquidus, which may make a more significant difference due to the difference between melt and solid having such a large effect on viscosity (see §2.4.3).A full investigation of this is beyond the scope of this work, but given the whole mantle passes below the critical melt fraction fairly early, such differences are likely unimportant over the planets' full lifetimes.We also assumed that the melting is in equilibrium, meaning temperature and pressure directly determine the melt fraction.Non-equilibrium effects are potentially important for compositional evolution during the early magma ocean phase (e.g., Solomatov 2007), but once again, since it passes through the critical melt fraction point early on, these effects are unlikely to be important for long term evolution.
Another simplification we made is in the treatment of the iron core.Since we are most interested in the mantle's evolution, the core's precise state is not important, only the energy it provides to the mantle.In reality, the core will start molten and solidify over time, which generates extra energy through the latent heat of fusion and through gravitational potential energy as the iron density increases through the phase transition, resulting in the contraction of the core.However, we do not include such a phase transition.
For our example in Figure 4, the melting point of iron at the core-mantle boundary is around 2450K (Sluiter 2012).This means that the planet must be between 10 8 − 10 9 years old, and the mantle is already essentially crystallised by the time any of these extra heat sources start to come into effect.
The total energy released as latent heat may be estimated as EL ∼ Mc∆LF e, where ∆LF e is the latent heat of fusion of iron, and Mc the mass of the core.The energy from gravitational contraction can be estimated with Eg ∼ GM 2 c /Rc( 1 3 ∆ρc/ρc), with Rc the core radius, ρc the density of iron and ∆ρc the solid melt density contrast.Typical values for, for instance, a 0.15M⊕ planet with a core mass fraction of 0.3 give total values of several 10 28 J for these two energy sources.This is of order the amount of energy lost by the core by cooling, as calculated by our models.Thus these additional energy sources may have some effect on the mantle at late times.However, it seems unlikely that the luminosity of the core should increase at this point, as the energy would instead re-melt the core.Furthermore, the mantle itself typically has energies an order of magnitude above that of the core, both in thermal energy and radioactive elements.Therefore, it seems unlikely that including these core freezing effects would affect the fact that the mantle has crystallised first.
Another energy source we have neglected is tides.The short orbits of these planets mean that tidal forces are potentially strong, although heating may not be long-lived because circularisation timescales are consequently short.Jackson et al. (2008), suggest that heating rates can be high even for small but non-zero eccentricities, perhaps even high enough to keep the mantle molten.Considering these effects would require consideration of the orbital evolution, which is beyond the scope of this work.
Orbital evolution could also occur due to angular momentum exchange between the outflow and the planet, thus changing the planet's irradiation level.This might result in the planet spiralling outwards or inwards, cutting off mass loss or causing a runaway.Neglecting this effect is justified if the material from the planet is expelled from the system by the star, retaining its own angular momentum, meaning no angular momentum is exchanged with the planet.This fits with the current understanding of the dusty tails (e.g., van Lieshout & Rappaport 2018).Furthermore, Perez-Becker & Chiang (2013) showed that there is little orbital evolution even when trying to maximise the amount of angular momentum exchanged.

OCCURRENCE RATE OF LOW MASS PLANETS
Our evolutionary models demonstrate that there is only a small period of time, relative to the main sequence lifetime of stars, that catastrophically evaporating planets have detectable mass loss, but have not yet completely evaporated, as also noted by Perez-Becker & Chiang (2013).The three systems that were found in the Kepler/K2 data must be in this phase.In §5.1 we estimate the occurrence rate of planets that evaporate at detectable rates, at some point in their lifetimes, needed to explain the number of detections.
The planets with detectable mass loss must be part of a wider distribution of low-mass planets, which is unlikely to be peaked at the temperature and mass that gives detectable mass loss.In §5.2 we assume this distribution is a power law in order to extend our occurrence rate to a wider range of lowmass planets.This allows comparison of the sub-Earth mass planet population, which we can probe but is generally undetectable by conventional exoplanet surveys (e.g.Christiansen et al. 2014;Fulton et al. 2017), with observed occurrence rates for higher mass planets.
The inferred occurrence rates will depend on the mass loss model used since it changes the duration that planets are observed and the threshold of detectability.We present calculations with the same mass loss model we have used throughout (Booth et al. 2023).

Occurrence rate of planets with detectable mass loss
We can crudely estimate the number of planets that must be detectable at some point in their lifetimes, in order to produce the number of observed evaporating systems, using the equation Here N detect is the number of planets detected, np is the unknown number of planets per star, N * is the number of stars surveyed and fevap is the fraction of the host star's lifetime that the planet is observable.Ptrans is the probability that the system is transiting in the right orientation to be observed from the Earth and is given by (Borucki & Summers 1984) where R * is the stellar radius and a the orbital separation.
To find fevap one must identify the phase in a planet's lifetime for which it is detectable.We address this over the following two subsections and then use it to estimate the occurrence rate from Equation 37.

Fitted function to the mass loss rate
Whether a planet is observable is dependent on its mass loss rate, and so to find fevap we must determine the mass loss rates of planets uniformly over temperature and mass.Our mass loss models were only run for a finite number of temperatures and initial masses, so we require a method of interpolating between them.We choose to fit the evolution of our models in mass-mass loss rate space (M − Ṁ ) using the function where Minit is the initial mass of the planet and all other parameters are fit for and depend on substellar temperature.The physical reasoning behind this form is as follows.The mass loss rate at any time, for a given substellar temperature, depends on the instantaneous mass and core mass fraction.In turn, the core mass fraction depends only on the current mass and initial mass, if the initial core mass fraction is assumed to be the same for all planets.Thus the mass loss rate is just a function of M and Minit.At low masses, the mass loss rate increases with mass, which we assume to act as a power law (first term in curly brackets, with a > 0.) At high masses, the mass loss decreases with increasing mass, approximately exponentially, hence the second term in curly brackets.These low and high mass behaviours are then combined using the smoothing parameter γ.In principle, more dependence on M or Minit could be introduced, for instance, the first term could include an Minit dependence, but we found that the above form could capture the behaviour without introducing extra degeneracy.The advantage of using this physically motivated form is that it not only allows interpolation but also extrapolation to higher core mass fractions with more confidence than an interpolation method such as a spline.The mass loss models become more numerically difficult to run at higher core mass fraction, and therefore we did not run them to a core mass fraction of 1 when the entire mantle has evaporated.
For each substellar temperature, we fit Equation 39 to a selection of our numerical outputs with different initial masses and at different evolutionary stages using scipy curve_fit.These are shown in Table 2. To find the mass and mass loss rate evolution for a given initial mass and substellar temperature, we integrate this function numerically.We find that this can reproduce our full calculations to better than 8%, even for initial masses that the function was not fitted with.For the purposes of our occurrence rate calculation this is sufficient given the other uncertainties and assumptions.To interpolate between substellar temperatures we use spline interpolation of the parameters in Equation 39/Table 2.

Length of time that systems are observable
To work out the times when the planets are observable, we solve the ODE in Equation 39 over a range of substellar temperatures and initial masses to find the mass loss rate evolution using two definitions of the detectable mass loss rate.Firstly, an optimistic case (henceforth Case 1) where we consider planets detectable if their mass loss rate exceeds 0.1M⊕Gyr −1 , which is considered a lower limit on detectability by Perez-Becker & Chiang (2013).In contrast, we also consider the threshold to be the region of temperature and current mass where dust is produced in the models of Booth et al. (2023) (henceforth Case 2), shown in Figure 3.This should, in principle, be the more physical case.We do not consider the evaporation of the core after the mantle has evaporated because the composition of the dust tails is thought to be inconsistent with iron (van Lieshout et al. 2014(van Lieshout et al. , 2016)).Therefore, we are technically only considering the progenitors of planets with observed evaporating mantles.That being said, iron likely has a much higher mass loss rate (Perez-Becker & Chiang 2013) due to its higher vapour pressure, and so a planet's evaporation should accelerate once the iron core is exposed.Therefore the lifetimes we calculate are likely reflective of planets having both their mantles and cores evaporated.
The times at which systems start and cease to be observable as well as the length of time that they are observable are shown in Figures 11 and 12 for Case 1 and 2 respectively.In both cases for lower-mass planets, the mass loss is observable from birth (see top panels), and so the length of time a planet is observable is limited by the time it takes either to evaporate or for the mass loss rate to decrease below detectability.As a result, the length of time that the planet is detectable increases with initial mass (lower panels).
Once planets are massive enough, they take time to become observable, and the detectability is now limited by the length of time it takes to become observable, which increases with initial mass.In both Case 1 and 2 the time it takes to become observable and the time it takes to lose the entire  11.The time period of observable mass loss for planets with core mass fractions of 0.3, mass loss rates calculated using the method of Booth et al. (2023) and assuming that a planet becomes observable if its mass loss rate exceeds 0.1 M ⊕ Gyr −1 (Perez-Becker & Chiang 2013).Top panel: The dotted line shows the age at which a planet of a given substellar temperature gains a mass loss rate high enough to be observable.The dashed line shows the time when the whole of the planet's mantle evaporates, while the dashdot shows when the mass loss rate becomes unobservable, which mostly coincides with the mantle evaporation line.Bottom panel: The length of time that a planet has an observable mass loss rate as a function of initial mass and substellar temperature.mantle seem to converge (top panels), which is conceptually because at high masses it is possible to lose the entire mantle without ever reaching a high enough mass loss rate to be observed.The two lengths of time are only actually equal before the age of the universe for the lowest temperatures.The peak length of time occurs when a planet starts small enough that it has a detectable mass loss from birth but it is large enough that it evaporates for a long time.Planets that start with higher masses spend less time observable because at the point when they reach this mass, the mass of the mantle is lower since all planets start with the same core mass fraction.For Case 1 the peak length of time moves to higher masses as the temperature increases.This is because the higher temperatures allow larger-mass planets to have observable mass loss rates.This effect is drowned out by other temperaturedependent effects for Case 2, which we will discuss below.
A major difference between Case 1 and 2 is that for Case 1 the timescales that systems are observable for increase with increasing temperature (Figure 11), whereas for Case 2 the timescales decrease for high temperatures (Figure 12).This is  due to the different assumptions behind the two cases.Firstly, in Case 2 dust is no longer produced at higher temperatures.
The physical reason for this is that the dust condensation rate is too low compared to the flow rate.Hence dust simply doesn't have time to form before the gas leaves the planet, as discussed further in Booth et al. (2023).Furthermore, for higher temperatures, mass loss rates increase, meaning planets spend less time in the observable region (see Figure 3).These two effects result in the timescale decrease with temperature for Case 2, seen in Figure 12.In contrast in Case 1 the observable timescale increases with temperature.This is because higher temperatures allow more massive planets to have observable mass loss.Larger mass planets take longer to fully evaporate, so the overall effect is that the timescales increase.This also means that the range of initial masses that can be detectable is much higher for Case 1, where there is no high-temperature cut-off.
A final difference for Case 2 is that the observable timescales are overall shorter since the mass loss rate required to be detectable is, at its lowest, similar to but generally higher than the fixed value for Case 1.
For the calculations we show, we assume that the progenitor planets all have core mass fractions of 0.3.If there were a range of core masses the amount of time the mantle is evaporating would change.It is not immediately clear whether it would increase or decrease because a larger mantle has more material to evaporate, but evaporates quickly as it has a lower gravity for the same mass.We briefly investigated this by using models with a very low core mass fraction of 0.1.We found that these two effects almost cancel out and that the timescales of observable evaporation are very similar, although the lower core mass fraction examples become observable earlier, which would slightly decrease their observability around old stars.Consequently, we do not believe this will greatly affect our analysis in the following sections.
5.1.3Basic estimate of the number of planets that will have detectable mass loss fevap in Equation 37 may be estimated as the ratio of the observable period to a typical stellar age.Taking a typical age to be ∼ 5 Gyrs, and a typical observable time (Figures 11 and 12) to be around 200 (80) Mrys for Case 1 (2) gives fevap ∼ 0.04 (0.016).The primary Kepler survey observed ∼ 180, 000 stars.Of these ∼ 145, 000 are main sequence stars, among which two catastrophically evaporating systems have been identified (Kepler 1520b and KOI-2700b -we disregard K2-22b as it was part of the separate K2 survey).Therefore, we take N detect = 2 and N * = 145, 000.We further use a value of Ptrans = 0.25, for the archetypal system Kepler 1520b (properties can be found in Table 3).Equation 37 then implies that the required number of planets per star that undergo observable evaporation at some point in their lifetimes, in order to reproduce the detected number of these planets, is np = 0.14% (0.3%) for Case 1 (2).This is essentially the same as the estimate of 0.1-1% in Perez-Becker & Chiang (2013) since the detectable timescales are very similar.

Probability of detection across the stellar population
The calculation above uses just a single value for Ptrans, which is a function of stellar properties.Here we improve our estimate by taking into account the full variation of stars that were observed.Kepler's primary mission observed stars with a distribution of masses (M * ) and ages (t) which we will write as n * (M * , t), defined as the number of stars in the range [M * , M * + dM * ] and [t, t + dt].
One can work out the probability of detecting a planet catastrophically evaporating around the average star, as a function of the planet's initial mass and substellar temperature (Equation 29), by integrating over the stellar population: Here n * (M * , t) is the normalised distribution of stellar ages and masses in the sample, i.e., n * (M * , t)/N * , where N * is again the total number of stars.Ptrans(a, R * ) is as defined in Equation 38.t obs,min/max (Tss, Mp) denote the earliest and latest times that the evaporation of a planet is observable, which are naturally functions of the planet's properties and were calculated in §5.1.2(Figures 11 and 12).This integral is a generalisation of the product fevapPtrans in Equation 37. To find the distribution of stellar masses and ages in the Kepler sample (N * (M * , t)) we use the stellar properties derived from GAIA in the GKS sample (Berger et al. 2020) and apply the colour cut in Fulton et al. (2017) to remove giant stars.To calculate R * and L * we use the MIST stellar isochrones (Dotter 2016;Choi et al. 2016;Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015) ) with solar metallicity.While this neglects metallicity variation, our simple model does not warrant this complication.In order to evaluate Equation 40 we binned the stars in the GKS sample according to stellar age and first integrated over the distribution of stellar masses for each age bin before integrating over the age distribution.Each distribution was normalised, thus giving the probability per star.
P detect (Mp, Tss) is plotted in Figure 13, for both the cases we consider.In the top panel, we see the result of the trends identified for Case 1 in §5.1.2.Higher-mass planets are observable at higher temperatures and it takes longer for the higher-mass planets to evaporate, thus they are more likely to be observed.However, the contours show two ridges for a given temperature.The ridge at lower initial mass is due to the peaks in the length of time observed, as seen in Figure 11.
However, stars that are several Gyrs old are more common, which means at a population level it is more likely to observe a larger initial mass planet because they are more long-lived, even if the time they are observable might be shorter.This creates a second ridge at higher initial mass.The upper value of P detect ∼ 0.01 also corresponds to the rough estimate in §5.1.3since P detect ∼ fevapPtrans ∼ 0.04 × 0.25.The bottom panel of Figure 13, meanwhile, reflects the trends for Case 2 ( §5.1.2,Figure 12), with detectability decreasing for higher temperatures.The consequence of this is that there are cuts of observability at both low and high temperatures.In this case, even the peak probability of detection, ∼ 3 × 10 −4 , is lower than that roughly estimated in §5.1.3 of P detect ∼ fevapPtrans ∼ 4 × 10 −3 .This is because the observable planets are limited to cooler temperatures, and are thus further away and so less likely to be transiting (Equation 38).
One can use these refined probabilities to redo the calculation in §5.1.3.For Case 1 the value of np is unchanged at 0.14%; for Case 2, taking a typical value of P detect ∼ 2×10 −4 , gives a requirement of ∼ 6.8% of stars having a planet which can have observable catastrophic evaporation in its lifetime.These values of np, by definition, refer just to planets that have observable mass loss in their lifetimes, i.e., just the nonzero probability regions in Figure 13.Thus for Case 1 the value refers approximately to planets of 0.1 < Mp/M⊕ < 0.3 and 2200 K < Tss < 2600 K and for Case 2 it refers to 0.04 < Mp/M⊕ < 0.08 and 1950 K < Tss < 2200 K.This will be an important consideration in the next section.

General occurrence rate of low mass, short-period planets
Currently, our calculated occurrence rates only encompass planets that will drive observable mass loss.However, it is unlikely that low-mass planets only exist at the temperature and sizes that produce observable mass loss.Thus, the observed catastrophically evaporating planets are likely to be part of a wider population of low-mass planets.Therefore, we can use our models of when planets do and do not produce observable mass loss to estimate the general occurrence rate of low mass, short-period planets7 of low-mass, short-period planets, assuming their orbits remain unchanged.. Since not all low-mass, short-period planets will give rise to detectable mass loss, this general occurrence rate is necessarily larger than our previously calculated values.At the most basic level, this increase in planetary occurrence rate from just the observable population to the general population will scale with the ratio of parameter space volumes occupied by the general population compared to the observable population.Since the parameter space volume occupied by the observable planets is small (i.e.we require a narrow range of planet masses and substellar temperatures), the occurrence of the general population will be significantly larger than just the observable catastrophically evaporating population.

Calculation
In order to deduce the general occurrence rate of low-mass planets from the catastrophically evaporating planets, one needs to know the region of parameter space that gives rise to detectable mass loss and how this "observable" region compares to the underlying population.We have already worked out the shape of the observable region for two different cases in §5.1.4(Figure 13).However, the distribution of the underlying population is unknown, but it is probably fair to assume it is reasonably smooth and not peaked exactly at the region where mass loss is observable.
We define the number distribution of the general low mass, short-period planet population, np(Mp, Tss), as the number of planets per star with masses and substellar temperatures in the ranges [Mp, Mp + dMp] and [Tss, Tss + dTss].Following both observed populations (e.g., Petigura et al. 2022) and population models (e.g., Makino et al. 1998), we choose to model it as a power law: where C is a normalisation constant, and α and β are indices that we will investigate different values of.This is implicitly assumed to be universal over stellar type and age.We choose this parameterisation because Mp and Tss are independent variables in our model.In a generalisation of Equation 37, the number of detections can be written as np(Mp, Tss)P detect (Mp, Tss) dMpdTss . ( where we highlight now that P detect can be zero when a planet with a given mass and surface temperature never produces observable mass loss.Thus, if we assume some α and β, we can evaluate this integral and infer the normalisation factor, C, that is needed to match the number of detections, N detect , allowing us to infer the number of underlying planets as a function of the assumed shape of the distribution. As explained at the start of this section, the idea is to extrapolate from planets in the region of parameter space with observable mass loss to the general population, so the detectable mass loss region must be well defined.This is not the case for Case 1, where the detectable region continues to extend to higher substellar temperatures (Figure 13, upper panel).However, if the stellar irradiation is too high, dust should never be able to condense.Therefore, to perform the calculation for Case 1, we impose a maximum sub-stellar temperature, above which it is unlikely for dust particles to exist.In reality, this cut-off should be compositionally dependent, but we simply set it to be 2500 K.It is important to note that our occurrence rate calculation for Case 1 explicitly depends on this choice.
It is instructive to convert our inferred number distribution into a measure that can be compared with other occurrence rate calculations.We choose to compare our results to calculations in Petigura et al. (2022), which is an up-to-date analysis of the Kepler data, using GAIA derived stellar properties.They present occurrence rates as a function of planetary irradiation, Firr, which is equivalent to substellar temperature through Equation 29 (Firr ≡ Firr(0)), and also give an analytic fit to their occurrence rate.Their occurrence rate is defined as 0.5 ∂Np/∂ log Firr for planets between 1-1.7 R⊕, which is roughly a 0.23 dex bin in radius.Thus, we choose to integrate over the same logarithmic range in radius, but centred on 0.2M⊕ (∼ 0.6R⊕) for Case 1 and 0.06M⊕ (∼ 0.4R⊕) for Case 2. These values are chosen because they are close to the peaks of P detect for each case (Figure 13).To convert from a radius range to a mass range, we assume a simple mass-radius relation of M ∝ R 4 (e.g., Valencia et al. 2007), which means that the Petigura et al. (2022) range of 0.23 dex in planet radius corresponds to 0.92 dex in mass8 .
The resultant values for the occurrence rate are shown as a function of the power law indices in Figure 14.The axes of Figure 14 attempt to cover a range that should encompass the physical reality.The x-axis, α, varies the dependence on initial planetary mass, and we ensure that we encompass values of ∼ −8/3 suggested for runaway growth through planetesimal accretion (e.g., Makino et al. 1998).The y-axis shows the dependence on substellar temperature, Tss.Petigura et al. ( 2022) measure the dependence on incident flux to be ∼ F −2.6 irr to F −3.1 irr corresponding to β ∼ −7.4 to −9, so we extend to values beyond this.
A very notable feature of Figure 14 is the strong variation in the occurrence rates with α and β for Case 1, but not for Case 2. The reason for this is the area of Mp-Tss parameter space that the two cases cover.The occurrence rate, 0.5 ∂Np/∂ log Firr, is essentially just a scaling of np evaluated at the peak of P detect .Since the non-zero region of P detect is so small for Case 2 (bottom panel of Figure 13), P detect effectively acts like a delta-function.This means that the integral Equation 42reduces to approximately N detect ∝ np, with np evaluated at the peak of P detect .Thus, the occurrence rate has a weak α and β dependence.In contrast, for Case 1, the non-zero area of P detect is much larger, meaning the variation of np over that area is much more important, so the occurrence rate has a much stronger dependence on α and β.
As expected, the general occurrence rates (Case 1: 0.5-20%, 2: 60-120%) are much higher than those calculated in §5.1.4(Case 1: 0.14%, 2: 6.8%), since we are now considering the number of planets per star in a larger region of parameter space.The values in §5.1.4are only true for the non-zero P detect region (Figure 13), whereas our more general measure, 0.23 dex in radius and 0.5 dex in Firr corresponds to 0.07-0.58M⊕and 2115-2693K for Case 1 and 0.02-0.17M⊕and 1775-2367K for Case 2, which is a region at least an order of magnitude larger than the region populated by observable catastrophically evaporating planets.

Comparison to other occurrence rate values
In Figure 15 we compare our occurrence rates with that found by Perez-Becker & Chiang (2013) for the catastrophically evaporating planets; with super-Earth (1-1.7 R⊕, ∼ 1-8M⊕) observations from Petigura et al. (2022); and with sub-Earth mass planets of radii 0.5-0.75R⊕ (0.6-0.3 M⊕) from Hsu et al. (2019).We convert the values from Hsu et al. (2019) (their Table 2) to the same measure as Petigura et al. (2022), assuming a 1M⊙ star (the median stellar mass for the GKS sample), with solar metallicity, as they only present results as a function of orbital period.The value of 0.1-1% derived by Perez-Becker & Chiang ( 2013) is only true for planets that will have detectable mass loss at some point.Therefore, for a valid comparison, we need to define what the detectable region of Tss-Mp space is, just as we have done for our own values above.We scale the Perez-Becker & Chiang (2013) values assuming the Case 1 observable region.Specifically, we scale them according to the ratio between the parameter space areas used by Petigura et al. (2022) (0.23 dex in radius, 0.5 dex in Firr) and that where P detect > 10 −4 (Figure 13, upper panel, limited again to Tss < 2500 K.) Comparing firstly our general occurrence rates to that found by Perez-Becker & Chiang (2013), we see in Figure 15 that for Case 1, the estimate is very similar.This is because we have simply scaled the values for planets that will have observed mass loss found in §5.1.4,which were very similar, to a larger area of parameter space.The only difference is the large variation caused by the uncertainty in the power law distribution, which we do not take into account for the Perez-Becker & Chiang (2013) values.Case 2, on the other hand, is around an order of magnitude higher.This was addressed in §5.1.4,where we also found the number of planets per star that will undergo catastrophic evaporation to be an order of magnitude higher (7% compared to 0.14%).This was mainly due to the fact that the dust model we use suggests detectability only at lower temperatures, which have lower transit probability, which means the occurrence rate must be higher to produce the same number of detections.
Compared to the observationally derived occurrence rates for both super-Earths and sub-Earths, the rates derived by both Perez-Becker & Chiang (2013) and our Case 1 are a little higher.However, the difference is insignificant and would suggest that the catastrophically evaporating planets are just the low-mass analogues of the same population.For Case 2, however, our inferred occurrence rate is significantly higher than that observed.
This very high occurrence rate prediction prompts the question of whether, if it is true, this population of low-mass planets should already have been detected by transits.It is difficult to detect sub-Earth-sized planets, but not impossible.Using the methodology of Fulton et al. (2017), which considers the signal to noise over the Kepler sample and includes the geometric probability of a transit (Equation 38), we estimate the chances of finding a planet of 0.4R⊕ to be ∼ 3×10 −3 at our typical substellar temperatures of ∼ 2000 K.According to our models, many of these planets should have been destroyed within 0.1 Gyrs, with the exact fraction depending on the shape of the number distribution.However, even taking the most extreme values of α and β, the number of planets that should be detectable in transit is only reduced by ∼70%.If our Case 2 occurrence rate is true, this should result in at least 120 of the ∼ 140, 000 main sequence stars in the Kepler sample having detectable Mercury-sized planets at around 1 day (where the catastrophically evaporating planets lie).However, there are no confirmed planets in the Exoplanet Archive (Akeson et al. 2013), excluding the catastrophically evaporating planets themselves, within the temperature and radius range.Since the Case 1 and Perez-Becker & Chiang (2013) estimates are both around 10-100 times lower, they do not imply such a strong inconsistency.That being said, the Case 1 model is clearly not perfect as it predicts many more planets at higher temperatures than lower, which is contrary to what is observed, which we return to in §5.3.
There are two main ways of reconciling the high value implied by Case 2. Firstly, it could be that the limited range of temperatures and masses that dust is produced at in the models is too small.While the models of Booth et al. (2023) are an advance on previous models of dust production in the systems' outflows, there remain various simplifications, such as a singular dust grain size.
On the other hand, if we take the models as true, the high occurrence rate can be explained by the observed evaporating planets forming further out, then being scattered into their observed orbits or reaching their current orbital period via tidal decay.Under our assumptions above, the observed catastrophically evaporating planets form at their current locations and at the same time as their host stars.Our occurrence rate refers to the progenitor population, and so if the orbits have changed, the planets can, firstly, have originated at a larger range of temperatures.Furthermore, they can also have started with lower initial masses since, under the current assumptions, low-mass planets are destroyed before their mass loss can be observed.This means that the overall parameter space that the progenitor population exists over would be increased, meaning the inferred occurrence rate would not be so high.We leave the investigation of scattering and tidal inspiral to future works.

The known systems in the context of our calculations
In Table 3, we show the properties of the three known catastrophically evaporating planets and their host stars.We also plot their substellar temperatures in Figure 13.As ought to be the case, our model suggests they all have detectable mass loss (P detect > 0).However, for Case 1, the probability of systems being detectable is strongly peaked towards high temperatures, whereas no such effect can be seen in the observed systems.This may, therefore, act as support towards the dust production models of Booth et al. (2023) being closer to the physical picture since such models predict a higher chance of detecting lower temperature systems, so they more closely resemble the data.The driving effect is the cut-off of dust production at high temperatures, which is probably independent of the details of the calculations.
This explanation for the temperatures of the systems is, of course, not unique.At least one alternative explanation could be that planets at very high substellar temperatures are intrinsically rare.This could be because hotter stars host fewer rocky planets, but the period distribution is the same.Alternatively, it could be because the occurrence rate decreases at lower orbital periods, thus higher substellar temperatures.This fits with the observed period distribution of Super-Earths (Sanchis-Ojeda et al. 2014;Petigura et al. 2022), which has a drop off at lower orbital periods, and makes sense theoretically as short period planets can inspiral onto their stars due to tidal interactions (e.g. Lee & Chiang 2017).
Another potentially interesting feature of the observed systems is that they are all K or M-dwarfs.The chances of pick-  41) for Case 1 and 2. We define our occurrence rate in the same way as Petigura et al. (2022) as 0.5 ∂Np/∂ log F irr for a 0.23 dex bin in radius (see text).Np is the number of planets per star and F irr the stellar irradiation at the planet's orbital radius.Note that for Case 1 the scale is logarithmic, while for Case 2 it is not.Our measure is evaluated at Tss = 2400 K, Mp = 0.2M ⊕ for Case 1 (∼ 0.7R ⊕ , F irr ∼ 1400F ⊕ ) and Tss = 2050 K, Mp = 0.06M ⊕ (∼ 0.5R ⊕ , F irr ∼ 740F ⊕ ) for Case 2, which are near the peaks of the detection probability distributions (Figure 13).
ing two stars with M * < 0.8M⊙ at random out of the initial Kepler sample (Kepler 1520b and KOI-2700b) is 2.4%, based on the mass distribution.This is a low probability, but not so low as to justifiably rule out it being a coincidence.If there is a bias towards smaller stars, it could be because of the fact that, at the same planet temperature, a planet around a smaller star is closer to the star and thus more likely to be transiting (Equation 38).Furthermore, the transit depth will be larger due to the lower stellar radius, and the period will be shorter, which may make it easier to identify the systems.These are both effects we did not take into account.Thus, the fact that the systems are around low-mass stars might not be surprising.

SUMMARY
In this work, we have presented a model for the thermal evolution of catastrophically evaporating planets.The model computes heat transport through conduction and convection, including the effects of melting, alongside hydrostatic equilibrium, allowing the mass of the model to be evolved self consistently.We have used this model to show that the catastrophically evaporating planets are likely almost entirely crystallised, other than a shallow molten region on the dayside of the planet.This result is robust to the details of the redistribution of heat from the star, which is challenging to model in full.Therefore, we suggest that the composition of the observed dusty tails samples only the particular region of the mantle that the planet has evaporated to.Future work should include modelling of the chemical evolution of the mantle in order to find the composition of the outer layer at a particular point in time.Additionally, further modelling of the chemical evolution of the lava pool, as well as the outflow itself, may be Table 3. Properties of the observed systems.For Kepler 1520b and KOI-2700b, we take stellar properties from the GKS catalogue (Berger et al. 2020) and orbital periods from Rappaport et al. (2012) and Rappaport et al. (2014) respectively.Properties for K2-22b are taken from Sanchis-Ojeda et al. (2015).Errors on the orbital period are small (∼ 1s).
required to make a full link between the mantle composition and that of the dust.
We also use our model to investigate the occurrence rate of the progenitors of the catastrophically evaporating planets and low-mass planets in general, assuming they have remained at their current locations since birth.When using a simple cut-off for the mass loss rate required to be observable, we obtain an estimate of 0.14% of stars having a planet that will have observable evaporation, in line with the estimates of Perez-Becker & Chiang (2013).However, if we instead consider the theoretical predictions of dust production from Booth et al. (2023), the value is an order of magnitude higher at ∼ 7%.Extending this to general planets in the temperature and mass region close to those which can have detectable mass loss rates, by assuming a power law distribution, the simple cut-off model implies that 0.5%-20% have a planet in the range 0.07-0.58M⊕and 2115-2693 K.This is consistent with observed planet demographics.However, for the dust production model, we find around 1 planet per star in the range 0.02-0.17M⊕and 1775-2367 K.This is up to 100 times more common than Super-Earths and sub-Earths in the same substellar temperature range, derived by Kepler, and is likely inconsistent with the lack of Mercury radius planets observed.This may be evidence that the catastrophically evaporating planets were moved onto their current orbits at later times or imply that the narrow range of temperatures and planet masses that the Booth et al. ( 2023) models predict dust is produced at is too small.However, when considering the detectability of systems in general, we found that the substellar temperatures of the detected systems are better explained by a dust production model, such as that of Booth et al. (2023), than a simple model where planets become observable above a certain mass loss rate.In a simpler model, the hottest planets are most observable as they have the highest mass loss rates.However, in the model of Booth et al. (2023), dust production is reduced at high temperatures.Therefore, our findings motivate further theoretical work to better determine what range of planetary parameters can produce observable dusty outflows.This, along with an improved understanding of any orbital evolution of these systems, will allow us to approach the true occurrence of these systems.

DATA AVAILABILITY
The code and simulated data underlying this article will be shared on reasonable request to the corresponding author.Any observational or experimental data used are available at the references in the main text.
is simply when the viscous and inviscid velocities (Equations 14a and 14b) are equal.
Thus, at any angle with an outward heat flux F (θ) we can work out the temperature-pressure structure using Equations B5 and B6.
The statement that conduction is not important in the inviscid region can be demonstrated as follows.The domination of conduction or convection is governed by the Peclet number where a high Peclet number corresponds to convection dominating heat transport.The second equality comes from considering when the temperature gradient is much higher than adiabatic.
If one considers the ratio of the Peclet and Reynolds numbers in the viscous limit P e Re = νρCP k (B9) one sees that for any reasonable values the Peclet number is larger than the Reynolds number, meaning the fluid first becomes viscous and then conductive.

B2 Fitting of luminosity-temperature function
To aid the convergence of the Henyey scheme it is desirable for the T0(L) relation to be smooth.However, the calculation described in §2.7 can produce numerical variation, particularly in the near-constant T0 section (see Figure 2).We opted to fit a function to our result, motivated by limits to the equations' solutions.
The solutions to Equation B5 have a large temperature increase in the conductive region close to the surface, followed by a less steep, close to adiabatic increase once convection can dominate.The temperature T0 can thus be approximated as where F kρg is the conductive gradient, and Pcrit is the point when conduction dominates.We can define this point using the Peclet number (Equation B8), which we can write as (B12) for the case where Pcrit < P0.For low temperatures, P0 < Pcrit so we take with A1 a constant.Once the system becomes partially melted the exponential part of the viscosity η ∝ exp(−αηϕ) (see Equation 25) becomes the most important so to this region.Around ϕ ≈ ϕc the viscosity, and so the flux, is a very strong function of temperature, meaning that temperature is approximately constant for a large range of fluxes, this region we approximate as a constant Tcrit.
For high temperatures when the fluid is molten η ∼ η l ϕ−ϕc 1−ϕc 2.5 . Since ϕ is linear in T (for a given P , Equation 7) this is not far from a power law in T so we take We combine all these into a master formula of where α, β, γ are smoothing constants and f1, f2 and f3 are given by Equations B13, B15 and B16.We fit the coefficients including the smoothing values using scipy curve_fit.The only value we do not fit is Tcrit which we set to 1778.6.Coefficients for different substellar temperatures are shown in Table B1.
As can be seen, a1 and b1 are very close to our predicted values of 1 and -1.Meanwhile, a3 and b3, show a not insignificant deviation from our predicted values of 0.75 and -0.25.This is not surprising, as we entirely neglected any viscosity dependence.The values of α, β, γ are high, reflecting sharp transitions between regions.
Our fit produces fractional errors on the calculated T0 of ≲ 1.5% which is probably well within errors in the physical parameters and heat transport assumptions.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

θFigure 1 .
Figure1.Schematic of our boundary model ( §2.7).The inner part of the planet is considered spherically symmetric, and its structure is solved for by our Henyey scheme.Its edge is at a pressure P 0 , which has a spherically symmetric temperature, T 0 .The planet's surface has a day-to-nightside temperature gradient, Ts(θ), due to the stellar irradiation determined by Equations 26 and 27.For the outer layer (P < P 0 ), we consider radial fluxes F (θ) such that the flux integrated over the surface is equal to the luminosity for the interior, L.

Figure 2 .
Figure 2. Temperature T 0 at a pressure P 0 = 1 GPa, for boundary conditions with no redistribution and Tss = 2320K, as a function of the mean flux from the planet (luminosity divided by surface area), shown for three different surface gravities g.See §2.7.In this figure, we display the smooth fitted function described in Appendix B2.

2)Figure 3 .
Figure3.Contours of gas mass loss rate for planets of different mass and substellar temperatures and fixed bulk density of 0.67ρ ⊕ , similar to that of Mars, as calculated byBooth et al. (2023).It differs slightly from their Figure5, which is calculated for a fixed mass-radius relation.Dust mass loss rates of 10 −6 , 10 −4 and 10 −2 M ⊕ Gyr −1 are shown with the dotted, dashed and solid black contours.The irregular shape is due to the finite number of data points.

Figure 4 .
Figure4.Evolution of the internal structure of a 0.15 M ⊕ planet with core mass fraction 0.3, and a substellar temperature of 2320K, but no mass loss.Note that the right-hand panels show only the mantle, so the pressure scale is different.The black solid lines on the temperature plot mark the liquidus (upper line) and solidus (lower line), which only apply to the mantle.Time snapshots are roughly spaced logarithmically in time but with extra snapshots between 10 3 − 10 4 yrs to demonstrate the crystallisation of the magma ocean and between 10 9 − 10 10 yrs to show the late time evolution.

Figure 5 .
Figure5.Evolution of the mean surface flux (top panel), temperature (middle panel) and the radius and molten state (bottom panel) of a 0.15 M ⊕ planet with a core mass fraction of 0.3 and a substellar temperature of 2320K, but no mass loss.In the bottom panel, regions to the left of the liquidus are entirely liquid, to the left of the ϕ = ϕc line are partially molten, to the left of the solidus behave like a solid but have some partial melting and to the right of the solidus are fully solid.The region above the 1 GPa line is unresolved by the 1D model and will have melting on the dayside surface but be cold and solid on the night side.We also do not plot any partial melting at the core-mantle boundary, which does occur in our models, at certain times, due to the thermal boundary layer there (see Figure4) because we are more interested in the state towards the surface.
Figure6.Evolution of planets' mass and mass loss rates with different initial masses and substellar temperatures.Mass loss rates are calculated usingBooth et al. (2023).All planets start with a core mass fraction of 0.3.

Figure 7 .
Figure7.Evolution of the molten state of planets with different initial masses and substellar temperatures.Mass loss rates are calculated usingBooth et al. (2023).All planets start with a core mass fraction of 0.3.

Figure 8 .
Figure8.Mean (outward) energy flux and melt fraction at 1 GPa as a function of time for various initial planet masses and substellar temperatures for planets undergoing mass loss according to the models ofBooth et al. (2023).All masses have an initial core mass fraction of 0.3.The increase in melt fraction at late times for some simulations is due to the increase in flux, as discussed in the text.Simulations that end with a dot had full crystallisation up to 1GPa.Those that end with a cross reached the maximum bulk planet density that mass loss rates were computed to.

Figure 9 .Figure 10 .
Figure9.The depth of a dayside magma pool under different assumptions, and defining the boundary of the 1D portion of the code at different pressures, P 0 .For both cases, the planet has a mass of 0.2 M ⊕ , a substellar temperature of Tss = 2320K and undergoes mass loss according toBooth et al. (2023).Note the scale change in the y-axis.In the first case, without circulation (upper panel), a maximum depth is found by calculating the position of the critical melt fraction assuming conduction to the edge of the 1D portion of the code.This naturally depends strongly on the value of P 0 .We also show the depth of P 0 .In the second case (bottom panel), magma pool circulation is included, which produces a much shallower pool.See §4.1 for further details of the calculations.

Figure 12 .
Figure 12.The same as Figure 11, but instead planets are considered observable if the models of Booth et al. (2023) predict dust production.Note the difference in the range of initial masses compared to Figure 11.

Figure 13 .
Figure 13.The probability that a planet of a given initial mass would be observable around a star, marginalised over the stellar mass and age distribution of stars observed in the Kepler primary mission, P detect (Mp, Tss) in Equation 40.In the top panel, we consider a planet observable if its mass loss rate exceeds 0.1 M ⊕ Gyr −1 (Perez-Becker & Chiang 2013), Case 1 in the text.In the bottom panel, planets are considered detectable if the models of Booth et al. (2023) predict dust production (see Figure 3, Case 2).Note the different y-axis scales and range of P detect values.Solid lines show the substellar temperatures of the known systems, with associated error bars shown with dashed lines.Note that the best fit substellar temperature of K2-22b coincides with the upper error bar of KOI 2700b and that K2-22b has significantly higher errors, as the stellar parameters are less well constrained.

Figure 14 .
Figure 14.Planetary occurrence rate as a function of α and β (the power law indices of Mp and Tss in Equation41) for Case 1 and 2. We define our occurrence rate in the same way asPetigura et al. (2022) as 0.5 ∂Np/∂ log F irr for a 0.23 dex bin in radius (see text).Np is the number of planets per star and F irr the stellar irradiation at the planet's orbital radius.Note that for Case 1 the scale is logarithmic, while for Case 2 it is not.Our measure is evaluated at Tss = 2400 K, Mp = 0.2M ⊕ for Case 1 (∼ 0.7R ⊕ , F irr ∼ 1400F ⊕ ) and Tss = 2050 K, Mp = 0.06M ⊕ (∼ 0.5R ⊕ , F irr ∼ 740F ⊕ ) for Case 2, which are near the peaks of the detection probability distributions (Figure13).

Figure 15 .
Figure15.Comparison of our occurrence rate with those measured for super-Earths byPetigura et al. (2022) and sub-Earths byHsu et al. (2019), and that inferred for the catastrophically evaporating planets byPerez-Becker & Chiang (2013).The y-axis is the definition of occurrence rate used byPetigura et al. (2022), integrated over 0.23 dex in planet radius.We scale the occurrence rate found byPerez-Becker & Chiang (2013) according to the area of substellar temperature-initial mass parameter space where mass loss should be observable (P detect > 10 −4 ) for our Case 1 (see Figure13).Horizontal error bars approximate this flux range and are centred on the temperature that Perez-Becker & Chiang (2013) consider (2145 K), whereas for our values we use temperatures close to the peak of P detect .Note that while the logarithmic range of radius/mass is the same for all measurements, they are centred at different values.The mass range forPetigura et al. (2022) is ∼ 1M ⊕ -8M ⊕ , forHsu et al. (2019) it is ∼ 0.6M ⊕ -0.3M ⊕ and for ours it is ∼ 0.07-0.58M ⊕ for Case 1 and ∼ 0.02-0.17M⊕ for Case 2. The vertical error bars associated with our measurement are due to the range in considered power-law indices of the number distribution, as seen in Figure14.The vertical error bars on the measurement from Perez-Becker & Chiang (2013) are due to the uncertainty they incorporate.

Table 2 .
Best fit parameters for our boundary condition function Equation39

Table B1 .
Best fit parameters for our boundary condition function Equation B17.