Evaporation of Close-in Sub-Neptunes by Cooling White Dwarfs

Motivated by the recent surge in interest concerning white dwarf (WD) planets, this work presents the first numerical exploration of WD-driven atmospheric escape, whereby the high-energy radiation from a hot/young WD can trigger the outflow of the hydrogen–helium envelope for close-in planets. As a pilot investigation, we focus on two specific cases: a gas giant and a sub-Neptune-sized planet, both orbiting a rapidly cooling WD with mass M * = 0.6 M ⊙ and separation a = 0.02 au. In both cases, the ensuing mass outflow rates exceed 1014 g s−1 for WD temperatures greater than T WD ≳ 50,000 K. At T WD ≃ 18,000 K (22,000 K), the sub-Neptune (gas giant) mass outflow rate approaches 1012 g s−1, i.e., comparable to the strongest outflows expected from close-in planets around late main-sequence stars. Whereas the gas giant remains virtually unaffected from an evolutionary standpoint, atmospheric escape may have sizable effects for the sub-Neptune, depending on its dynamical history, e.g., assuming that the hydrogen–helium envelope makes up 1% (4%) of the planet mass, the entire envelope would be evaporated away so long as the planet reaches 0.02 au within the first 230 Myr (130 Myr) of the WD formation. We discuss how these results can be generalized to eccentric orbits with effective semimajor axis a′=a/(1−e2)1/4 , which receive the same orbit-averaged irradiation. Extended to a much broader parameter space, this approach can be exploited to model the expected demographics of WD planets as a function of their initial mass, composition, and migration history, as well as their potential for habitability.


INTRODUCTION
The systematic search for Earth-like biosignatures in habitable zone planets around G-type stars is currently prohibitive for the next two decades.In spite of their high activity levels, nearby low-mass stars have thus emerged as favorable targets.More recently, the discovery space has expanded to include planets around dead stars, such as white dwarfs (WDs).To date, a handful of exoplanet candidates have been discovered around WDs, including one within ∼0.02 AU (Vanderburg et al. 2020; see also Veras 2021 and references therein).Additionally, between 25-20% of WDs exhibit atmospheric metal contamination, likely from in-falling planetary debris (Farihi et al. 2009;Klein et al. 2010;Zuckerman et al. 2010;Koester et al. 2014).
With their small radii and low luminosities, WDs are challenging targets for planet searches with traditional techniques.However, the large ratio of planet to WD radius makes any such system a relatively inexpensive target for follow-up transit spectroscopy.Indeed, Kaltenegger et al. (2020) articulated the notion of a "window of opportunity" for detecting biosignatures around WD habitable planets, should they exist.More recently, Limbach et al. (2022) demonstrated how the near-featureless spectra of nearby WDs can be exploited to detect an infrared excess from cold Jovian exoplanets with modest exposures.Systematic searches for cold Jovian and (transiting) terrestrial planets around nearby WDs are currently under way with JWST and TESS, respectively.
At the same time, a robust theoretical framework for the formation and evolution of WD planets is still lacking.First generation systems beyond ∼5 AU may survive the red giant phase of the host star; some planets (or perhaps their moons) could then migrate inward via planet-planet scattering (Veras 2016; Chamandy et al. 2021), Kozai-Lidov oscillations (Stephan et al. 2017), or common envelope evolution (Lagos et al. 2021).Alternatively, second generation systems could form after the WD itself (Debes & Sigurdsson 2002;Zuckerman et al. 2010;van Lieshout et al. 2018), although the abil-ity of any such planet to actually develop and/or retain the conditions for life under these circumstances is far from obvious (Agol 2011;Fossati et al. 2012;Barnes & Heller 2013;Kozakis et al. 2020).
Absorption of high-energy stellar radiation in the upper atmosphere can significantly affect the evolution of close-in exoplanets.X-ray (5-100 Å) and EUV (100-912 Å) photonscollectively termed XUV-are absorbed at high altitude, in the lower density parts of the atmosphere, where they dissociate and/or ionize both molecules and atoms.The resulting free electrons heat the gas collisionally.These processes inflate the gas and drive hydrodynamic outflows that can potentially remove the primordial atmospheric envelope of hydrogen and helium, a process known as XUV-driven atmospheric escape (Vidal-Madjar et al. 2003;Yelle 2004;Tian et al. 2005;Koskinen et al. 2014).
Along with core-powered mass loss (Ginzburg et al. 2018;Gupta & Schlichting 2019), photoionization-driven atmospheric escape is thought to play a significant role in shaping the observed properties of the Kepler planet population (Fulton et al. 2017), likely contributing to carving the observed radius valley by stripping sub-Neptune-sized planets of their primordial envelopes, and turning them into rocky cores (e.g., Owen & Wu 2013, 2017;Rogers & Owen 2021;Affolter et al. 2023;Owen & Schlichting 2023).
Because of the vastly different regime of parameter space realized by WDs, in particular the different XUV/bolometric ratios and radii compared to FGK-type stars, which dominate the Kepler sample, the photoionizing radiation arising from young WDs is expected to affect the atmosphere of close-in planets in fundamentally different ways.Depending on when and where WD planets form and/or migrate inward, WD-driven outflows may yield different exoplanet properties compared to those of Kepler stars.The WDs can drive far more efficient mass loss rates (compared to young stars) at very early stages, and much less (to virtually none) after the first few hundred million years.This sequence of events affects both the expected planet radius distribution and their potential for habitability.Shedding any previously accreted hydrogen-helium envelopes is also a prerequisite for habitability, as the presence of a massive, light element envelope dramatically increases surface temperature and pressure, making it impossible for liquid water to exist, regardless of whether the planet resides in the habitable zone.
Beyond a simple energy-limited approach (e.g., Watson et al. 1981; see also below), the range of parameters pertaining to WD-driven atmospheric outflows is completely unexplored in the literature (see Schreiber et al. 2019, e.g., on the ultimate fate of the Solar System gas giants using energy-limited atmospheric escape).Addressing these ques-tions quantitatively requires an investigation of the physics of XUV-driven photoevaporation for the specific case of a WD spectrum.This Paper takes the first step in this direction by considering the idealized case of a sub-Neptune-sized planet and a gas giant in a circular orbit with radius a =0.02 AU around a rapidly cooling WD.We also show how these results can be generalized to planets on eccentric orbits subject to the same (orbit-averaged) irradiation.

FROM STELLAR TO WD-DRIVEN ATMOSPHERIC ESCAPE
Assuming that all of the stellar XUV flux that is absorbed by a planet atmosphere is converted into expansion work, the instantaneous ion mass loss rate can be expected to scale linearly with the ratio between the incident XUV flux at the planet's orbital separation, F XUV , and the average planetary mass density, ρ p (Erkaev et al. 2013).In practice, numerical work has shown conclusively that the validity of this approximation (known as energy-limited escape) is relatively narrow, as radiative losses can seldom be neglected (e.g.Murray-Clay et al. 2009;Owen & Jackson 2012;Erkaev et al. 2016;Salz et al. 2016;Salz et al. 2016;Caldiroli et al. 2022).This complication warrants a numerical treatment of the problem.

Numerical Approach
Our starting point is the publicly available code ATmospheric EScape1 (ATES).For a given planetary system architecture (planet mass, radius, equilibrium temperature, average orbital separation, stellar mass and input stellar spectrum) ATES computes the temperature, density, velocity, and ionization fraction profiles of a highly irradiated H-He atmosphere, and calculates the instantaneous, steady-state massloss rate.This computation is done by solving the monodimensional Euler, mass, and energy conservation equations in radial coordinates through a finite-volume scheme; the expression for the gravitational potential also accounts for the Roche potential (Erkaev et al. 2007), which may significantly affect the inferred mass loss rates.The hydrodynamics module is paired with a photoionization equilibrium solver that includes cooling via bremsstrahlung, recombination, and collisional excitation and ionization, whilst also accounting for advection of the different ion species.ATES has been validated through a systematic comparison against The PLUTO-CLOUDY Interface (TPCI; Salz et al. 2015), a publicly available interface between the magnetohydrodynamics code PLUTO (Mignone et al. 2012) and the plasma simulation and spectral synthesis code CLOUDY Figure 1.Left: evaporation efficiency (squares, left y-axis) and atmospheric mass loss rate (circles, right y-axis) as a function of XUV irradiation (bottom x-axis) and or WD age (top x-axis) for the case of a sub-Neptune-sized planet (8 M ⊕ , 2.7 R ⊕ ) orbiting a 0.61 M ⊙ WD at 0.02 AU.The symbol colors trace the WD cooling history (top axis), with T WD varying between 69,000 K (deepest red) and 5,000 K (deepest blue).Right: same as left, but for the case of a gas giant (1.24 M J and 1.19 R J ). (Ferland et al. 1998) (the photoionization equilibrium equations are solved "on the fly" in ATES, which results in a factor 80 speed-up of the computing time compared to TPCI).In its current, publicly available version, ATES recovers stable, steady-state solutions for ϕ p < ∼ 12.9 + 0.17 log F XUV (in cgs units), where ϕ p is the planet's gravitational potential energy.The code is well suited to model photoionization-driven atmospheric escape from planets (i) with radii in excess of > ∼ 1.7R ⊕ (i.e., sub-Neptune-sized, or larger), for which a sizable hydrogen-helium envelope may be retained (Lopez & Fortney 2014a), and (ii) at distances within ∼0.5 A.U., beyond which atmospheric mass loss is typically negligible (Jeans escape limit), and within which the upper thermosphere comprises only atomic (rather than molecular) hydrogen.
Hereafter, we define the evaporation efficiency2 η as the ratio between the mass outflow rate inferred numerically ( Ṁ) and the well-known energy-limited approximation for the case where all of the absorbed XUV flux is converted into abiabatic expansion: where K accounts for the host star gravitational pull (Erkaev et al. 2007).The evaporation efficiency is set primarily by the planet's gravity.Specifically, atmospheric mass loss is always far less efficient (η ≪ 1) than the energy-limited benchmark for gas giants (more specifically, for high-gravity planets, with ϕ p > ∼ 14 × 10 12 erg g −1 ).As demonstrated in Section 3 of Caldiroli et al. (2022), this inefficiency arises because, after accounting for photo-ionization losses, the mean energy that remains available to heat the ions in the atmosphere is lower than the ions' gravitational binding energy.For a fixed value of ϕ p , the evaporation efficiency of highgravity planets increases with XUV irradiation, as the fractional contribution of adiabatic cooling increases with flux at a faster pace than Lyα losses, which are the dominant radiative cooling channel.
Conversely, for lower gravity planets (i.e., Neptune-sized and smaller, with ϕ p < ∼ 8 × 10 12 erg g −1 ) the mass loss rate is largely independent of ϕ p , and is regulated primarily by F XUV .These planets can be expected to exhibit energylimited outflows only below a critical irradiation level (η ≃ 1 for F XUV < ∼ 10 5 erg cm −2 s −1 ).As the irradiation increases, however, the evaporation efficiency decreases, by up to an order of magnitude (even though the mass loss rate increases).This trend is due to a progressive increase in the fractional contribution of advective cooling, followed by a sharp surge in Lyα cooling, with respect to purely adiabatic losses.
This framework enables a reliable, physically motivated prediction of the photoionization-driven mass outflow rate For the latter, the non-monotonic behavior of the abiabatic cooling is responsible for the non monotonic trend of the evaporation efficiency shown in the right panel of Figure 1.
for a given known or theoretical planetary system in close orbit around a late main sequence star (whilst foregoing other mechanisms for driving atmospheric outflows, such as corepowered mass loss and/or boil-off).We note that the above results are fairly insensitive to the assumed X-ray and EUV spectral shape, so long as the combination of shape and normalization preserve the total number of photons per unit time in each band.
We aim to replicate a similar investigation for the parameter space pertaining to planets orbiting a hot WD.To model a cooling WD, we make use of BaSTI3 (Salaris et al. 2022) stellar evolution models.These consider Carbon-Oxygen WDs with masses of 0.54, 0.61, 0.68, 0.77, 0.87, 1.0, and 1.1 M ⊙ , a range of metallicities4 , hydrogen and helium or pure helium atmospheres.Throughout, we consider a 0.61 M ⊙ WD with Z=0.017, and mixed atmosphere.
Unlike for late main sequence stars, for high temperature/young ages the bulk of the WD flux is emitted in the EUV portion of the spectrum, that is, the energy range that is primarily responsible for driving hydrodynamic outflows.Such extreme ionizing fluxes can be expected to have dramatic consequences on the outflow properties.The ensuing numerical challenges are two-fold.Firstly, we can expect the development of an extremely sharp ionization front close to the inner boundary of the computational domain.Underresolved ionization fronts translate into sharp discontinuities in the heating rate, which cannot be handled by the code hydrodynamic solver.In addition, non-physical shock waves may develop if the initial and boundary conditions are not carefully chosen.Specifically, a moving shock can arise, e.g., if too much energy is deposited into the initial state.These challenges compromise the stability of ATES' numerical scheme.
To address them, we implement an adaptive grid refinement routine within the domain of the simulation; this modification enables us to handle sharp, stationary ionization fronts.Next, we modify the code so that the atmosphere is irradiated gradually over the first few thousands time steps; this feature avoids the formation of strong shock waves.

Circular Orbit Results
The modifications and additions outlined above enable the code reach steady-state solutions up to F XUV ≃ 10 11 erg cm −2 s −1 , which is about five orders of magnitude higher than the highest XUV irradiation considered by Caldiroli et al. (2022) for input stellar spectra (albeit assuming different planets' equilibrium temperatures).For this pilot study we examine a low-gravity pl anet (8 M ⊕ and 2.7 R ⊕ , i.e., a typical sub-Neptune) and a high-gravity planet (1.24 M J and 1.19 R J , i.e., a typical gas giant), both on a circular orbit at 0.02 AU of a 0.61 M ⊙ WD, where the orbital separation is chosen specifically to match WD 1856+534 b's ( Vanderburg et al. 2020).
The main results are summarized in Figure 1.The initial WD surface temperature is 69,000 K (dark red); the corresponding planet equilibrium temperature is T eq ≃ 3, 200 K.After 600 Myr, the WD has cooled down to 10,000 K (darkest blue), yielding T eq ≃ 380 K.For the sub-Neptune case (left panel of Figure 1), the mass loss rate exceeds 10 12 g s −1 for T WD > ∼ 18,000 K, i.e., within the first ∼100 Myr.Mass outflow rates as high as 10 9−10 g s −1 are still driven after 600 Myr, when T WD < ∼ 10, 000 K. The evaporation efficiency increases steeply, and monotonically, with decreasing F XUV , i.e., as the WD cools down.Nevertheless, the outflow is never quite energy-limited (i.e., η < ∼ 1 at all times).As shown in the left panel of Figure 2, this behavior can be understood in terms of the fractional increase in the adiabatic with respect to the (total) radiative cooling rate going towards lower F XUV .
The gas giant (right panel of Figure 1) shows very low evaporation efficiencies, with η < 10 −2 at all WD ages/temperatures.Nevertheless, the ensuing mass outflow rates range between 10 8 and 10 15 g s −1 , with a peak approaching 5 M ⊕ /Myr.The evaporation efficiency increases slightly overt the first ∼30 Myr, as F XUV decreases from 10 11 down to ≃ 10 8 erg s −1 cm −2 .Only at later times/lower F XUV does η decrease with irradiation.This behavior can be understood again by examining the relative contribution of adiabatic vs. (total) radiative cooling, shown in the right panel Figure 2.For XUV irradiations below about < ∼ 10 8 erg s −1 cm −2 , the ratio of abiabatic to radiative cooling decreases with decreasing flux, leading to a decrease in efficiency.However, the ratio of abiabatic to radiative cooling increases with decreasing flux between 10 11 and 10 8 erg s −1 cm −2 , which means that the evaporation efficiency increases with decreasing flux over this range.
These results can be converted into approximate atmospheric evaporation timescales.Ignoring momentarily how changes to the planet's atmospheric mass fraction affect its radius, the instantaneous mass loss rates in Figure 1 can be used to estimate how long it would take for a close-in planet which has retained its light element envelope to shed it away.Figure 3 illustrates the fractional amount of mass that the sub-Neptune (left panel) and gas giant (right panel) planet would lose to hydrodynamic escape as a function of (i) the planet "arrival time" t start at given WD temperature (T WD,start ), and (ii) its "residence" time (t − t start ), respectively.
As for stellar irradiation, hydrodynamic escape will have sizable evolutionary effects only for the case of a close-in (0.02 AU) sub-Neptune-sized planet.Specifically, the dashed black lines in Figure 3 track the locus of the points where the planet would shed 0.1, 1, 2, . . . up to 10% of its mass (in the form of a hydrogen-helium envelope) via hydrodynamic escape.As an example, assuming that it makes up 1 [/4] per cent of the planet total mass, the entire envelope of the sub-Neptune would be evaporated in less than about 400 Myr so long as the planet arrives at 0.02 AU within the first 230 [/130] Myr of the WD system life.Note that, although the plot extends to nearly complete evaporation (i.e., ∼80% M p ), the hydrogen-helium envelope of realistic sub-Neptune-sized planets likely does not exceed mass fractions of a few percent (Lopez & Fortney 2014b).
Fractional mass losses are far less extreme for the gas giant; this system will only shed close to 1% of its initial mass if it starts to be exposed to the extreme XUV flux of an extremely young WD, i.e. within the first few Myr.
In summary, strong irradiation from a very young/hot WDs can be expected to drive significantly greater atmospheric mass outflow rates than a late main sequence stars.For the sub-Neptune-sized planet WD-driven atmospheric escape can shed the entire planet's light element envelope only if the planet reaches an orbital separation as close as 0.02 AU within the first few hundreds Mys of the WD formation, which in turn requires a very efficient mechanism for shrinking the orbit to sub-AU separations.The gas giant (at the same orbital separation) remains fairly unaffected, regardless of how fast it reaches such a very close orbit.

Generalization to Eccentric Orbits
So far, we have explored the idealized case of a sub-Neptune and a gas giant on a circular orbit, at 0.02 AU of a rapidly cooling WD, and shown that only in the former case is WD-driven atmospheric escape not negligible.We thus focus on the sub-Neptune case hereafter.There are several reasons why extending these calculations to the case of non-zero eccentricity might be highly desirable, starting with the fact that numerical work indicates that scattering events, leading to highly eccentric orbits, are likely the most efficient mechanism to drive this kind of planets within sub-AU periastron distances of young WDs (Veras 2016).
Whereas we defer a full exploration of the parameter space to a follow-up work, below we examine how the results presented above can be extended to orbits of arbitrary eccentricity e 0. To do so, we wish to derive an analytical expression that characterizes eccentric orbits with the same orbitaveraged irradiation as the circular orbit case explored above (i.e., e = 0 and a = 0.02 AU).We start by deriving the analytic expression of the flux ⟨F⟩ T , averaged over the orbital period T , that is received by a planet orbiting a WD (or a star) of luminosity L and mass M * ≫ M p : where r(t) is the time-dependent orbital separation, and r p and r a are the (t = 0) pericenter and (t = T/2) apocenter distance, respectively.The radial component of the orbital velocity dr/dt in Equation 2 can be expressed in terms of the semi-major axis a and of the orbital eccentricity e, through the specific energy of the orbit, E: where the tangential component of the orbital velocity is written in terms of the orbit's specific angular momentum L. Since E = −GM * /(2a) and L = GM * a(1 − e 2 ), tion (3) yields: Substituting for dr/dt above in Equation ( 2), with T = 2π a 3 /(GM * ), gives: where the limits written as p = a(1 − e) and r a = + e).Integrating Equation ( 5) by substitution (with y ≡ a/r) we find: consistent with the result of Adams & Laughlin (2006).The above derivation shows that, over a single eccentric orbit, a planet experiences an average irradiation which is a factor (1 − e 2 ) −1/2 higher compared to the case of a circular orbit with the same energy.A different way to express this result is that the same orbit-averaged irradiation that is experienced by a planet on a circular orbit of radius a will be experienced by a planet on a (lower energy) eccentric orbit with semimajor axis a ′ = a(1 − e 2 ) −1/4 .In order to translate an orbit-averaged flux into an orbitaveraged atmospheric mass loss rate we need to account for the functional dependence of the instantaneous mass loss rate on irradiation.If the scaling is linear (i.e., energy-limited escape), then the orbit-averaged mass loss rate coincides with that of a circular orbit with the same orbit-averaged flux.The scaling of Ṁ with irradiation becomes shallower at higher irradiation, as a progressively higher fraction of the absorbed WD radiation is used up in radiative processes as opposed to adiabatic expansion.In the limiting case where Ṁ scales as ∝ F 1/2 XUV (i.e., radiation-limited escape; Murray-Clay et al. 2009), it can be shown analytically (following similar steps as in Equations [1-6]) that the orbit-averaged mass loss rate is a factor (1 − e 2 ) 1/4 lower than the mass loss rate for a circular orbit with the same orbit-averaged flux.
The intermediate regime between energy-and radiationlimited warrants a numerical treatment.For three specific choices of eccentricity and semi-major axes (i.e., e = 0.9, 0.95, 0.98, and a ′ = 0.03, 0.036, 0.044 AU, respectively5 , such that the orbit-averaged flux equals that of a circular a = 0.02 AU orbit), we calculate r(t) by integrating Equation (4).From this, we derive the time-dependent irradiation F XUV (t) along the orbit.To obtain the time-dependent mass loss rate, and its average value, we interpolate the Ṁ(F XUV ) results in the left panel of Figure 1 through a cubic spline function.This confirms that progressively increasing the eccentricity (while reducing the orbit's energy) reduces the average mass loss rate along a single orbit, albeit the difference is a factor of about 2 at most.As expected, the differences are more pronounced at high irradiation.This finding is con-sistent with the analytic results above, whereby the circular to eccentric mass loss rate ratio is < ∼ (1 − e 2 ) −1/4 .We emphasize that the above estimates can only be taken as indicative, as they ignore the fact that the size of the system's Roche lobe varies with eccentricity.

SUMMARY AND FUTURE WORK
This work presents the first numerical investigation of planetary atmospheric escape driven by strong irradiation from a hot/young WD.We have examined the specific case of a sub-Neptune planet and a gas giant on a 0.02 AU circular orbit around a rapidly cooling 0.6 M ⊙ WD, and showed how the ensuing mass loss rates and evaporation time-scales can be extended to lower energy, eccentric orbits with the same orbit-averaged irradiation.In a forthcoming paper (Caldiroli et al., in prep.)we plan to vastly expand the parameter space, and to carry out proper evolutionary models, assuming different core compositions, envelope mass fractions, and using theoretical mass-radius relations to update the planetary parameters at each time step.
As shown by Adams (2011) and Owen & Adams (2014), any appreciable planetary magnetic field (e.g., 1 G or stronger) is likely to exert a greater magnetic pressure than the ram pressure of the hydrodynamic outflow over a large portion of the parameter space; this acts to limit the net mass outflow rate.Extending this kind of investigation to close-in planets around very young/hot WDs is necessary to attain a more realistic estimate of their cumulative mass loss.
Notwithstanding the limitations of this approach (we do not know when/where WD planets formed), the results can be reversed to place constraints on the residence time of a given planet at a given orbital distance, and to argue in favor or against the need for very efficient migration.The model presented here can be generalized to new WD planet candidates, or used to derive a theoretical maximum (present-day) WD planet mass as function of the parent mass for different evolutionary histories.For the specific case of planets within the outer edge of the "continuously habitable zone" (Agol 2011), it can be used to estimate the relevant timescales over which any primordial/residual hydrogen-helium envelope can be photo-evaporated away-a prerequisite for habitability.Whether or not WDs indeed host a sizable planet population, the observational community has already chosen to invest massively in this opportunity.As with M dwarf hosts, virtually every candidate planet is going to be follow-up for spectroscopy.Parallel efforts on the modeling side are necessary to ensure a physically motivated interpretation of WD planets, inform future observational campaigns, and attempt to predict/model their demographics.

Figure 2 .
Figure2.Fractional contributions of the total radiative cooling rate (solid blue line) and abiabatic cooling rate (solid green line) with respect to the heating rate (solid red line), shown as a function of the XUV irradiation/WD temperature for the sub-Neptune (left) and gas giant (right).For the latter, the non-monotonic behavior of the abiabatic cooling is responsible for the non monotonic trend of the evaporation efficiency shown in the right panel of Figure1.

Figure 3 .
Figure 3. Left: Expected evaporation timescales for the sub-Neptune-sized simulated in the left panel of Figure 1.The color bar on the right shows the fractional amount of mass lost to hydrodynamic escape as a function of the planet arrival time at 0.02 AU (t start ), and of the total residence time at this separation (t − t start The dashed black lines track the locus of the points planet would shed 1, to 10% its the form of a H-He envelope-via hydrodynamic escape.Same as left, but for the gas giant case in the right panel of Figure 1.