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

Our starting point is the publicly available code Atmospheric Escape ⁹ (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 mass-loss 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 while also accounting for advection of the different ion species. ATES has been validated through a systematic comparison against The PLUTO-CLOUDY Interface (TPCI; Salz & Mignone 2015), a publicly available interface between the magnetohydrodynamic code PLUTO (Mignone et al. 2012) and the plasma simulation and spectral synthesis code CLOUDY (Ferland et al. 1998; the photoionization equilibrium equations are solved "on the fly" in ATES, which results in a factor 80 speedup of the computing time compared to TPCI). In its current, publicly available version, ATES recovers stable, steady-state solutions for ${\phi }_{p}\lesssim 12.9+0.17\mathrm{log}{F}_{\mathrm{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.7 R_⊕ (i.e., sub-Neptune sized or larger), for which a sizable H–He envelope may be retained (Lopez & Fortney 2014), and (ii) at distances within ∼0.5 au, 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 efficiency ¹⁰ η as the ratio between the mass outflow rate inferred numerically ($\dot{M}$) and the well-known energy-limited approximation for the case where all of the absorbed XUV flux is converted into adiabatic expansion:

$$\begin{eqnarray}&&\eta =\dot{M}{\left(\displaystyle \frac{3{F}_{\mathrm{XUV}}}{4{GK}{\rho }_{p}}\right)}^{-1},\end{eqnarray} \tag{ 1 }$$

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¹² erg g⁻¹). As demonstrated in Section 5 of Caldiroli (2022), this inefficiency arises because, after accounting for photoionization 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 high-gravity 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¹² erg g⁻¹) the mass-loss rate is largely independent of ϕ_p and is regulated primarily by F_XUV. These planets can be expected to exhibit energy-limited outflows only below a critical irradiation level (η ≃ 1 for F_XUV ≲ 10⁵ erg cm⁻² s⁻¹). 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 a given known or theoretical planetary system in close orbit around a late main-sequence star (while foregoing other mechanisms for driving atmospheric outflows, such as core-powered 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 BaSTI ¹¹ (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 metallicities; ¹² and 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 twofold. First, 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, nonphysical 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's 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 thousand time steps; this feature avoids the formation of strong shock waves.

The modifications and additions outlined in Section 2.1 enable the code to reach steady-state solutions up to F_XUV ≃ 10¹¹ erg cm⁻² s⁻¹, which is about 5 orders of magnitude higher than the highest XUV irradiation considered by Caldiroli (2022) for input stellar spectra (albeit assuming different planets' equilibrium temperatures). For this pilot study we examine a low-gravity planet (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 ≃ 3200 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¹² g s⁻¹ for T_WD ≳ 18, 000 K, i.e., within the first ∼100 Myr. Mass outflow rates as high as 10^9–10 g s⁻¹ 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 toward lower F_XUV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The gas giant (right panel of Figure 1) shows very low evaporation efficiencies, with η < 10⁻² at all WD ages/temperatures. Nevertheless, the ensuing mass outflow rates range between 10⁸ and 10¹⁵ g s⁻¹ , with a peak approaching 5 M_⊕ Myr⁻¹. The evaporation efficiency increases slightly over the first ∼30 Myr, as F_XUV decreases from 10¹¹ down to ≃10⁸ erg s⁻¹ cm⁻². Only at later times/lower F_XUV does η decrease with irradiation. This behavior can be understood again by examining the relative contribution of adiabatic versus (total) radiative cooling, shown in the right panel Figure 2. For XUV irradiations below about ≲10⁸ erg s⁻¹ cm⁻², the ratio of adiabatic to radiative cooling decreases with decreasing flux, leading to a decrease in efficiency. However, the ratio of adiabatic to radiative cooling increases with decreasing flux between 10¹¹ and 10⁸ erg s⁻¹ cm⁻², 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 that 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 a given WD temperature (T_{WD, start}) and (ii) its "residence" time (t − t_start), respectively.

Zoom In Zoom Out Reset image size

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 H–He envelope) via hydrodynamic escape. As an example, assuming that it makes up 1% (4%) 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 Myr (130 Myr) of the WD system life. Note that, although the plot extends to nearly complete evaporation (i.e., ∼80% M_p ), the H–He envelope of realistic sub-Neptune-sized planets likely does not exceed mass fractions of a few percent (Lopez & Fortney 2014).

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 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 of Myr 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.

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 nonzero 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 planet 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 in Section 2.2 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 orbit-averaged 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, which is received by a planet orbiting a WD (or a star) of luminosity L and mass M_* ≫ M_p :

$$\begin{eqnarray}&&\langle F{\rangle }_{T}=\displaystyle \frac{2}{T}{\int }_{0}^{T/2}{dt}\,\displaystyle \frac{L}{4\pi {r}^{2}}=\displaystyle \frac{L}{2\pi T}{\int }_{{r}_{p}}^{{r}_{a}}\displaystyle \frac{{dt}}{{dr}}\displaystyle \frac{{dr}}{{r}^{2}},\end{eqnarray} \tag{ 2 }$$

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 semimajor axis a and of the orbital eccentricity e, through the specific energy of the orbit, ${ \mathcal E }$:

$$\begin{eqnarray}&&{ \mathcal E }=\displaystyle \frac{1}{2}\left[{\left(\displaystyle \frac{{dr}}{{dt}}\right)}^{2}+{\left(\displaystyle \frac{{ \mathcal L }}{r}\right)}^{2}\right]-\displaystyle \frac{{{GM}}_{* }}{r},\end{eqnarray} \tag{ 3 }$$

where the tangential component of the orbital velocity is written in terms of the orbit's specific angular momentum ${ \mathcal L }$. Since ${ \mathcal E }=-{{GM}}_{* }/(2a)$ and ${ \mathcal L }=\sqrt{{{GM}}_{* }a(1-{e}^{2})}$, Equation (3) yields

$$\begin{eqnarray}&&\displaystyle \frac{{dr}}{{dt}}=\sqrt{\displaystyle \frac{{{GM}}_{* }}{a}}\sqrt{-1+2\displaystyle \frac{a}{r}-\displaystyle \frac{{a}^{2}}{{r}^{2}}(1-{e}^{2})}.\end{eqnarray} \tag{ 4 }$$

Substituting for dr/dt above in Equation (2), with $T=2\pi \sqrt{{a}^{3}/({{GM}}_{* })}$, gives

$$\begin{eqnarray}&&\langle F{\rangle }_{T}=\displaystyle \frac{L}{4{\pi }^{2}a}{\int }_{a(1-e)}^{a(1+e)}\displaystyle \frac{{dr}}{r}\displaystyle \frac{1}{\sqrt{-{r}^{2}+2{ar}-{a}^{2}(1-{e}^{2})}},\end{eqnarray} \tag{ 5 }$$

where the integration limits are written as r_p = a(1 − e) and r_a = a(1 + e). Integrating Equation (5) by substitution (with y ≡ a/r) we find

$$\begin{eqnarray}&&\langle F{\rangle }_{T}=\displaystyle \frac{L}{4{\pi }^{2}{a}^{2}}{\int }_{1/(1+e)}^{1/(1-e)}{dy}\,\displaystyle \frac{1}{\sqrt{-1+2y-{y}^{2}(1-{e}^{2})}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{L}{4\pi {a}^{2}\sqrt{1-{e}^{2}}},\end{eqnarray} \tag{ 7 }$$

consistent with the result of Adams & Laughlin (2006). The above derivation shows that, over a single eccentric orbit, a planet experiences an average irradiation that 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^{\prime} =a{(1-{e}^{2})}^{-1/4}$.

In order to translate an orbit-averaged flux into an orbit-averaged 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 $\dot{M}$ 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 $\dot{M}$ scales as $\propto {F}_{\mathrm{XUV}}^{1/2}$ (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 radiation-limited warrants a numerical treatment. For three specific choices of eccentricity and semimajor axes (i.e., e = 0.9, 0.95, 0.98, and $a^{\prime} =0.03,0.036,0.044\,\mathrm{au}$, respectively, ¹³ 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 $\dot{M}({F}_{\mathrm{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 consistent with the analytic results above, whereby the circular-to-eccentric-mass-loss-rate ratio is $\lesssim {(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.