Evolution of atmospheric escape in close-in giant planets and their associated Ly$\alpha$ and H$\alpha$ transit predictions

TL;DR

This study models atmospheric escape in close-in giant exoplanets, predicting Ly-alpha and H-alpha transit signatures over their evolution, revealing how mass loss impacts planetary atmospheres and observable signals.

Contribution

It introduces a hydrodynamic escape model combined with spectroscopic transit simulations to analyze atmospheric evolution and escape signatures in giant exoplanets.

Findings

Younger planets exhibit higher escape rates due to increased irradiation.

Less massive planets lose significantly more mass over time.

Ly-alpha absorption remains significant at older ages, H-alpha diminishes after 1.2 Gyr.

Abstract

Strong atmospheric escape has been detected in several close-in exoplanets. As these planets consist mostly of hydrogen, observations in hydrogen lines, such as Ly-alpha and H-alpha, are powerful diagnostics of escape. Here, we simulate the evolution of atmospheric escape of close-in giant planets and calculate their associated Ly-alpha and H-alpha transits. We use a one-dimensional hydrodynamic escape model to compute physical properties of the atmosphere and a ray-tracing technique to simulate spectroscopic transits. We consider giant (0.3 and 1M_jup) planets orbiting a solar-like star at 0.045au, evolving from 10 to 5000 Myr. We find that younger giants show higher rates of escape, owing to a favourable combination of higher irradiation fluxes and weaker gravities. Less massive planets show higher escape rates (1e10 -- 1e13 g/s) than those more massive (1e9 -- 1e12 g/s) over their evolution. We estimate that the 1-M_jup planet would lose at most 1% of its initial mass due to escape, while the 0.3-M_jup planet, could lose up to 20%. This supports the idea that the Neptunian desert has been formed due to significant mass loss in low-gravity planets. At younger ages, we find that the mid-transit Ly-alpha line is saturated at line centre, while H-alpha exhibits transit depths of at most 3 -- 4% in excess of their geometric transit. While at older ages, Ly-alpha absorption is still significant (and possibly saturated for the lower mass planet), the H-alpha absorption nearly disappears. This is because the extended atmosphere of neutral hydrogen becomes predominantly in the ground state after ~1.2 Gyr.

Tables

Table 1. Table 1: Coefficients to fits given in Equations ( 19 ) and ( 20 ) for the three considered host star rotations at various ages, as shown in Figure 3 , for the 1- M jup subscript 𝑀 jup M_{\rm jup} planet.

age range [Myr]	$a$	$b$	$c$	$d$
Slow
10 – 30	13.8	1.58	1.71	0.107
30 – 250	12.2	0.532	1.66	0.0902
250 – 5000	14.4	1.37	1.79	0.175
Intermediate
10 – 30	13.1	0.714	1.78	0.128
30 – 250	12.9	0.621	1.71	0.0991
250 – 5000	15.5	1.69	1.75	0.116
Fast
10 – 30	12.2	-0.0918	1.94	0.174
30 – 250	12.7	0.287	1.70	0.0695
250 – 5000	17.4	2.24	1.67	0.0414

Table 2. Table 2: The same as in Table 1 , but for the 0.3- M jup subscript 𝑀 jup M_{\rm jup} planet.

age range [Myr]	$a$	$b$	$c$	$d$
Slow
10 – 30	14.4	1.52	1.91	0.104
30 – 250	12.7	0.342	1.82	0.0758
250 – 5000	14.3	0.977	1.95	0.160
Intermediate
10 – 30	13.7	0.630	1.91	0.106
30 – 250	13.4	0.488	1.86	0.0897
250 – 5000	15.3	1.26	1.93	0.125
Fast
10 – 30	12.9	-0.138	2.02	0.136
30 – 250	13.4	0.227	1.85	0.0662
250 – 5000	17.0	1.76	1.87	0.0716

Equations

$$u\frac{du}{dr}=-\frac{1}{\rho }\frac{dP}{dr}-\frac{G{M}_{\mathrm{pl}}}{r^2}+\frac{3G{M}_{\ast }r}{a^3},$$

$$\rho u\frac{d}{dr}\left[\frac{k_{B}T}{(\gamma -1)m}\right]=\frac{k_{B}T}{m}u\frac{d\rho }{dr}+Q-C,$$

$$\frac{n_n{F}_{\mathrm{EUV}}{e}^{-\tau }{\sigma }_{{\nu }_0}}{e_{\mathrm{in}}}=n_{+}^2{\alpha }_{\mathrm{rec}}+\frac{1}{r^2}\frac{d}{dr}(r^2{n}_{+}u),$$

$$\frac{d(r^2\rho u)}{dr}=0.$$

$$R_{\mathrm{Roche}}=\frac{a}{R_{\mathrm{pl}}}{\left(\frac{M_{\mathrm{pl}}}{3M_{\ast }}\right)}^{1/3}.$$

$$E_{\mathrm{kinetic},\mathrm{output}}=E_{\mathrm{irradiation},\mathrm{input}}$$

$$\frac{\dot{M}{u}_{\mathrm{term}}^2}{2}=ϵF_{\mathrm{EUV}}(\pi R_{\mathrm{pl}}^2).$$

$${\dot{M}}_E=ϵ\frac{F_{\mathrm{EUV}}(\pi R_{\mathrm{pl}}^2)}{G{M}_{\mathrm{pl}}/R_{\mathrm{pl}}}.$$

$$K=1-\frac{3R_{\mathrm{pl}}}{2R_{\mathrm{Roche}}}+\frac{R_{\mathrm{pl}}^3}{2R_{\mathrm{Roche}}^3}.$$

$${\dot{M}}_E^{\mathrm{corr}}=ϵ\frac{F_{\mathrm{EUV}}(\pi R_{\tau =1}^2)}{KG{M}_{\mathrm{pl}}/R_{\mathrm{pl}}}.$$

$${\tau }_{\nu }=\int n_n\sigma {\varphi }_{\nu }dz,$$

$$\sigma =\frac{\pi e^{2}f}{m_{e}c}$$

$${\varphi }_{\nu }=\frac{{\lambda }_0}{\sqrt{\pi }{u}_{\mathrm{th}}}\frac{\xi }{\pi }{\int }_{-\infty }^{\infty }\frac{e^{-w^2}}{{\chi }^2+(\Delta u/u_{\mathrm{th}}-w{)}^2}dw$$

$$\frac{I_{\nu }}{I_{\star }}=e^{-{\tau }_{\nu }}.$$

$$\Delta F_{\nu }=\frac{\int \int (1-e^{-{\tau }_{\nu }})dxdy}{\pi R_{\star }^2}.$$

$$\Delta F_{\nu }=\frac{\pi (R_{\mathrm{pl},\nu }^{\mathrm{eff}}{)}^2}{\pi R_{\ast }^2}.$$

$$t_{\mathrm{dur}}\approx \frac{P_{\mathrm{orb}}{R}_{\star }}{\pi a}=\frac{2R_{\star }}{\sqrt{G{M}_{\ast }/a}}=2.75\mathrm{h},$$

$$W_{\mathrm{H}\alpha }={\int }_{v_i}^{v_f}(\Delta F_{\nu }-\Delta F_{\mathrm{geom}})dv,$$

$$\frac{\dot{M}}{[g/s]}=10^a{\left(\frac{\mathrm{age}}{[\mathrm{Myr}]}\right)}^{-b}$$

$$\frac{u_{\mathrm{term}}}{[km/s]}=10^c{\left(\frac{\mathrm{age}}{[\mathrm{Myr}]}\right)}^{-d}$$

Full text

††pubyear: 2019††pagerange: Evolution of atmospheric escape in close-in giant planets and their associated Ly$\alpha$ and H$\alpha$ transit predictions–A

1 Introduction

Following the first detection of an exoplanetary atmosphere around the hot-Jupiter HD209458 b (Charbonneau et al., 2002), there has been a plethora of studies aimed at better understanding the atmospheres of exoplanets. Whether a planet can retain an atmosphere for a significant time period is essential for its habitability, as atmospheres regulate the planetary surface temperature, as well as absorb harmful solar or cosmic rays (e.g., Kasting et al., 1993; Lammer et al., 2009). Consequently, studies regarding the stability of planetary atmospheres are critical in the ongoing searches for habitable planets.

A planet that is too close to its host star, however, receives intense high-energy (i.e., X-rays and ultraviolet, UV) irradiation, which can cause atmospheres to heat, expand and ‘evaporate’. For many years now, the Hubble Space Telescope (HST) has been a key instrument for observing escaping planetary atmospheres in the UV (e.g. Vidal-Madjar et al., 2003). Major advancements in the field can be expected with new missions soon to be launched, such as the Colorado Ultraviolet Transit Experiment (Fleming et al., 2017), and the James Webb Space Telescope. Predicting prime candidates for future observations of escaping atmospheres is therefore important both to guide these observations as well as for interpreting their findings.

Strong atmospheric escape has been detected in several close-in exoplanets so far using transmission spectroscopy in the UV (e.g., Vidal-Madjar et al., 2003; Fossati et al., 2010; Kulow et al., 2014; Ehrenreich et al., 2015; Bourrier et al., 2016). In particular, as the composition of gas giants consists predominantly of hydrogen, their atmospheres show strong absorption in the Ly$\alpha$ line. Unfortunately, the centre of this line is contaminated in our geocorona and can be significantly absorbed by neutral hydrogen in the interstellar medium (Wood et al., 2002; Ben-Jaffel, 2008). This renders the line core unusable and the signatures of escaping atmospheres can only be observed in the wings of the Ly$\alpha$ line.

Another promising line for conducting transmission spectroscopy is H$\alpha$, as this line can be observed with ground-based spectrographs and it does not suffer contamination at the line centre, thus allowing us to probe the accelerating material at the bottom of the escaping atmosphere. However, as the number of hydrogen atoms initially in the first excited state ($n=2$) can be much lower than the population at ground level ($n=1$), transit depths in H$\alpha$ are expected to be smaller than in Ly$\alpha$ transits (Jensen et al., 2012; Cauley et al., 2015). Another complication is that an active host star shows variability in H$\alpha$, making it harder to disentangle signatures from the star and from the planetary atmosphere (Cauley et al., 2017). Recently, escaping atmospheres were detected using the metastable helium triplet at $\sim 10830$Å (Nortmann et al., 2018; Allart et al., 2019; Alonso-Floriano et al., 2019), opening up new possibilities for probing planetary atmospheres with ground-based instrumentation in the infrared, such as with CARMENES (Quirrenbach et al., 2014) and SpIRou (Cloutier et al., 2018).

The reason why atmospheres of close-in giants are more susceptible to strong evaporation relies on two factors: 1) the intense high-energy irradiation they receive from their host stars, due their close proximity to the host, and 2) on their relatively low gravitational potential, required for a planet to ‘hold on’ to their atmospheres. Young exoplanets are more susceptible to escape, because young stars have higher X-ray and UV fluxes (Ribas et al., 2005; Johnstone et al., 2015b) and because planets at young ages are more ‘puffed up’, and thus have a lower gravitational potential. As the system evolves and the high-energy flux decreases with time, escape rates are expected to decrease.

In this paper, we model the evolution of atmospheric escape on close-in giant planets. In particular, we focus on two different types of gas giants: a Jupiter-mass planet and a Saturn-mass planet (1 and 0.3 $M_{\rm jup}$, respectively). In both cases, we consider these planets to be at an orbital distance of 0.045 au or, equivalently, at a $3.5$ day orbit around a solar mass star. This distance is just inside the short-period ‘Neptunian desert’, which shows a lack of planets with masses between 0.03 and 0.3 $M_{\rm jup}$ at orbital periods $\lesssim 5$–$10$ days (Mazeh et al., 2016). The paucity of detected planets within this desert is suggestive of significant mass loss in planets with masses $\lesssim 0.3~{}M_{\rm jup}$. This will indeed be demonstrated in our results.

Our paper is presented as follows. Section 2 presents the two main inputs required in our study of the evolution of atmospheric escape: a prescription for how the high-energy flux of the host star evolves during the main-sequence phase and a prescription describing the contraction of the planet after its formation. These two factors are then considered in our atmospheric escape models, which are described in Section 3, along with a description of the ray tracing technique used to simulate the spectroscopic transits. Our results are shown in Section 4: for each age ranging from 10 to 5000 Myr, we calculate the physical properties of the escaping atmospheres and predict their Ly$\alpha$ and H$\alpha$ transit profiles, which can be compared to observations. Section 6 presents a discussion of our results and our conclusions.

2 Evolution of planetary radii and stellar EUV Flux

The evolution of atmospheric escape of a close-in planet depends on two important factors:

as the host star evolves, its activity declines due to spin down, resulting in declining fluxes in the extreme ultraviolet (EUV) and 2. 2.

as the planet evolves, cooling causes it to contract with time.

Therefore, to simulate the evolution of a close-in giant, we couple evolving predictions of EUV flux (Johnstone et al., 2015c) and planetary radius (Fortney & Nettelmann, 2010) with our hydrodynamic escape model, which will be described in Section 3. The evolution of these input parameters are shown in Figure 1 and described next.

The EUV emitted from the host star originates from its magnetically heated, chromospheric and coronal plasma (e.g., Jardine et al., 2006). Given that magnetic activity is linked to rotation, the EUV flux depends on the host star’s rotational evolution. In rotational evolution studies, it is common to separate stars in different rotational tracks: those representing stars born as slow, intermediate or fast rotators (Gallet & Bouvier, 2013; Johnstone et al., 2015a). Due to the rotation-activity relations, stars that are born as slow rotators have overall lower EUV flux throughout their lives than those that were born as fast rotators (Tu et al., 2015). As the host star evolves, its magnetic activity declines due to spin down (Vidotto et al., 2014). At $\sim$ 1 Gyr, the rotational evolution paths for solar-mass stars converge to a unique solution (Skumanich, 1972; Gallet & Bouvier, 2013). This convergence is thus also predicted to occur in the evolution of the EUV flux ($F_{\rm EUV}$).

Here, similarly to Johnstone et al. (2015c), we consider three host star rotational evolutions, resulting in three differing EUV flux evolutions for fast, intermediate and slow rotators (Figure 1, top). This allows us to study the influence of initial conditions of the star’s rotation on atmospheric escape of an orbiting close-in giant.111The values predicted by Johnstone et al. (2015c) results in $F_{\rm EUV}=$ 6.7 erg s${}^{-1}{\rm cm}^{-2}$ at 1 au assuming a solar-like star at the approximate solar age. To match the observed flux of 1 erg s${}^{-1}{\rm cm}^{-2}$ from Ribas et al. (2005), we normalise the curves in Johnstone et al. (2015c) by dividing them by 6.7. The flux of 1 erg s${}^{-1}{\rm cm}^{-2}$ is also the one adopted in the fiducial case of Murray-Clay et al. (2009), on which our hydrodynamic model is based.

A planetary evolution model is required to calculate the change in planetary radius $R_{\rm pl}$ with time. Here, we use the results from Fortney & Nettelmann (2010) to describe the evolution of a 1-$M_{\rm jup}$ planet (magenta line in the bottom panel of Figure 1) and of a 0.3-$M_{\rm jup}$ planet (orange line). These curves represent the calculations with no core in Fortney & Nettelmann (2010), for planets orbiting a solar-like star at 0.045 au, and they do not account for the mass loss. Ideally, as the atmosphere is removed from the planet, one would also need to take into account the reduction in the total mass of the planet, and consequently the corresponding change in radius. For example, this could be done by coupling atmospheric losses to planetary evolution models in a self-consistent way. This is left for a future work. Although mass loss is not self-consistent considered in the adopted radius evolution curves, it is unlikely that this would affect the 1-$M_{\rm jup}$ planet, as total mass lost represents a small fraction of the total planetary mass, as we will see further on. However, this correction is more relevant for the 0.3-$M_{\rm jup}$, as lower-mass planets have stronger evaporation. We will discuss this further in Section 4.

3 Atmospheric escape and transit calculations

3.1 Hydrodynamic escape model

Because of the strong mass-loss rates inferred in observations ($>10^{9}$g/s, Ehrenreich & Désert 2011), atmospheric escape in close-in planets can be described using a fluid approximation. For escape to be considered in this regime, the atmosphere has to be collisional, whose indicator is given by the Knudsen number, $K_{n}={\lambda_{\rm mfp}}/{h}$, which is the ratio between the mean free path of an atmospheric particle $\lambda_{\rm mfp}$ and the atmospheric scale height $h$. The motion of a planetary atmosphere for which $K_{n}\ll 1$ is best described by a fluid, in an escape process known as hydrodynamic escape. This escape process occurs when heating in the collisional region of an atmosphere causes an upward pressure gradient force, which drives a bulk, radial outflow (Catling & Kasting, 2017).

Our numerical model is based on the model presented in Murray-Clay et al. (2009). We treat the escaping atmosphere as a fluid, solving the equations of fluid dynamics in a co-rotating frame. These coupled differential equations are ionisation balance as well as the conservation of mass, energy and momentum. To achieve convergence of the wind solution, we use a shooting method based on the model of Vidotto & Jatenco-Pereira (2006).

In steady state, the momentum equation in spherical symmetry can be written as

[TABLE]

where $r$ is the radial coordinate from the centre of the planet, $\rho$ and $u$ represent the atmospheric mass density and velocity, respectively, and $P$ is the thermal pressure. The first term on the right side of Equation (1) is the thermal pressure gradient while the second term represents attraction due to gravity. This equation is analogous to the momentum equation for a stellar wind (Parker, 1958) with an additional term on the right side due to tidal effects. The tidal term is the sum of the centrifugal force and differential stellar gravity along the ray between the planet and star (García Muñoz, 2007). In this paper, we use the terms ‘hydrodynamic escape’ and ‘planetary wind’ interchangeably.

Conservation of energy requires that

[TABLE]

where $k_{B}$ is the Boltzmann constant, T the temperature, $\gamma={5}/{3}$ is the ratio of specific heats for a monatomic ideal gas and $m$ is the mean particle mass, which ranges from $0.5m_{H}$ to $1m_{H}$, for a fully ionised to fully neutral atomic hydrogen plasma, respectively. Here, $m_{H}$ is the mass of atomic hydrogen. The term on the left indicates the change in the internal energy of the fluid. The first term on the right represents cooling due to gas expansion, while the second ($Q$) and third ($C$) are heating and cooling terms. This equation is, again, typically used in stellar wind models, although with different heating and cooling terms (Vidotto & Jatenco-Pereira, 2006).

Our planetary wind model assumes that the EUV flux is concentrated at one photon energy $e_{\rm in}=20$ eV, thus the volumetric heating rate [erg/s/cm3] can be written as $Q=\epsilon F_{\rm EUV}e^{-\tau}\sigma_{\nu_{0}}n_{n}$, where $F_{\rm EUV}$ is the EUV flux received by the planet, $\tau$ the optical depth to ionising photons, $n_{n}$ the number density of neutral hydrogen and $\sigma_{\nu_{0}}$ is the cross section for the photoionisation of hydrogen given by $\sigma_{\nu_{0}}=6\times 10^{-18}\left({e_{\rm in}}/{13.6{\rm eV}}\right)^{-3}=1.89\times 10^{-18}{\rm cm}^{2}$ (Spitzer, 1978). Of the total received energy flux $F_{\rm EUV}$, a fraction $\epsilon$ is converted to thermal energy to heat the atmosphere. In our simulations, we follow the approach by Murray-Clay et al. (2009) and assume $\epsilon=1-{13.6~{}\rm{eV}}/{e_{\rm in}}=0.32$.

We assume that the escaping atmospheric gas cools by radiative losses resulting from collisional excitation, also known as Ly$\alpha$ cooling. The volumetric cooling rate is $C=7.5\times 10^{-19}n_{+}n_{n}\exp[{-1.183\times 10^{5}}/{T}]$, where $n_{+}$ is the number density of protons, which is equivalent to the number density of electrons in a purely hydrogen plasma.

The ionisation balance is given by

[TABLE]

where the rate of photoionisations (first term) balances the sum of the rates of radiative recombinations (second term) and of the advection of ions (third term). Here, $\alpha_{\rm rec}=2.7\times 10^{-13}(T/10^{4})^{0.9}$ is the case B radiative recombination coefficient for hydrogen ions (Storey & Hummer, 1995; Osterbrock & Ferland, 2006). The ionised fraction is given as $F_{\rm ion}={n_{+}}/{n_{H}}$, where $n_{H}=n_{+}+n_{n}$ is the total number density of hydrogen.

Finally, the conservation of mass requires that

[TABLE]

Our simulated planetary outflow originates from the sub-stellar point of the planet, which is the point closest to the star. We then apply our calculated solution over 4$\pi$ steradians rendering it an upper limit to atmospheric escape (Murray-Clay et al., 2009; Johnstone et al., 2015b). Therefore, from Equation (4), we have that the mass-loss rate of the escaping atmosphere is $\dot{M}=4\pi r^{2}\rho{u}$.

Similarly to stellar wind theory, the initially subsonic flow is accelerated transonically until it reaches an asymptotic speed, known as the terminal velocity $u_{\rm term}$. Convergence is reached once the values of the mass-loss rate and terminal velocity between two subsequent runs are below 1%. We then calculate the Knudsen number $K_{n}$ for each run of our simulations, so as to ensure that the atmosphere remains collisional throughout all our simulations.

Figure 2 shows typical solutions of our simulations, where we show radial profiles of various atmospheric parameters for $1~{}M_{\rm jup}$ (magenta) and $0.3~{}M_{\rm jup}$ (orange) planets orbiting a solar mass star at 0.045 au. The adopted radii for these illustrative simulations are $1.1~{}R_{\rm jup}$. For both these planets, we assume an $F_{\rm EUV}=480$ erg cm*-2* s*-1*. In addition to these parameters, which determine the physical properties of the planetary systems, our hydrodynamic models require values of the temperature $T_{0}$ and density $\rho_{0}$ at the base of the escaping atmosphere (here, assumed at 1 $R_{\rm pl}$). For all of our simulations, we assume these parameters are $T_{0}=1000$ K and $\rho_{0}=4\times 10^{-13}$ g cm*-3*, similar to values adopted in Murray-Clay et al. (2009). Note that, although $T_{0}$ and $\rho_{0}$ are free parameters of the model, Murray-Clay et al. (2009) demonstrated that large variations in these values had negligible effect on the resulting simulated escape. We confirmed their results with our models as well.

The left panels of Figure 2 show profiles of planetary wind velocity $u$, temperature $T$, and total mass density $\rho$, while the right panels show fraction of ionised hydrogen $F_{\rm ion}$, optical depth to ionising photons $\tau$ and thermal pressure $P$. In each of the panels, the cross represents the Roche lobe boundary

[TABLE]

Atmospheric particles are no longer gravitationally bound to the planet for $r>R_{\rm Roche}$, with stellar tidal gravity dominating over planetary gravity. Note that each of our calculations extend out to 10 $R_{\rm pl}$, which is above $R_{\rm Roche}$. The circle indicates the sonic point $R_{\rm sonic}$, which is the distance at which the planetary wind reaches the sound speed $c_{s}=\sqrt{{5k_{B}T}/{3m}}$. It is evident from the velocity profile that our simulated wind becomes transonic within the Roche lobe boundary. This is also the case for our simulated planets featured in the subsequent sections during all evolutionary stages.

Figure 2 shows a steep drop in density and pressure with simultaneous sharp rises in $T$ and $F_{\rm ion}$, occurring below $1.5R_{\rm pl}$. These strong variations in the atmospheric structure are the result of strong absorption of stellar EUV radiation at this distance (region with high optical depth $\tau$), which leads to the photoionisation of atmospheric hydrogen, as seen in the dramatic increase of $F_{\rm ion}$ from $\simeq 0$% to $75\%$ ($30\%$) at $1.5R_{\rm pl}$ for a 1-$M_{\rm jup}$ (0.3-$M_{\rm jup}$) planet. Ejected free electrons then heat the planetary atmosphere through collisions, causing the sharp temperature increase. The high temperatures cause the atmosphere to expand and escape hydrodynamically. At large distances, the temperature drops due to the expansion of the atmosphere (adiabatic cooling). This energy sink dominates over heating once the atmosphere becomes sufficiently ionised. Ly$\alpha$ cooling also acts against the heating and is at its strongest at higher atmospheric temperatures due to increased collisional excitation (Salz et al., 2016b).

Comparing the curves for the 1-$M_{\rm jup}$ and 0.3-$M_{\rm jup}$ planets in Figure 2, we see that the wind of the more massive planet has lower velocities due to its higher gravitational potential. Its atmosphere is also more rarefied, as seen by the density and optical depth curves. Its temperature, nevertheless, reaches higher values, leading to a more ionised atmosphere. The mass-loss rates of the 1-$M_{\rm jup}$ planet is $1.7\times 10^{9}$ g s*-1*, 26 times lower than the mass-loss rate of $4.4\times 10^{10}$ g s*-1*, for the 0.3-$M_{\rm jup}$ planet. Altogether, the lower gravitational potential of the less massive planet results in stronger planetary winds, with higher velocities and mass-loss rates.

3.2 Comparison to energy-limited escape

The energy limited approximation (Watson et al., 1981) is often adopted for calculating mass-loss rates of atmospheric escape in place of more computationally heavy hydrodynamic models. This approximation can be derived as follows. We first assume energy balance between the final kinetic energy of the planetary wind and the input energy due to irradiation

[TABLE]

For the input energy, a fraction $\epsilon$ of the stellar high-energy flux $F_{\rm EUV}$ intercepted by a planet with a cross section of $\pi R_{\rm pl}^{2}$ is assumed to drive a flow that reaches a kinetic energy $\dot{M}u_{\rm term}^{2}/2$, where $u_{\rm term}$ is the terminal velocity.

[TABLE]

where $\dot{M}$ is the mass-loss rate of the planetary wind. We further assume that the terminal velocity $u_{\rm term}$ is on the order of the surface escape velocity $v_{\rm esc}=(2GM_{\rm pl}/R_{\rm pl})^{1/2}$, we have thus that the mass-loss rate derived in the energy-limit approximation is

[TABLE]

Using an efficiency of $\epsilon=0.32$, for the 1-$M_{\rm jup}$ and 0.3-$M_{\rm jup}$ planets presented in Figure 1, we find energy limited mass-loss rates of $\dot{M}_{E}=1.9\times 10^{9}$ g s*-1* and $6.2\times 10^{9}$ g s*-1*, respectively. The energy-limit approximation overestimates the mass-loss rate computed from hydrodynamical models for the 1-$M_{\rm jup}$ planet by only a factor of 1.1. However, it underestimates by a factor of 1/7 the mass-loss rate of the 0.3-$M_{\rm jup}$ planet. Examining Equations (7) and (8), we expect that $\dot{M}_{E}>\dot{M}$, when the terminal velocities of the wind are $u_{\rm term}<v_{\rm esc}$. Indeed, for the 1-$M_{\rm jup}$ planet, terminal velocities are $\sim 20$ km/s, compared to an escape velocity of $\sim 60$ km/s. On the other hand, we expect that $\dot{M}_{E}<\dot{M}$, when the terminal velocities of the wind are $u_{\rm term}>v_{\rm esc}$. This is the case for the 0.3-$M_{\rm jup}$ planet, where the terminal velocities are above 32 km/s and escape velocity is $\sim 30$ km/s.

Recent studies showed that the energy-limited approximation can be a poor indicator of the mass lost in atmospheric escape processes (e.g. Salz et al., 2016a; Kubyshkina et al., 2018) and can underestimate the mass-loss rate of low density, highly irradiated exoplanets planets in the boil-off regime (Kubyshkina et al., 2018). To better reconcile the analytical expression with hydrodynamic models, a correction factor $K$ is often considered in the denominator of Equation (8) to account for stellar tidal forces that reduces the effect of the planetary gravitational potential (Erkaev et al., 2007)

[TABLE]

Additionally, the cross section of the planet is assumed to be $\pi R_{\tau=1}^{2}$, where $R_{\tau=1}$ refers to the distance where optical depth unity to ionising photons is reached. These two additional factors can reduce the discrepancy between the hydrodynamic calculations and the analytical limit, hence

[TABLE]

However, to calculate $R_{\tau=1}$, one needs to compute the density and optical depth profiles through hydrodynamic simulations or from a derived prescription based in such simulations. For our examples presented in Figure 2, using $R_{\tau=1}$ from our simulations, and the tidal correction factor $K$, we find $\dot{M}_{E}^{\rm corr}$ that are factors of 1.6 and 0.25 the values from our hydrodynamical models. In these cases, the ‘correction’ made the discrepancy smaller for the 0.3-$M_{\rm jup}$ planet (from 1/7 to 1/4), but increased the discrepancy for the 1-$M_{\rm jup}$, when we adopted Equation (10) instead of (8). In the comparisons done further on in this paper, we compute energy-limited mass loss from Equation (10).

3.3 Transit calculations

After solving for the hydrodynamic equations for the planetary outflow, we use a ray tracing model to simulate planetary transits as observed in both Ly$\alpha$ and H$\alpha$. This model is described in Vidotto et al. (2018), here modified for Ly$\alpha$ and H$\alpha$ transitions. In the ray tracing model presented here, we upgrade from a Doppler line profile to a Voigt line profile. As the ray tracing model from Vidotto et al. (2018) is three dimensional, we create a three-dimensional Cartesian grid with 201 evenly spaced cells in each direction. This grid is centred on the planet and each cell is ‘filled’222Similar to a solid of revolution. with the one dimensional hydrodynamic calculations of $u$, $T$, and $\rho$, taking into account that Doppler shifts result from the projection of the wind velocities along the line-of-sight. One of the axis of the grid is aligned along the observer-star line, so that the grid seen in the plane-of-the-sky is a square of $201\times 201$ cells.

To calculate a frequency-dependent transit, we use 51 velocity channels from $u_{\rm channel}=-500$ km/s to $+500$ km/s, using increments of 50 km/s for the outer velocities ($|u_{\rm channel}|=100$ to $500$ km/s), while we improve the resolution at line centre (35 channel increments within $\pm$ 100 km/s). With this, the frequency-dependent optical depth along a single ray in the direction connecting the observer to the star-planet system is

[TABLE]

where $\sigma$ is the absorption cross-section at the line centre, $\phi_{\nu}$ the Voigt line profile. The subscript ‘$\nu$’ indicates that the line profile and the optical depth are dependent on the frequency (or velocity channel) of the observation. We calculate $\sigma$ for both Ly$\alpha$ and H$\alpha$ transitions using

[TABLE]

where $f$ is the oscillator strength of the transition, $m_{e}$ the electron mass, $e$ the electron charge and $c$ the speed of light. Our values of $f$ for each of these transitions were obtained from the NIST Catalogue333 https://www.nist.gov/pml/atomic-spectra-database, which gives $f=0.416410$ for Ly$\alpha$ and $f=0.64108$ for H$\alpha$. The Voigt line profile $\phi_{\nu}$ is a convolution of a Gaussian and Lorentzian line profile accounting for both Doppler and natural broadening

[TABLE]

where $\lambda_{0}$ is the wavelength at line centre ($\lambda_{0}=1215.67$Å for Ly$\alpha$; $\lambda_{0}=6562.79$Å for H$\alpha$), $u_{\rm th}=(2k_{B}T/m_{H})^{1/2}$ and $m_{H}$ is the mass of atomic hydrogen. $\chi={A_{ji}\lambda_{0}}/({4\pi u_{\rm th}})$ is the damping parameter, where $A_{ji}$ is the transition rate ($A_{ji}=4.69860\times 10^{8}$ s*-1* for the Ly$\alpha$ transition; $A_{ji}=4.4101\times 10^{7}$ s*-1* for the H$\alpha$ one; values from NIST). The velocity offset from the line centre is $\Delta{u}={u_{\rm channel}}-{u}_{\rm LOS}$, where ${u}_{\rm LOS}$ is the line of sight flow velocity of the escaping wind. We make use of IDL’s inbuilt voigt function to calculate $\phi_{\nu}$.

In computing the transmitted spectrum, we neglect centre-to-limb variations in the stellar disc and assume that, for a given frequency, the stellar disc emits a uniform specific intensity $I_{\star}$. The fraction of ‘transmitted’ specific intensity at a given frequency during transit is hence

[TABLE]

$1-I_{\nu}/{I_{\star}}$, therefore, represents the fraction of absorbed specific intensity as a result of the absorption from both the planet disc and its atmosphere. To calculate the total absorption, we then integrated over all rays. Given our three-dimensional grid construction, the grid projected in the plane of the sky has $201\times 201$ cells. We therefore shoot 51 frequency-dependent rays through each of the $201^{2}$ grid elements. Given that the grid is larger than the projected area of the stellar disc, rays emitted from pixels outside the stellar disc are assigned a zero specific intensity.

Therefore, integrating over all these rays and dividing by the ‘flux’ of star ($\pi R_{\star}^{2}$) allows us to calculate the transit depth

[TABLE]

where $dx$ and $dy$ are the sizes of each cell. By construction, $dx=dy=2\times 10R_{\rm pl}/201\simeq 0.1R_{\rm pl}$ (the factor of 2 comes from the ‘solid of revolution’ method used and $10R_{\rm pl}$ is the extension of the hydrodynamical computations). In practice, $\Delta F_{\nu}$ represents the ratio of a frequency-dependent effective area of the planet $\pi(R^{\rm eff}_{{\rm pl},\nu})^{2}$, i.e., planet disc plus its absorbing atmosphere, by the stellar area

[TABLE]

Because of the nature of our spherically symmetric model, the effective area of the planet is that of a circle, but in 3D models, this area could take a different form (e.g. Vidotto et al., 2018). Here, we assume a transit along the centre of the stellar disc resulting in an impact parameter $b=0$. Assuming circular orbit, this gives a transit duration of

[TABLE]

where $P_{\rm orb}$ is the orbital period.

4 Evolution of Atmospheric Escape in close-in giants

To probe the evolution of the atmospheric escape, we perform escape simulations for two planet masses (0.3 and 1 $M_{\rm jup}$) orbiting stars that were born as slow, intermediate and fast rotators. This results in six different evolutionary tracks. For each track, we have more than 50 models sampling ages from 10 Myr to 5 Gyr. The ages sampled are shown as diamonds in Figure 1. In total, we ran more than 300 models to compute the evolution of atmospheric escape. For each of these models, we also predict Ly$\alpha$ and H$\alpha$ transit signals. The stellar mass ($M_{*}=1M_{\odot}$), radius ($R_{*}=1R_{\odot}$) and orbital distance ($a=0.045$ au) are fixed in our study. In Figure 3, we present the evolution of three atmospheric properties predicted using our models: the mass-loss rate (top panels), terminal velocities (middle) and the position of the sonic point (bottom). The results for our 1-$M_{\rm jup}$ and 0.3-$M_{\rm jup}$ planets are shown on the left and right panels, respectively.

The evolution of $\dot{M}$ for the three considered host star rotations are shown as solid lines in the top panels of Figure 3. We see that close-in giants experience greater levels of mass loss during younger evolutionary stages and that their mass-loss rates decrease with time. This is the result of the lowered EUV flux incident on the planetary atmosphere as well as the contracting planetary radius with evolution (Figure 1). The lowered EUV flux leads to a reduced level of heating in the atmosphere due to a lower number of photoionisations. This ultimately weakens the hydrodynamic wind as the planet evolves, resulting in the declining curves of $\dot{M}$.

The contracting planetary radius with evolution results in a strengthening planetary gravitational force on atmospheric particles. This acts to further reduce the atmospheric escape experienced with evolution by improving the planet’s ability to retain its atmospheric particles. A weaker gravitational potential is directly responsible for the stronger escape experienced by the lower mass 0.3-$M_{\rm jup}$ planet, as compared to the 1-$M_{\rm jup}$ planet. Both its lower mass and larger radii are responsible for this.

Although not considered in our models, the decreasing planetary mass resulting from atmospheric escape will soften the effect of the growing planetary gravity with age. We do not include this decrease in mass in our simulations. As a result, our mass-loss rates increasingly underestimate the true mass-loss rate as the planet ages due to a growing overestimation of planetary gravity. Fortunately, we can assume that the resulting inaccuracy is of minor importance as the effect of the neglected mass variation on the planetary gravity is negligible in comparison to that resulting from the radius variation during evolution (Section 4.1).

We used the evolutionary curves from Fortney & Nettelmann (2010), which did not consider the presence of a core. By considering a 25-$M_{\oplus}$ core, Fortney & Nettelmann (2010) predicted planets with smaller radii than those without a core. In our models, this smaller radius would lead to stronger planetary gravity and consequently to a somewhat weaker mass loss. However, this would not change the trends we present in Figure 3.

The significantly steeper $\dot{M}$ curves for the close-in giant orbiting a fast-rotating star (blue) further demonstrates the significance of stellar activity on atmospheric escape. Stellar rotation determines the EUV emission, with faster rotation producing higher emission, until a saturation plateau is reached for a fast rotating star (Figure 1). The $\dot{M}$ curves for the three rotations become indistinguishable after $\sim 1$ Gyr, owing to a convergence of the received EUV flux at this age.

The dashed lines in the upper panels of Figure 3 show the mass-loss rates estimated using the energy limit approximation (Equation 10). For the 1-$M_{\rm jup}$ planet, the energy limited approximation overestimates the mass loss when orbiting a fast rotating host star (i.e., in the case of high $F_{\rm EUV}$), while it underestimates mass-loss rates in the case of planets orbiting slowly rotating hosts. The results for the hydrodynamic case and the energy-limit case agree quite well after $\sim 500$ Myr. For the case of the 0.3-$M_{\rm jup}$ planet, the situation is quite different, in which we see that the estimates of $\dot{M}$ provided by the energy-limited approximation are very different from the hydrodynamical calculation.

The middle panels of Figure 3 show the terminal velocity ${u}_{\rm term}$ of the planetary outflow, where we see that at younger ages, larger terminal velocities are reached. We define ${u}_{\rm term}$ here as the velocity of atmospheric particles leaving the planetary atmosphere at our maximum simulated distance of 10 $R_{\rm pl}$. As with the curves of $\dot{M}$, we see a decline of the ${u}_{\rm term}$ as the planet ages. However, while mass-loss rates vary by more than 3 orders of magnitude from 10 to 5000 Myr, the terminal velocities only vary by a factor of $\sim 2$. Both $\dot{M}$ and ${u}_{\rm term}$ decline with age due to a combination of an increasingly weaker atmospheric heating and an improved ability of the planet to retain its atmospheric particles causing weaker, slower hydrodynamic winds at later evolutionary stages. The velocity of escaping neutral hydrogen atoms is of particular importance for Ly$\alpha$ observations as only absorption in line wings are observable (Section 5.1).

The bottom panels of Figure 3 show the evolution of the sonic point $R_{\rm sonic}$, which climbs to higher positions of the atmosphere with evolution. The sonic point occurs when the wind velocity reaches the speed of sound, $u=c_{s}$. With evolution, there is a slower rise in the wind velocity with distance (i.e., lower acceleration). As a result, the wind must travel a greater distance before it reaches the sound speed $c_{s}$. With age, $c_{s}$ will also be slightly lowered due to overall lower atmospheric temperatures, but this only partially counteracts the effects of the lower wind acceleration. As a result, the variations in wind velocity profiles dominate in setting the growth of $R_{\rm sonic}$ with evolution.

We draw attention to the striking resemblance of each of the curves of Figure 3 to the curves of the received $F_{\rm EUV}$ included in Figure 1, demonstrating a strong dependence of $\dot{M}$ and ${u}_{\rm term}$ on the irradiation level of the planetary atmosphere. We present analytical fits to these results in Appendix A.

4.1 Total mass lost during evolution

Figure 4 shows the total mass lost calculated over the simulated life of the 1- (top panel) and 0.3-$M_{\rm jup}$ (bottom) planets. We calculate the total mass lost by integrating mass-loss rates given in Figure 3 with respect to time. Figure 4 shows the largely differing mass losses resulting from the three considered rotations of the host star. The total mass lost is clearly dependant on stellar rotation and mass loss is greatest for planets orbiting faster rotating hosts. The mass loss experienced by each planet levels off with age due to the slowing mass-loss rates with evolution shown in Figure 3.

Past studies predicted that atmospheric escape does not play an influential role in the overall evolution of hot Jupiters (Murray-Clay et al., 2009; Owen & Wu, 2013; Ehrenreich et al., 2015). In spite of the large mass-loss rates observed, the planetary mass is still very large, constituting essentially an unending atmospheric reservoir. Indeed, over its entire evolution, we predict that the 1-$M_{\rm jup}$ planet would lose at most $1\%$ of its initial mass. We therefore expect hydrodynamic escape to have a negligible effect on mass of highly irradiated Jupiter-mass planets, owing to their large mass reservoir.

However, for the 0.3-$M_{\rm jup}$ planet, our results indicate the opposite – mass lost due to hydrodynamic escape is considerably large, and could potentially shape the planet’s evolution. In this case, the total mass lost can reach more than 20%. For orbital distances even shorter than the 0.045 au considered in this work, mass-loss rates are expected to increase even more, giving support to the idea that close-in rocky planets might have begun their evolution as Neptune-like planets, but eventually evaporated the entirety of their atmospheres (Ehrenreich et al., 2015).

Regardless of the planetary mass, orbiting a faster rotating host star leads to greater losses, with a slow rotator host producing $\sim$ 4 to 5 times less lost mass than a fast rotating host. The dashed lines in Figure 4 show the total mass lost calculated using the energy-limit approach. This approximation overestimates the total mass lost when orbiting a fast rotating host star, while it underestimates the loss for planet orbiting a slow host.

5 Spectroscopic transits in Hydrogen lines

We use the solutions of the hydrodynamic calculations presented in the previous Section to calculate spectroscopic transits in Ly$\alpha$ and H$\alpha$. To model these transits, we need to know the density of neutral hydrogen in the ground state $n=1$ (hereby $N_{1}$) and in the first excited state $n=2$ ($N_{2}$), respectively. For that, we solve the statistical equilibrium equation in the coronal-model approximation (Del Zanna & Mason, 2018). In this approximation, the spontaneous radiative decay balances the excitation process, and collisional de-excitation is neglected. Only direct excitations from the ground state $n=1$ are included. To perform this calculation, we use the CHIANTI software (v 9.0.1, Dere et al., 2019)444We use sun_photospheric_2015_scott.abund for the photospheric abundance file and ion equilibrium file chianti.ioneq. written in IDL (more specifically, we use the ch_pops procedure). In practice, the population calculation requires the temperature and electron density, which are outputs from our hydrodynamic calculations. Note that for all our simulations, the ratio $N_{2}/N_{1}$ is quite low, with values not exceeding $1$ – $2\times 10^{-5}$ for the young planets, and maximum values well below that for the older planets.555We also computed the ratio $N_{2}/N_{1}$ using the Boltzmann equilibrium equation, and found ratios larger than the ones found by solving the statistical equilibrium equation. The Boltzmann equation predicted larger H$\alpha$ transits, therefore representing an upper limit for the transit calculations. In general, though, these transit depths were not larger by more than a factor of 2 of the ones shown in the present paper. However, as hydrogen is so abundant, we still predict a reasonably large transit in H$\alpha$, as we will show in Section 5.2.

5.1 Ly$\alpha$ transits

Figure 5 shows our predicted transit signatures in the Ly-$\alpha$ line for the 1-$M_{\rm jup}$ planet. For each block of three panels in Figure 5, colour indicates the planetary age as given by the relevant colour bar. In the three upper plots (a-c), we show the line profiles of hot Jupiters orbiting a slowly rotating host; the central plots (d-f) consider an intermediate rotator, while lower plots (g-i) correspond to a fast rotating host star. The lightcurves resulting from the geometric transit of the planet itself are shown in panels a, d and g. Therefore, these plots simply reveal the transit effect of the evolving planetary radius, where the transit depths can be converted in $(R_{\rm pl}/R_{*})^{2}$.

Panels b, e and h of Figure 5 display predicted transit lightcurves as observed at the Ly$\alpha$ line centre. These lightcurves feature considerably deeper transits than those resulting from the geometric transit, as neutral hydrogen in the extended atmospheres absorbs strongly in Ly$\alpha$. In addition, the typical Ly$\alpha$ transit duration is approximately 2 hours longer than the geometric transit, with young planets producing even longer durations. These plots show the evolution of the atmospheric size: as the atmospheres of younger close-in giants are inflated further than those of older planets, they occult more of the stellar disc, thus resulting in a greater transit duration and depth in Ly$\alpha$.

As Ly$\alpha$ observations are not possible at the line centre, in panels c, f and i of Figure 5, we show the predicted Ly$\alpha$ line profile at mid-transit. In these panels, we show the occulted flux $\Delta F_{\nu}$ (Equation 15) as a function of wavelength and Doppler velocity. These profiles show saturation of Ly$\alpha$ absorption at line centre and larger absorption in the line wings for younger planets. This is a direct result of the stronger atmospheric escape experienced by these planets, resulting in higher velocities of the escaping neutral hydrogen atoms. The symmetric line profiles are a result of our hydrodynamical models being spherically symmetric. A more realistic simulation, including for example interaction with stellar wind and acceleration due to radiation pressure, would present asymmetries in blue and red wings of the line profile (e.g., Kislyakova et al., 2014; Bourrier et al., 2016; Villarreal D’Angelo et al., 2018; Vidotto et al., 2018; McCann et al., 2019; Esquivel et al., 2019).

Across the three considered host-star rotations, a faster rotating host leads to longer and deeper Ly$\alpha$ transits (line centre) for a given planetary age below $\sim 800$ Myr. These younger planets also show wider line profiles at mid-transit, when orbiting fast-rotating stars. This is because planets orbiting more active stars receive more $F_{\rm EUV}$, which then creates stronger evaporation.

In Figure 6, we present the same analysis as that shown in Figure 5, only now we consider the 0.3-$M_{\rm jup}$ planet. It is clear that each of the Ly$\alpha$ transit profiles exhibit more extreme characteristics for this less massive planet: longer duration, deeper transits and wider line profiles. These are due to both a greater inflation of the atmosphere and faster moving neutral hydrogen atoms.

5.2 H$\alpha$ transits

Similar to our predictions of Ly$\alpha$ transits, here we present our estimates of H$\alpha$ transits. An advantage of observing in H$\alpha$ is that the line centre is not contaminated as in the Ly$\alpha$ line and, perhaps, most importantly, they can be conducted with ground-based spectrographs.

Figure 7 shows the H$\alpha$ transits for the 1-$M_{\rm jup}$ and 0.3-$M_{\rm jup}$ planets, left and right panels respectively, in a similar design as Figures 5 and 6. We omit here the intermediate rotator for brevity and only show results for the slow (top panels) and fast (bottom) host-star rotations. The trends in the H$\alpha$ transits are qualitatively similar to the ones seen in the Ly$\alpha$ ones: the H$\alpha$ transits exhibit greater depths and duration than the geometric lightcurves due to absorption of H$\alpha$ photons by neutral hydrogen in the first excited state in the evaporating atmosphere. Likewise, we find that younger planets produce deeper, longer and broader H$\alpha$ transits compared to their older counterparts. For a given age below $\sim 800$ Myr ($\log({\rm age})=2.9)$, the obscuration is larger for planets orbiting fast rotators, due to the higher $F_{\rm EUV}$ of these stars. However, Ly$\alpha$ and H$\alpha$ transits are quantitatively different.

Contrary to the saturation often seen in Ly$\alpha$ for several ages and host star rotations, H$\alpha$ absorption is considerably weaker. This is due to the significantly lower density of hydrogen in the first excited state $n=2$ compared to the ground state $n=1$. As a result, the transit depths we predict in our models reach, at line centre, at most $\Delta F_{\nu}\sim$5 – 8% (a $\sim$ 3 – 4% excess from the geometric transit) at very young ages reducing, at older ages, to fractions of a percent above the geometric transit, for the 1-$M_{\rm jup}$ planet, or no excess absorption at all for the 0.3-$M_{\rm jup}$ planet. The left panels in Figure 8 show the evolution of the H$\alpha$ excess flux, i.e., above the geometric transit, for the planets simulated here.

Our results, therefore, indicate that H$\alpha$ transits could be more easily detected in planets orbiting young and/or active stars. An extended H$\alpha$ transit has indeed been observed for HD189733b, which orbits a relatively active star (Jensen et al., 2012). These authors report a planet size in H$\alpha$ that is $\sim 1.4R_{\rm pl}$, or $\Delta F_{\nu}\simeq 4.7\%$ at line centre at mid-transit. These values are actually not too dissimilar to our models of the 1-$M_{\rm jup}$ planet at younger ages. H$\alpha$ was also detected in the atmosphere of KELT-9b (Yan & Henning, 2018; Cauley et al., 2019). Although the host is not an active star, it has a strong EUV flux due to its early spectral type. As a result, KELT-9b shows a wide H$\alpha$ line profile ($\sim$ 50 km/s) and a transit depth of 1.8% (1.15% in excess of the geometric transit of 0.68%, Yan & Henning 2018). Conversely, an H$\alpha$ extended transit was not detected in the less active KELT-3 system (Cauley et al., 2017).

Our models predict that lower gravity giants, like our 0.3-$M_{\rm jup}$ case, would show a larger H$\alpha$ obscuration than the more massive planets at young ages (compare top and bottom panels in Figure 8). However, at older ages, the situation is reversed, with the 0.3-$M_{\rm jup}$ showing no excess H$\alpha$ flux above $\sim$1.2 Gyr. This is because the extended atmosphere of neutral hydrogen becomes almost entirely in the ground state at older ages. Our result agrees with the lack of absorption observed in the H$\alpha$ transits of GJ436b (Cauley et al., 2017), in spite of its impressively long and extended transit in Ly$\alpha$ (Kulow et al., 2014; Ehrenreich et al., 2015; Lavie et al., 2017). This situation can occur when temperatures are not high enough to excite hydrogen. Indeed, GJ436b has a relatively small EUV flux ($F_{\rm EUV}=10^{3}$ erg/cm2/s, Ehrenreich et al. 2015). In our simulations, this flux is reached at an age of $\sim$ 3Gyr, when we already see a lack in H$\alpha$ absorption.

We also calculate the equivalent width of our H$\alpha$ transmission spectra as

[TABLE]

where ${v_{i}}$ and ${v_{f}}$ represent the velocity range over which the integration is performed and $\Delta F_{\rm geom}$ is the geometric transit depth ($R_{\rm pl}^{2}/R_{*}^{2}$), which is wavelength-independent. Equation (18) is thus the integral of the ‘excess’ transit over a velocity (or wavelength) range. We integrate over our entire simulated velocity range of $\pm 500$ km/s or, equivalently, $\pm 10.9$Å. However, we note that decreasing the range of integration to $\pm 50$ km/s ($\pm 1$Å) gives the same equivalent width – this is because the bulk of our predicted absorption occurs within $\pm 50$ km/s.

The right panels of Figure 8 show the evolution of $W_{{\rm H}\alpha}$ for the 1- (top) and 0.3-$M_{\rm jup}$ (bottom) planets, respectively. Equivalent widths are larger at younger ages for both mass planets, reaching 20 – 40 mÅ for planets orbiting active, young stars. Our values are lower than, but still similar to, the $\sim 65$ mÅ observed for the hot Jupiter HD189733b (Jensen et al., 2018).

Both the excess transit and the equivalent width of the H$\alpha$ line (Figure 8) show trends that are similar to the $F_{\rm EUV}$ evolution with rotation of the host (Figure 1 top): the curve for fast rotator lies above the curve for slow rotator. However, contrary to the saturation of flux seen at early ages, which is also reflected in the curves of $\dot{M}$ (Figure 3), the H$\alpha$ line shows neither a plateau-like saturation in transit depth nor in the equivalent width (Figure 8). This can be seen in Figure 9, were we plot the excess transit versus the escape rates. At younger ages, from 10 to 400 Myr (between crosses and squares), we see relatively small variations in $\dot{M}$ for planets around fast rotating stars (in the saturated regime), but a significant drop in excess absorption.

6 Discussion and conclusions

In this work, we studied the evolution of atmospheric escape in close-in giants. For that, we modelled escape of a higher gravity planet (1.0 $M_{\rm jup}$) and a lower gravity one (0.3 $M_{\rm jup}$). Both planets were assumed to be at 0.045 au, orbiting stars whose EUV fluxes evolve from early ages (10 Myr) to about 5 Gyr. Given that the host star’s fluxes evolve differently depending on their initial rotation, we considered three evolutionary paths for the EUV flux, by following host stars that initiated their lives as slow, intermediate or fast rotators (Johnstone et al., 2015c). In our models, we also took into account contraction of the planet during its evolution (Fortney & Nettelmann, 2010).

To probe the entire evolution, for the two different planet masses, and three different stellar rotations, we run more than 300 simulations of atmospheric escape, using a hydrodynamic escape model. We found that younger planets experience higher levels of atmospheric escape, owing to a favourable combination of higher irradiation levels and weaker planetary gravity (due to their larger radii). Lesser mass planets experience even higher escape rates than the more massive ones. Altogether, the lower gravitational potential of the less massive planet results in stronger planetary winds, with higher velocities and mass-loss rates.

In spite of the relatively large escape rates ($10^{9}$ – $10^{12}$ g/s), we estimated that, over 5 Gyr, the 1-$M_{\rm jup}$ planet would lose at most 1% of its initial mass due to escape. However, for the 0.3-$M_{\rm jup}$ planet ($10^{10}$ – $10^{13}$ g/s), up to 20% of its mass could evaporate through hydrodynamic escape. This strong escape in the lower mass planet is in line with observations of the short-period ‘Neptunian desert’, which shows a lack of planets at $\lesssim 5$–$10$-day orbit, with masses between 0.03 and 0.3 $M_{\rm jup}$ (Mazeh et al., 2016). The desert is suggestive of significant mass loss in planets with masses $\lesssim 0.3~{}M_{\rm jup}$, which is indeed confirmed by our models.

Using the results of our hydrodynamic models, we then computed spectroscopic transits in both Ly$\alpha$ and H$\alpha$ lines. The Ly$\alpha$ line is saturated at line centre and at mid-transit (we assume no impact parameter for simplicity), while H$\alpha$ transits produce transit depths of at most 3 – 4% in excess of their geometric transit at younger ages. Planets older than 200 Myr (or 600 Myr, if orbiting around fast rotating star) would have H$\alpha$ absorption $\lesssim$1% in excess of the geometric transit. While at older ages ($\gtrsim$1.2 Gyr), Ly$\alpha$ absorption is still significant (and the line could even be saturated, as in the case of the lower mass planet), the H$\alpha$ absorption nearly disappears. This is because the extended atmosphere of neutral hydrogen becomes almost entirely in the ground state at such ages.

Our models, however, predict line profiles that are symmetric, in contrast to asymmetric profiles derived from observations (e.g., HD189733b (Lecavelier des Etangs et al. 2012); HD209458b (Vidal-Madjar et al. 2003; Ben-Jaffel 2008); GJ436b (Ehrenreich et al., 2015)). This is because our escape models are spherically symmetric. To overcome this limitation, we would need to consider additional physical ingredients in our models, such as charge-exchange between stellar wind protons and atmospheric neutrals, acceleration of neutral atoms due to radiation pressure and interaction of planetary flows with a stellar wind. Several multi-dimensional simulation studies have moved into this direction (e.g., Kislyakova et al., 2014; Khodachenko et al., 2015; Bourrier et al., 2016; Villarreal D’Angelo et al., 2018; Vidotto et al., 2018; McCann et al., 2019; Esquivel et al., 2019; Daley-Yates & Stevens, 2019; Debrecht et al., 2019; Carolan et al., 2019). These physical mechanisms are also required to reproduce the large velocities seen in Ly$\alpha$. While these model provide more physical intuition in the details of the atmospheric escape, they can be computationally expensive.

The advantage of our models, versus more complex ones, is thus the ability of computing a large number of simulations due to their low computational cost. As a result, we can make general predictions that could help guide observations. We showed here that younger systems show deeper transits as well as larger widths in both H$\alpha$ and Ly$\alpha$ lines. In theory, this indicates that atmospheric escape is easier to detect during early ages. However, in practice, at younger ages, stars are more active, generating variability that could mimic signatures of atmospheric transits (Barnes et al., 2016; Jensen et al., 2018). This effect makes escape observations at young ages and/or around active stars more difficult to interpret and, thus, models, like ours, can assist in the interpretations of transit signatures.

Acknowledgements

We acknowledge funding from the Irish Research Council Consolidator Laureate Award 2018. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), University of Cambridge (UK) and NASA Goddard Space Flight Center (USA). The authors thank C. Villarreal D’Angelo and A. Esquivel for their insightful comments and discussion.

Appendix A Analytical fits to mass-loss rates and terminal velocities

Here, we present power-law fits of the evolutionary curves of mass-loss rate and terminal velocity presented in Figure 3. These quantities take the form

[TABLE]

[TABLE]

where $a$, $b$, $c$ and $d$ are shown in Tables 1 and 2. These expressions give errors of up to 1.5% for $\dot{M}$ and 1.6% for uterm.

Equation 19 offers an alternative to the widely used energy-limited approximation for $\dot{M}$. In addition, a priori knowledge of the effective radius $R_{\tau=1}$ included in the energy limited approximation (Equation 10) is not required by our expression for mass-loss rates.