Issue 
A&A
Volume 619, November 2018



Article Number  A151  
Number of page(s)  21  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201833737  
Published online  16 November 2018 
Grid of upper atmosphere models for 1–40 M_{⊕} planets: application to CoRoT7 b and HD 219134 b,c
^{1}
Space Research Institute, Austrian Academy of Sciences,
Schmiedlstrasse 6,
8042
Graz,
Austria
email: daria.kubyshkina@oeaw.ac.at
^{2}
Institute of Computational Modelling of the Siberian Branch of the Russian Academy of Sciences,
660036 Krasnoyarsk, Russia
^{3}
Siberian Federal University,
660041,
Krasnoyarsk,
Russia
^{4}
Institute for Astronomy, University of Vienna,
Türkenschanzstrasse 17,
1180
Vienna,
Austria
^{5}
Institute of Physics/IGAM, University of Graz,
Universitätsplatz 5,
8010
Graz,
Austria
Received:
28
June
2018
Accepted:
6
September
2018
There is growing observational and theoretical evidence suggesting that atmospheric escape is a key driver of planetary evolution. Commonly, planetary evolution models employ simple analytic formulae (e.g. energy limited escape) that are often inaccurate, and more detailed physical models of atmospheric loss usually only give snapshots of an atmosphere’s structure and are difficult to use for evolutionary studies. To overcome this problem, we have upgraded and employed an existing upper atmosphere hydrodynamic code to produce a large grid of about 7000 models covering planets with masses 1–39 M_{⊕} with hydrogendominated atmospheres and orbiting latetype stars. The modelled planets have equilibrium temperatures ranging between 300 and 2000 K. For each considered stellar mass, we account for three different values of the highenergy stellar flux (i.e. low, moderate, and high activity). For each computed model, we derived the atmospheric temperature, number density, bulk velocity, Xray and EUV (XUV) volume heating rates, and abundance of the considered species as a function of distance from the planetary centre. From these quantities, we estimate the positions of the maximum dissociation and ionisation, the massloss rate, and the effective radius of the XUV absorption. We show that our results are in good agreement with previously published studies employing similar codes. We further present an interpolation routine capable to extract the modelling output parameters for any planet lying within the grid boundaries. We used the grid to identify the connection between the system parameters and the resulting atmospheric properties. We finally applied the grid and the interpolation routine to estimate atmospheric evolutionary tracks for the closein, highdensity planets CoRoT7 b and HD 219134 b,c. Assuming that the planets ever accreted primary, hydrogendominated atmospheres, we find that the three planets must have lost them within a few Myr.
Key words: hydrodynamics / planets and satellites: atmospheres / planets and satellites: dynamical evolution and stability / planets and satellites: individual: CoRot7 b / planets and satellites: individual: HD 219134 b / planets and satellites: individual: HD 219134 c
© ESO 2018
1 Introduction
The results of the NASA Kepler mission have revealed the presence of a large variety of planetary systems, with structures and geometries often very different from the solar system. The detection of a large number of extrasolar planets (hereafter exoplanets) with masses and radii in between those of Earth and Neptune is a striking example (e.g. Bonfils et al. 2013; Mullally et al. 2015).
SuperEarths and miniNeptunes, absent in the solar system, are extremely common and are easier to detect and characterise compared to Earthmass planets. Therefore, these planets are raising great interest and are among the primary targets for planetfinding and characterisation missions such as CHEOPS (Broeg et al. 2013), TESS (Ricker et al. 2015), CUTE (Fleming et al. 2018), JWST (Gardner et al. 2006; Deming et al. 2009), PLATO (Rauer et al. 2014), and ARIEL (Tinetti et al. 2017).
SuperEarths and miniNeptunes can have a large variety of average densities ranging from being consistent with rocky planets up to planets with thick hydrogendominated envelopes (e.g. Weiss & Marcy 2014; Lopez & Fortney 2014; Howe et al. 2014; Wolfgang et al. 2016; Cubillos et al. 2017a). Assuming planets were formed inside the protoplanetary disk, and thus accreted a gaseous envelope, the rocky planets most likely lost their primordial hydrogenrich envelope through escape, while the lowdensity planets still retain their primordial atmosphere. Fulton et al. (2017) revealed the presence of a dichotomy in the radius distribution of the superEarths and miniNeptunes discovered by the Kepler mission (see also Van Eylen et al. 2017; Fulton & Petigura 2018), which Owen & Wu (2017) and Jin & Mordasini (2017) interpreted as being the result of atmospheric escape processes occurring during the first few hundred million years following the dispersal of the protoplanetary disk (see Ginzburg et al. 2018, for an alternative explanation).
These works (see also e.g. Lundkvist et al. 2016) clearly showed that atmospheric escape is likely to play a major role in shaping the currently observed exoplanet population and massradius distribution. Atmospheric escape is also becoming more relevant in the characterisation of lower atmospheres: for example, Cubillos et al. (2017b) showed that the penetration depth in the planetary atmosphere of the highenergy stellar radiation (XUV: EUV + Xray) can be used to constrain the lower pressure levels for the presence of clouds.
In this study, we have focussed on 1–40 M_{⊕} planets that have accreted a primordial hydrogendominated envelope while forming inside the protoplanetary nebula (see e.g. Stökl et al. 2016). Once released from the protoplanetary nebula, planets experience a short phase of extreme hydrodynamical thermal escape, caused mostly by their high temperature and low gravity (Stökl et al. 2015; Owen & Wu 2016; Ginzburg et al. 2016; Fossati et al. 2017). This socalled boiloff phase is followed by a much longerphase in which the hydrodynamic atmospheric escape is driven by the incident stellar XUV flux (e.g. Lammer et al. 2003). Usually, this type of escape is called blowoff and the atmospheric escape rates can be estimated using the energy or recombinationlimited formulas (Watson et al. 1981; Lecavelier des Etangs et al. 2004; Erkaev et al. 2007; Lammer et al. 2009; Ehrenreich & Désert 2011; Salz et al. 2016; Chen & Rogers 2016). Both these escape conditions are different from classical Jeans escape, in which only the fraction of particles lying close to or above the exobase with velocities larger than the planetary escape velocity leave the planet.
Overall, the energylimited formula well reproduces the escape rates obtained through detailed hydrodynamic upper atmosphere modelling, particularly for closein gas giants with atmospheres in blowoff (e.g. Lammer et al. 2009; Fossati et al. 2015; Salz et al. 2016; Erkaev et al. 2016, 2017). Because of its analytical form, hence allowing for fast computations, the vast majority of planetary evolution and population synthesis models employ the energy and recombinationlimited formalisms to model atmospheric escape for a wide range of planets subject to (very) different stellar irradiation levels (e.g. Jackson et al. 2012; Batygin & Stevenson 2013; Lopez & Fortney 2013; Jin et al. 2014; Owen & Wu 2017; Jin & Mordasini 2017; Lopez 2017). However, it has also been shown that in many cases, particularly for highly irradiated lowmass planets and for planets with hydrostatic atmospheres, the energylimited formula tends to significantly over or underestimate the escape rates for low or high mass planets respectively (e.g. Erkaev et al. 2015, 2016; Lammer et al. 2016; Salz et al. 2016; Owen & Mohanty 2016; Fossati et al. 2017, 2018).
In this work, we have followed and expanded on the approach of Johnstone et al. (2015a), who computed a small grid of upper atmosphere hydrodynamic models and extracted the massloss rates by interpolating between the grid cells to model the possible evolution of the atmosphere of earlyEarth and to avoid the assumptions connected with the use of analytical formalisms. This approach enables more reliable planetary evolution computations, appropriately accounting for boiloff, blowoff, and Jeans escape, and smoothly transitioning among the different escape regimes, without significantly affecting the computational time.
We present here a large grid of upper atmosphere hydrodynamic models computed for a wide range of parameters for 1–40 M_{⊕} planets. We also present an interpolation routine we developed to extract model output parameters, such as atmospheric temperatures, velocities, densities and hydrogen species abundances, and resulting escape rates, for any planet contained within the grid boundaries. The model grid and interpolation routine can quickly produce the results of a full hydrodynamic upper atmosphere computation for planets covered by the grid, without the need to actually run a model. This enables faster, yet more accurate, interpretation and characterisation of planetary atmospheres in comparison to the results provided by, for example, the energylimited formula. This has now become particularly important to understand the massradiusperiod distribution of the large number of planets expected to be discovered in the near future by allsky surveys such as TESS and PLATO (Rauer et al. 2014; Barclay et al. 2018). Furthermore, a grid approach enables accurate planetary evolution and population synthesis computations and the thorough exploration of trends in the characteristics of planetary upper atmospheres as a function of system parameters.
This paper is organised as follows. In Sect. 2, we present the hydrodynamic model used to compute the grid and a comparison with the literature, while in Sect. 3 we describe the grid boundaries and structure. Section 4 gives an overview of the results and provides a description of the interpolation routine. In Sect. 5, we discuss the results and present an application of the grid to the case of the lowmass, closein planets CoRoT7 b and HD 219134 b,c. In Sect. 6, we summarise our conclusions.
2 Upper atmosphere modelling
2.1 Hydrodynamical model
The construction of a large grid requires a hydrodynamic model that satisfies two basic criteria: it has to reliably compute upper atmosphere profiles within a short time and it has to be able to cover a wide range of stellar, orbital, and planetary parameters. The first point is critical because the need to cover a large parameter space requires the computation of numerous models (i.e. >1000). These criteria are well matched by the onedimensional hydrodynamic upper atmosphere model described by Erkaev et al. (2016), which has been successfully tested for a very wide range of planetary systems (e.g. Lammer et al. 2013, 2016; Erkaev et al. 2013, 2014, 2015, 2016, 2017; Fossati et al. 2017; Cubillos et al. 2017a,b).
To simplify and speed up the calculation of a large number of models in the grid, we have implemented a new computational scheme, which provides an automatic selection of an initial atmospheric profile for each planet (i.e. each point in the grid). Our hydrodynamic code includes Xray heating and cooling, which are relevant for some of the planets close to our grid boundaries. The addition of Xray heating also provides us with a further important degree of freedom relevant for young, closein planets, which are subject to strong blowoff (e.g. Kubyshkina et al. 2018). We provide below a detailed description of the modelling scheme.
We set the lower boundary of the atmospheric profile at the photospheric radius (R_{pl}), at which weconsidered the planetary atmosphere to have a temperature equal to the equilibrium temperature (T_{eq}; see Fossati et al. 2017) assuming zero Bond albedo and full energy redistribution. The upper boundary was set at the Roche radius (1)
where M_{pl} and M_{*} are the planetary and stellar masses, respectively, and d_{0} is the orbital separation. The boundary conditions at the upper limit were set to be free (the radial derivatives of the computed quantities become zero). We assumed a pure hydrogen atmosphere and that at the lower boundary the atmosphere is composed exclusively of molecular hydrogen. Following Fossati et al. (2017), for each planet we computed the pressure at the lower boundary of the atmosphere assuming solar abundances.
The chemical network implemented in the code accounts for hydrogen dissociation, recombination, and ionisation. In addition, the code accounts for Lyα cooling, XUV heating, and cooling. In the literature, the height averaged heating efficiency (η), which is the fraction of absorbed stellar XUV radiation converted into thermal energy of the atmosphere, ranges between 10 and 60% (e.g. Watson et al. 1981; Yelle 2004; MurrayClay et al. 2009; CecchiPestellini et al. 2009; Owen & Jackson 2012; Shematovich et al. 2014; Salz et al. 2016). Salz et al. (2016) showed that for Earth to Jupitermass planets η varies approximately between 10 and 25%. The implementation of a selfconsistent calculation of the heating efficiency would have made the hydrodynamic code too slow to allow the computation of a large grid. For this reason, we decided to follow the considerations of Erkaev et al. (2016) adopting for all planets a constant η value of 15% at all wavelengths.
The code solves the equations for mass, momentum, and energy conservation
where ρ, v, and T are, respectively, the mass density, bulk velocity, and temperature as a functions of the radial distance from the planetary centre r. The quantity (5)
is the planetary gravitational potential accounting for the Roche lobe effect (Erkaev et al. 2007). In Eq. (5), U_{0} = GM_{pl}∕R_{pl}, δ = M_{pl}∕M_{*}, λ = d_{0}∕R_{pl} (where d_{0} is the orbital separation), and ζ = R∕R_{pl}. The term (6)
in erg cm^{−1} s^{−1}, is the thermal conductivity of the neutral gas (Watson et al. 1981). After testing the assumption of thermal conductivity of electrons and ions, we concluded that this effect leads to negligible variations to the results and so have not included it. In particular, the largest effect is on the temperature maximum for highly irradiated, highdensity planets, which decreases by as much as 7%, with a typical decrease lying around 1–2%.
The terms P and E are the atmospheric pressure and thermal energy, which are defined as (7)
Finally, Q_{XUV}, Q_{Ly α}, and are the volume heating/ cooling rates for XUV heating, Lyα cooling, and cooling, respectively.
The spectral dependence of the stellar XUV flux varies significantly from star to star. Since we aim at computing a grid of models valid for a wide range of system parameters, it is impossible to account for the full spectral dependence of the stellar XUV flux, though the code would in principle allow it. For this reason, we assumed that the whole stellar EUV flux is emitted at a single wavelength of 60 nm (MurrayClay et al. 2009). To account for Xray heating, we assumed that the stellar Xray photons are all emitted at a wavelength of 5 nm, roughly in the middle of the Xray wavelength band.
The XUV heating function Q_{XUV} is therefore composed of two terms, Q_{EUV} and Q_{X}, which describe the heating by the EUV and Xray stellar flux, respectively. These two functions are constructed in the same way, except for the absorption crosssections and absorption functions of the stellar flux inside the planetary atmosphere that are defined at 5 and 60 nm. The total heating function thus becomes Q_{XUV} = Q_{EUV} + Q_{X}. Each heating function takes the form of (9)
where m stands for either EUV or X, σ_{m} is the absorption crosssection for the specific wavelength, and ϕ_{m} is the flux absorption function (10)
In Eq. (10), J_{m}(r, θ) is a functionwith spherical coordinates describing the spatial variation of the EUV, or Xray, flux due to atmospheric absorption (Erkaev et al. 2015)and r, in this case, corresponds to the radial distance from the planetary centre.
We defined the absorption crosssection as σ = σ_{0} , where σ_{0} = 6 × 10^{−18}, E_{i} = 13.6 eV is the hydrogen ionisation energy, and E_{λ} is the photon energy in a specific wavelength range (E_{λ} = 20 eV in the EUV and 248 eV in the Xray domain). It follows that the EUV flux absorption crosssections are 2 × 10^{−18} and 1.2 × 10^{−18} cm^{−3} for atomic and molecular hydrogen, respectively (Spitzer 1978). The Xray absorption crosssection is approximately three orders of magnitude smaller than the EUV one. This implies that the stellar Xray photons penetrate deeper into the planetary atmosphere than the EUV photons, thus heating the atmosphere closer to R_{pl}. For this reason, despite that stellar Xray fluxes are significantly smaller than the EUV fluxes, Xrays can still cause significant atmospheric heating.
We implemented Lyα cooling by adding the following function to the energy conservation equation (Watson et al. 1981) (11)
To implement cooling, we followed Miller et al. (2013), and included the following heating function (12)
where C_{n} are the temperaturedependent coefficients listed in Table 5 of Miller et al. (2013). Heating rates in Eqs. (9), (11), and (12) are given in erg cm^{−3} s^{−1}.
Since we consider nonmagnetic planets, we did not include electric conduction due to ionised components. If large enough, conduction prevents the penetration of the interplanetary magnetic field inside the ionosphere, which results in the formation of a magnetopause separating the stellar wind protons from the atmospheric ions. In addition, in case of a strongly magnetised planet, the hydrodynamic flow of the escaping ionised gas can produce electric currents in the ionosphere, which would generate a resisting force against the escaping hydrodynamic flow.
The complete list of chemical reactions and the relative crosssections (ν_{H}, , α_{H}, , ν_{diss}, γ_{H}, ν_{Hcol}, , , ) considered in the model are listed in Appendix A. The continuity equations connected with the chemical reactions are (13) (14) (15) (16)
Here, the electron density is defined as (17)
while the total hydrogen number density is the sum of the number density of all species. Finally, the mass density is (18)
For computational convenience (e.g. simplification of the continuity equations), we apply the set of normalisations presented in Appendix B. The numerical solution is based on the finite differential McCormack scheme (predictorcorrectormethod; see Erkaev et al. 2016, for more details).
2.2 Comparison with previous results
To test the modelling results, we compared the massloss rates obtained for a sample of previously (observationally and/or theoretically) studied planets with those present in the literature (Table 1). Of the four planets considered in this comparison, only GJ 436 b and Kepler11 b fall within the grid boundaries. The inclusion in the comparison of the two classical hot Jupiters, HD 209458 b and HD 189733 b, is due to the fact that these are the best studied systems in terms of atmospheric escape. For our calculations, we employed the stellar XUV fluxes and masses given by Guo & BenJaffel (2016).
We find good agreement between our values and those published in the literature, in particular for HD 209458 b, GJ 436 b, and Kepler11 b. We note that for Kepler11 b Kislyakova et al. (2014) considered mostly nonthermal escape, which is significantly smaller than the XUV driven escape, while Lammer et al. (2013) adopted a completely different lower boundary condition, which led to a significant underestimation of the massloss rates (see Lammer et al. 2016, for more details). In case of HD 189733 b, our estimation lies within the interval given by Bourrier & Lecavelier des Etangs (2013), but significantly below that of Guo & BenJaffel (2016), which appears to be an outlier compared to other estimations, and it is an order of magnitude smaller than that given by the energylimited formula. The reason may be thatEq. (11) possibly overestimates the cooling for hot Jupiters, which are optically thick to Lα in the region where the cooling peaks, preventing radiation from escaping efficiently. This was addressed in detail by Menager et al. (2013) and Koskinen et al. (2013).
The works that most closely resemble our are those of MurrayClay et al. (2009), Guo & BenJaffel (2016), and Salz et al. (2016), with which we find good agreement. We remark that none of the comparison massloss rates was computed with the energylimited approximation.
Comparison between the massloss rates obtained from our hydrodynamic modelling (Col. 6), from the energylimited formula (Col. 7), and from the literature (Col. 8).
3 Model grid
We designed the grid to model superEarths and miniNeptunes orbiting mainsequence stars. The computations were made considering the following system parameters: planetary mass M_{pl}, planetary radius R_{pl}, equilibrium temperatureT_{eq}, orbital separation d_{0}, stellar mass M_{*}, and the stellar XUV flux at the planetary orbit F_{XUV} = F_{EUV} + F_{X}. As mentioned above, we have considered the planetary radius R_{pl} to be equal to the photospheric radius assuming a clear hydrogendominated atmosphere and solar abundances.
The stellar temperature and radius change along the mainsequence phase of evolution, defined as in Yi et al. (2001). Consequently, each T_{eq} value corresponds to a range of possible orbital separations defined by the possible range of changes in stellar parameters. By fixing stellar equilibrium temperature and radius, this range of orbital separations corresponds to T_{eq} variations of the order of 10–20 K. To save computation time, we adopted a single value of the orbital separation, namely that at the centre of the range, for each T_{eq} value. Therefore, d_{0} is derived from the stellar mass and equilibrium temperature. This implies that just five input parameters of the grid are independent.
We computed models for planets with masses ranging between 1 and 39 M_{⊕} (i.e. up to twice the mass of Neptune or one tenth of Jupiter’s), with a variable step size that increases logarithmically with mass for a total of 14 planetary mass values. The planetary radius ranges between 1 and 10 R_{⊕} (i.e. up to one Jupiter radius and 2.5 times Neptune’s), in regular steps of 1 R_{⊕} (i.e. total of 10 planetary radius values). The equilibrium temperature of the planets in the grid ranges between 300 and 2000 K with regular steps of 400 K (i.e. a total of five temperature values). The cooler boundary was set to cover planets orbiting in the habitable zone, while the hotter boundary was set to ensure that our assumption of the composition of the atmosphere at the lower boundary (i.e. H_{2} dominated) holds true (Koskinen et al. 2010). Our focus is on planets orbiting early M to late Ftype stars, thus we considered stellar masses between 0.4 and 1.3 M_{⊙} for a total of five different stellar masses. We plan to extend the grid to lower mass stars, which are primary targets for various planetfinding facilities, such as CARMENES (Quirrenbach et al. 2010) and TESS (Ricker et al. 2015). Table 2 lists the values of stellar mass, equilibrium temperature, planetary radius, and planetary mass considered for the computation of the grid.
We set the range of orbital separations covered by the grid on the basis of the stellar mass and planetary equilibrium temperature, thus stellar radius (R_{*}) and effective temperature(T_{eff}). R_{*} and T_{eff} were derived considering the range of radii and effective temperatures covered by a star of each considered mass along the mainsequence on the basis of stellar evolutionary tracks (Yi et al. 2001). This leads to the fact that only a limited range of orbital separations had to be considered for each given stellar mass, saving computation time. Considering all stellar masses, the orbital separation ranges between 0.002 and 1.3 AU.
For the XUV stellar fluxes, we considered three distinct values corresponding roughly to a chromospherically active star, a moderately active star, and a quiet star. To set the high XUV flux value, we considered that the Xray saturation threshold observed for mainsequence latetype stars lies at roughly L_{X} ∕L_{bol} = 10^{a}, where L_{X} is the Xray luminosity, L_{bol} is the bolometric luminosity at the zero age main sequence (Yi et al. 2001), and a ranges between −2.5 (e.g. Reiners et al. 2014) and −3.1 (e.g. Wright et al. 2011). We therefore set the maximum Xray luminosity as L_{Xmax}∕L_{bol} = 5 × 10^{−3} and the minimum Xray luminosity as L_{Xmin}∕L_{bol} = 10^{−7}. The EUV stellar luminosity was then derived from the Xray luminosity following (SanzForcada et al. 2011) (19)
The specific Xray and EUV luminosities adopted for each stellar mass are listed in Table 3.
To avoid spending time calculating planets that probably do not exist in nature, we restricted the computations to planets with an average density larger than 0.03 g cm^{−3} (equal to the lowest known measured density; Masuda 2014) and a restricted Jeans escape parameter Λ smaller than 80 (where atmospheres are presumably stable), where (Jeans 1925; Chamberlain 1963; Öpik 1963; Fossati et al. 2017) (20)
Λ is the value of the Jeans escape parameter calculated at the observed planetary radius and mass for the planet’s equilibrium temperatureand considering atomic hydrogen, independent of the atmospheric temperature profile. We further excluded planets for which the Roche lobe is closer than 0.5 planetary radii from the surface. This cut is most relevant for the hottest planets (>1500 K) orbiting stars less massive than about 0.8 M_{⊙}. As an example, Fig. 1 shows the positions of the modelled planets in the massradius diagram at two different equilibrium temperatures.
To summarise, our grid consists of five data points for stellar mass and planetary equilibrium temperature, each, ten data points for planetary radius, 14 data points for planetary mass, and three data points for stellar XUV luminosity. This leads to a total of 10 500 models. However, because of the restrictions described above, the total number of models in the grid is reduced to 6700.
Stellar masses, equilibrium temperatures, planetary radii, and planetary masses considered for the computation of the grid.
Xray and EUV luminosities adopted for each stellar mass.
Fig. 1 Position of some of the modelled planets (black crosses) in the massradius diagram. All planets orbit a 1 M_{⊙} star and have an equilibrium temperature of 300 K (top panel) and 2000 K (bottom panel). The blue and magenta solid lines indicate the boundaries of the grid set by the cut on Λ and on bulk density, respectively. Since Λ depends on T_{eq}, the top boundary is different in the different panels. The green solid line indicates the boundary of the grid set by the cut on the Roche lobe. The position of this boundary depends on the orbital separation, thus on T_{eq}. For reference, the red dashed line indicates Earth’s density. 
4 Results
For each modelled planet, we computed the main atmospheric parameters as a function of the radial distance from the planetary centre. These include the atmospheric temperature, number density, bulk velocity, Xray and EUV volume heating rates, and abundance of the considered species (H_{2}, H, , H^{+}, , e). From these quantities, we estimated the positions of the maximum dissociation and ionisation (the distances corresponding to the maximum of the number densities of atomic and ionised hydrogen, respectively), the massloss rate Ṁ, and the effective radius of the XUV absorption that is defined as (21)
where is the XUV flux as it travels through the planetary atmosphere along the starplanet direction and is mostly determined by the density n(r). The massloss rate is computed as the product of the atmospheric density and velocity at the upper boundary. To account for the fact that we employ a onedimensional model, this value is then multiplied by the surface area of a sphere with radius equal to R_{roche}. For most planets, hydrogen dissociation occurs in a relatively narrow range of distances, which is typically smaller than one planetary radius.
Although the atmospheric parameters vary significantly across the grid, there are some common characteristics. One of the most important of these common characteristics is that the atmospheric behaviour strongly depends on Λ. For planets with low Λ values (i.e. ten or less), the atmosphere is weakly bound to the planet and experiences strong boiloff. The energy budget of these planets is determined by adiabatic cooling, and the massloss rates are not significantly affected by variations in the stellar XUV flux. With increasing Λ, the role of planetary gravity in the atmospheric dynamics decreases and the atmosphere gradually switches to being controlled by the stellar XUV heating. We find that the border between these two regimes lies at Λ values between 10 and 30, depending on the system parameters, in agreement with Fossati et al. (2017).
As an example, Fig. 2 compares the atmospheric density, temperature, velocity, and atomic hydrogen abundance profiles for two planets with Λ equal to 4.8 and 66.7. The two planets orbit a 1 M_{⊙} star, have an equilibrium temperatureof 1100 K (i.e. d_{0} = 0.075 AU), a radius of 3 R_{⊕}, and are subject to an incident XUV flux of 92.6 erg cm^{−2} s^{−1}. The planet with the lower Λ has a mass of 2.1 M_{⊕}, while that with the higher Λ has a mass of 29.1 M_{⊕}. For the lessmassive planet, we derive a value of the effective XUV absorption radius (R_{eff}) of 5.5 R_{pl}, a Roche radius of 6.5 R_{pl}, and a massloss rate of 1.1 × 10^{14} g s^{−1}. For the moremassive planet, we derive a R_{eff} value of 1.2 R_{pl}, a Roche radius of 17.1 R_{pl}, and a massloss rate of 4.0 × 10^{7} g s^{−1}. For the less massive planet, we also find that the velocity of the atmospheric particles becomes supersonic close to the Roche radius (at 6.2 R_{pl}), while for the more massive planet the particles become supersonic well below the upper boundary (at 9.1 R_{pl}).
Figure 2 shows how the density decreases with increasing distance from the planetary surface; the decrease is steeper for the more massive planet, because it hosts a more compact atmosphere (because of stronger gravity). The temperature profiles show that the stellar XUV flux efficiently heats the more massive planet, so inducing a temperature peak at the thermospheric level which reaches its maximum close to R_{eff}, where the model also indicates the presence of strong H_{2} dissociation. At higher altitudes, the atmosphere is composed fully of atomic hydrogen and is dominated by adiabatic cooling, which is caused by the atmosphere’s expansion and dominates over XUV heating. In contrast, for the less massive planet, the XUV stellar flux does not penetrate deep enough into the planetary atmosphere to cause thermospheric heating, thus the atmosphere expands adiabatically, driven by its internal heat and by the low planetary gravity. In general, the profiles of planets with low Λ values do not develop steep gradients (see e.g. Kubyshkina et al. 2018), making the definition of R_{eff} and of the position of the maximum dissociation and ionisation ambiguous. We discuss this point in more detail in Sect. 5.1.2.
Because of its relevance, for example, in understanding planetary evolution, the massloss rate is one of the key output parameters of the modelling. Figure 3 shows a few examples of how the massloss rates depend on planetary mass and radius for different T_{eq} and F_{XUV} values. All planets shown in Fig. 3 orbit a 1.0 M_{⊙} star. Appendix C presents similar plots, but for planets orbiting stars that are more or less massive than 1.0 M_{⊙}. As expected, the highest massloss rates are found for the lowest gravity planets, whose atmospheres are in boiloff. With increasing Λ (thus gravity), the massloss rates decrease first steeply and then more gradually. The dependence of the massloss rates on the stellar XUV flux tends to strengthen with increasing planetary mass.
We further checked a posteriori the validity of the hydrodynamic equations to the modelled planets, namely whether the atmosphere remains collisional (i.e. with efficient energy redistribution) up to the sonic point. This condition is satisfied if the Knudsen number Kn = λ∕l < 1, where λ is the mean free path of the gas and l = is the characteristic length scale. In Fig. 3 and in the figures in Appendix C, we show lines corresponding to Kn = 1 and Kn = 10. This indicates that the results we obtained for the highestdensity planets in our grid should be taken with caution. These are nevertheless planets for which the high bulk density disfavours the presence of a hydrogendominated atmosphere.
In Fig. 4, we compare the massloss rates as a function of Λ (top) and planetary mass (bottom) obtained from the hydrodynamic modelling with those derived from the energylimited formula (22)
where the factor K accounts forRochelobe effects (Erkaev et al. 2007). By design, the energylimited approximation works best for planets for which the atmosphere is hydrodynamic and the escape is driven by absorption of the stellar XUV flux, in other words, in blowoff. This impliesthat Eq. (25) overestimates the massloss rates for planets with hydrostatic atmospheres (see e.g. Fossati et al. 2018). The top panel of Fig. 4 shows that, being in boiloff, the massloss rates for the lowergravity planets are much higher than those predicted by Eq. (25). We also find that the Ṁ /Ṁ_{en} ratio decreases steeply with increasing Λ, having all other parameters constant. The value of Λ at which the massloss rate computed with the hydrodynamic code becomes comparable to Ṁ_{en} is about 20, in agreement with Owen & Wu (2016) and Fossati et al. (2017). Figure 4 also shows that Eq. (25) overestimates the massloss rates for planets with large Λ values.
However, the plot in the top panel of Fig. 4 might appear to be counterintuitive. This is because at large Λ values the hotter planets, thus more likely to have a stable atmosphere, present massloss rates differing more from the energylimited approximation than the cooler ones. This can be explained by the fact that for a given value of Λ hotter planets have higher masses in comparison to cooler planets (see Eq. (23)), therefore in this plot the hotter planets are more likely to have a hydrostatic atmosphere. This is clarified by the bottom panel of Fig. 4, which shows that for the higher mass planets the difference between the massloss rates computed by the hydrodynamic model and with Eq. (25) is independent of T_{eq}.
Equation (25) assumes that the entire stellar XUV energy input goes into driving the escape, but in reality part of this energy goes into running the chemical reactions, mainly ionisation of atomic hydrogen and dissociation of molecular hydrogen. Also, Erkaev et al. (2015) showed that the energylimited formula neglects kinetic and thermal energy terms in the denominator, which can also be important for some star/planet parameters; that is, even without including detailed chemistry or ionisation the results of Eq. (25) should be higher than the hydrodynamic massloss rates for XUVdriven outflows. This explains why for most planets the energylimited approximation overestimates the massloss rates. However, this is not the case when there is a significant component of thermal escape. In this case mass loss is driven partially by the stellar XUV flux and partially by the intrinsic planetary thermal energy. Such massloss rates significantly exceed those predicted by the energylimited formula.
Grid interpolation. We developed a routine which interpolates the model results over the grid parameter space considering planetary mass, planetary radius, planetary equilibrium temperature, stellar mass, and stellar XUV flux. For any system that has parameters covered by the grid, the routine extracts the density and outflow velocity, the massloss rate, the value of maximum temperature, the effective radius of XUV absorption, and the position of maximum dissociation and ionisation.
The routine performs the interpolation in the following consecutive steps.
 1.
For planets with input parameters [,,, , ] the routine finds in the grid the two closest values of stellar mass and equilibrium temperature [M_{*1}, M_{*2}], [T_{eq1}, T_{eq2}].
 2.
For each of the four combinations of [,], with i, j = 1,2, the routine finds the two closest values of the stellar XUV flux at the planetary orbital separation , with k = 1,2 (i.e. eight F_{XUV} values).
 3.
For each set of [, , ] the atmospheric parameters depend therefore only on planetary radius and mass. Planets of the same mass have different atmospheric properties for different equilibrium temperatures, but we find similar atmospheric properties for planets with similar Λ values. Therefore, we substitute planetary mass with Λ. At this point, the routine interpolates the output parameters simultaneously over the pair [R_{pl}, Λ] for each of the eight sets of [, , ]. However, for planets beyond 0.1 AU, the simultaneous interpolation on R_{pl} and Λ is not necessary, thus we reduced it to an interpolation of Λ only.
 4.
The routine interpolates the output parameters over .
 5.
The same equilibrium temperature for different stellar masses occurs at rather different orbital separations (e.g. within the grid, a T_{eq} of 300 K corresponds to distances ranging from 5×10^{−3} to 1.52 AU). Our analysis (see Sect. 5 and also Kubyshkina et al. 2018) indicates that between T_{eq} and d_{0}, the latter has the larger influence on the results, thus for the interpolation we substitute T_{eq} with d_{0}. Therefore, the routine interpolates the output parameters over for the pair of values.
 6.
The routine interpolates the output parameters over .
We developed this routine keeping in mind that the size of the grid will increase in the future, and so an interpolation routine must be capable of quickly handling the addition of grid points. Because of this we avoided the use of complicated, multidimensional interpolation functions that would require recomputing every time a new model is added to the grid.
Because the output parameters behave differently as a function of the input parameters, for almost each interpolation step and almost each output parameter, we employed a different function. For the massloss rates and the density at the Roche radius the routine interpolates over R_{pl}, Λ, and F_{XUV} according to (23) (24)
where X is either the massloss rate or the density at the Roche radius and the coefficients depend on the other system parameters. For planets in boiloff, Eqs. (26) and (27) are not accurate enough, so we used a piecewise polynomial interpolation with input and output parameters in logarithmic scale. The interpolation of the massloss rates and the density at the Roche radius over the other input parameters (i.e. d_{0} and M_{*}) is done on the basis of a linear function.
For the outflow velocity, we performed the interpolation on the [R_{pl}, Λ] pair using a third order polynomial function, which becomes linear for Λ ≳ 20. For the interpolation of the outflow velocity over the other input parameters we use a linear function.
For the interpolation of the maximum temperature over each input parameter, we employ a linear function. The routine interpolates the position of the maximum dissociation and ionisation and of the R_{eff} value on the [R_{pl}, Λ] pair using a linear function for Λ ≳ 20, while for smaller Λ values we interpolate over Λ using a function ofthe form a∕(b + Λ), where the coefficients a and b depend on the otherinput parameters. For the interpolation of the maximum temperature, maximum dissociation and ionisation, and R_{eff} value over the other input parameters (i.e. F_{XUV}, d_{0}, and M_{*}) we employ a linear function.
We performed two tests to validate the interpolation routine. We first compared the results obtained from the models with those derived using the interpolation for 500 systems randomly distributed across the grid. We found an agreement of better than 5% in 95% of the cases, while for the remaining systems the agreement was better than 20% (Fig. 5).
The second test is dedicated to check the validity of substituting T_{eq} with d_{0} for the interpolation. We used real planets for this test, namely Kepler11 b, GJ 436 b, HATP 26 b, HD 97658 b, GJ 3470 b, HATP 11 b, and 55 Cnc e, which lie within our grid boundaries, and compared the massloss rates obtained with direct modelling and interpolation. Figure 6 shows the results of this comparison, indicating that we obtain an excellent agreement for all of them. This validates our choice of interpolating on the orbital separation rather than on the planetary equilibrium temperature. We run the same test, but this time interpolating on the planetary equilibrium temperature, instead of orbital separation, obtaining significantly larger discrepancies.
Fig. 2 From top to bottom: atmospheric profiles for density, temperature, velocity, and fraction of atomic hydrogen for planets with Λ equal to 4.8 (red) and 66.7 (blue). Both planets orbit a 1 M_{⊙} star, have T_{eq} = 1100 K and R_{pl} = 3 R_{⊕}, and are subject to an incident XUV flux of 92.6 erg cm^{−2} s^{−1}. The planet with the lower Λ has a mass of 2.1 M_{⊕}, while the other one has a mass of 29.1 M_{⊕}. The density (top panel) is normalised to its value at R_{pl}. The blue and red dashed lines show the effective radii of the XUV absorption. To allow comparing planets with significantly different Roche radii, the xaxis is in units of the planetary Roche radii, starting from the planetary surface. 
Fig. 3 Logarithm of the massloss rates (colour coded) as a function of planetary mass and radius. The adopted T_{eq} and F_{XUV} values are given on the top of each panel. The equilibrium temperature increases from top to bottom, while F_{XUV} increases from left to right. All planets orbit a 1.0 M_{⊙} star. For reference, the dashed lines mark constant Λ values of 8, 20, and 50 (from top to bottom). The red lines indicate planets for which the Knudsen number at the upper boundary is equal to one (solid line) and ten (dashed line). 
Fig. 4 Ratio of the massloss rates computed with the hydrodynamic model to the energylimited formula as a function of Λ (top panel) and M_{pl} (bottom panel). The green lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 300 K, F_{XUV} = 159.4 erg cm^{−2} s^{−1}. The blue lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 700 K, F_{XUV} = 4900 erg cm^{−2} s^{−1}. The red lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 1100 K, F_{XUV} = 30784 erg cm^{−2} s^{−1}. The horizontal dotted line indicates the equality between the two values. 
Fig. 5 Relative deviation of the interpolated massloss rates (Ṁ_{int}) from the computed ones (Ṁ) as a functionof Λ. The horizontal red lines indicate deviations of 1%, 5%, and 10%. 
5 Discussion
We discuss here in more detail how the results of the grid depend on the input parameters (Sect. 5.1) and briefly explore one of its possible future applications (Sect. 5.2). The reader not interested in the technicalities of the results can skip to Sect. 5.2, which can be understood independently of the content of Sect. 5.1.
5.1 Behaviour of the atmospheric parameters as a function of input parameters
The grid allows for the detailed description of how the atmospheric structure changes with respect to the input system parameters. The behaviour of the main output parameters can be separated into common patterns. We find a common behaviour between i) massloss rates and densities of the atmospheric species (Sect. 5.1.1); ii) effective radius of the stellar XUV absorption and position of the maximum dissociation and ionisation (Sect. 5.1.2); iii) outflow velocity and atmospherictemperature (Sect. 5.1.3).
5.1.1 Massloss rates and densities of the atmospheric species
One of the major parameters controlling the longterm evolution of a planetary atmosphere is the massloss rate Ṁ, which strongly depends on the planetary gravity and orbital separation, thus the equilibrium temperature. Figure 6 shows the dependence of Ṁ on Λ. Within the parameters covered by our grid, the massloss rate varies by several orders of magnitude, with planetary gravity and Roche lobe radius being among the main parameters controlling it.
For the planets with Λ values smaller than about 20, the escape rates reach extreme values of up to 10^{20} g s^{−1}, due to a combination of low planetary gravities and high equilibrium temperatures (i.e. boiloff). The atmospheres of these planets are therefore characterised by strong thermal escape and inefficient XUV heating. The escape rates for the majority of these planets lie above the predictions of the energylimited formula, and the energy budget of the atmosphere is dominated by adiabatic cooling. In first approximation, for a given stellar mass, T_{eq}, and stellar XUV flux, the dependence of the massloss rate on Λ can be described by Eq. (27). Such high escape rates would imply a rapid escape of the atmosphere (we provide a practical example of this in Sect. 5.2).
With increasing Λ, the efficiency of XUV heating increases with XUV penetration depth (see Fig. 2) and the massloss rates become strongly dependent on the stellar XUV flux. This further dependence of the massloss rates is partly responsible for the increased spread in massloss rates at large Λ values. We also note that the spread increases with decreasing stellar mass, due to the decreasing orbital separation for the same temperature. The dependence of the massloss rates on the stellar XUV flux can be roughly described by a linear function, in agreement with the energylimited formula.
Figure 6 shows also that the Roche radius plays an important role almost exclusively when it lies below about 15 R_{pl} (in Fig. 6, the dark blue colour corresponds to Roche radii ranging from 15 up to 400 R_{pl}). It is important to recall that the Roche radius is tightly related with the orbital separation and the smallest Roche radii can be reached just for the shortest starplanet distances. As an example, within our grid, the smallest Roche radii (<3 R_{pl}) are reached only for planets lying less than 0.06 AU from the host star, while Roche radii of 15 R_{pl} are found for planets orbiting up to 0.3 AU from the host star. The massloss rates increase with decreasing Roche radii because a smaller Roche radius moves the sonic point closer to the planet, that is, to regions of higher density, which leads to an increase in massloss rate. In addition, the Roche radius decreases with decreasing orbital separation (see Eq. (1)), thus increasing T_{eq} and XUV irradiation.
To better illustrate how the escape rates change with input parameters, we set eight test planets (called Pa1, Pa2, Pb1, Pb2, Pc1, Pc2, Pd1, Pd2), whose parameters are listed in Table 4. The numbers “1” and “2” separate planets by mass, where the planets identified by the number “1” have the lower mass. The difference between the “a” and “b” and “c” and “d” planets is the stellar mass, which is higher for the former, while the difference between the “a” and “c” and “b” and “d” planets is the stellar XUV flux, which is higher for the latter. To have the planetary mass alone controlling the value of Λ, the eight planets have the same equilibrium temperature (700 K) and planetary radius (3 R_{⊕}). Figure 6 indicates the position of the eight planets in the Ṁ vs. Λ plane.
As expected, an increase in the XUV stellar flux (i.e. Pa →Pb or Pc →Pd) or decrease in Λ (i.e. Px2 →Px1, where x is any of a, b, c, or d) leads to an increase in the massloss rates. Figure 6 shows also that a decrease in stellar mass (i.e. Pa →Pc or Pb →Pd) leads as well to an increase in massloss rates. This is because, to maintain the same T_{eq}, planets orbiting around the lower mass star lie at a closer distance, thus have a smaller Roche radius and for the reasons described above have a higher massloss rate.
For the Pb planets, Fig. 6 also compares the massloss rates with those predicted by the energylimited formula assuming two different values for the effective radius and a heating efficiency of 15%, as for the hydrodynamic calculations. Since with decreasing Λ the effective radius moves farther away from the planet, the distance between the two lines in Fig. 6 increases with decreasing Λ. The dashedline in Fig. 6 presents a clear bend at Λ ≈ 4, which is caused by the fact that at small Λ values the effective radius reaches the Roche radius. Figure 6 indicates that at low Λ values the energylimited approximation significantly underestimates the escape rates (between two and three orders of magnitude; comparison to the black circles), while at large Λ values the approximation overestimates the escape rates by about one order of magnitude (see also Fig. 4).
The atmospheric densities (at the Roche radius) behave similarly to the massloss rates. The only small difference is found for planets with large Λ values (densities decrease steeper). This is because the escape rates are calculated from the product of the atmospheric density and outflow velocity at the Roche radius, where the velocity increases with increasing Roche radius, which increases with increasing planetary mass.
The dependence of the escape rates and atmospheric densities on planetary mass is essentially the same as that on Λ, though slightly less pronounced. The dependence of these parameters on T_{eq} is similar to what is displayed by the colour code in Fig. 6 and it can be described by using a loglinear approximation of the form logX = c_{1} + c_{2} T_{eq}, where X is either the massloss rate or the atmospheric density, and c_{1} and c_{2} are coefficients, which depend on the system parameters. This follows the high Λ limit of the Parker wind problem.
Fig. 6 Atmospheric massloss rate Ṁ as a functionof Λ for all computed planets. The colour code indicates the planetary Roche radius in units of R_{pl}. The position of the test planets listed in Table 4 is shown by black squares (Pa1 and Pa2), black circles (Pb1 and Pb2), purple squares (Pc1 and Pc2), and purple circles (Pd1 and Pd2). The lines indicate the predictions obtained by using the energylimited formula for the Pb test planets varying planetary mass only and assuming the value of R_{eff} equal to the planetary radius (dashdotted line) or equal to the value derived from the grid (dashed line). The black crosses and plus signs indicate the escape rates estimated for some of the known transiting exoplanets using the interpolation routine and direct hydrodynamic calculations, respectively. 
Test planets considered for the discussion of the results.
5.1.2 Effective radius and position of the maximum dissociation and ionisation
We discuss here the behaviour of three closely related parameters: the effective radius (Eq. (24)) and the position of the maximum dissociation and ionisation. As defined in Sect. 4, the position of the maximum dissociation and ionisation corresponds to the position of the maximum of n_{H} and , respectively. Figure 7 shows the position of these three quantities as a function of Λ. For most planets, they lie close to each other, except for planets with small Λ values (i.e. ≲20), where the effective radius significantly exceeds the other two. At small Λ values the effective radius can be up to ten times larger than the position of the maximum dissociation or ionisation, which stays close to one another for the whole interval of parameters. For Λ values above ~50, the difference between the position of maximum dissociation or ionisation and the effective radii lies roughly within 30%. We also found that for these planets the difference between the three values decreases slightly with increasing stellar XUV flux and planetary mass, where the latter dependence is caused by the gradual compression of the atmosphere with increasing planetary mass.
The possible range of values found in the grid for the effective radius and the position of the maximum dissociation and ionisation increases significantly with decreasing Λ and it reaches the maximum close to Λ = 5. Here the effective radius reaches the Roche radius, which decreases with decreasing Λ. The position of maximum dissociation or ionisation reaches this artificial border at smaller Λ values, in other words, approximately equal to two.
We alsofind a clear dependence of the three quantities on orbital separation. To highlight this, Fig. 8 shows the effective radius as a function of the Roche radius and of the orbital separation. Again, for each given Roche radius, the effective radius presents an upper limit, which is clearly given by R_{roche} itself. The spread in effective radii also increases with increasing R_{roche} and d_{0}. This can be understood as follows. At short orbital separations, R_{roche} is generally small and therefore there is only a small range of possible effective radii. At large orbital separations, instead, depending on the planetary and stellar masses, there is a much wider range of possible Roche radii within which the position of maximum dissociation can lie. However, Fig. 8 shows that in our grid there are very few planets with large Roche and effective radii, namely, those with a rather lowdensity and orbiting far from the host stars. For these planets, R_{roche} is large and the stellar XUV flux is too weak for the effective radius to be close to R_{pl}.
Finally, we make few remarks regarding these parameters. Unlike the definition of the effective radius (given by Eq. (24)), the definitions of the position of maximum dissociation and ionisation are not univocal. In addition, in some cases, the dissociation and ionisation profiles are very smooth, which makes the position of maximum dissociation and ionisation somewhat dependent on their definitions. It is artificial to consider that the effective radius and/or the position of the maximum ionisation and dissociation are located at the Roche lobe if their position moves beyond it. It is therefore important to keep in mind that for some particular planets these values are indicative, rather than sharp results, and that comparisons with the literature and/or with future studies should consider these issues.
Fig. 7 Effective radius (black), position of the maximum dissociation (blue), and position of the maximum ionisation (red) as a functionof Λ, for all planets in the grid. 
Fig. 8 Effective radius R_{eff} as a function of R_{roche}. Each point is colourcoded with respect to the orbital semimajor axis. The black and purple circles and squares indicate the position of the eight test planets following the same symbols as in Fig. 6. The circle and square purple symbols largely overlap. 
5.1.3 Outflow velocity at the Roche radius and maximum atmospheric temperature
We define the outflow velocity V_{roche} as the velocity of the atmospheric particles crossing the Roche lobe. Figure 9 shows the velocity at the Roche radius as a function of Λ and orbital semimajor axis. For planets with Λ values of ten or more, the velocity of the escaping material grows linearly with increasing Λ. This is caused by the gradual increase of the Roche radius with Λ, which allows for longer acceleration distances and for a smaller planetary gravitational pull at the Roche radius. This lies partly at the origin of the connection between the velocity at the Roche radius and orbital semimajor axis.
For planets with very small Λ values, the velocity behaves exactly in the opposite way, it increases with decreasing Λ as a result of the low gravity and small Roche radius. We found that the planets for which V_{roche} is larger than 2 km s^{−1} have a Roche radius located closer than 5R_{pl}, while for more extreme planets with V_{roche} greater than 4 km s^{−1}, the Roche radius is always smaller than 1.5 R_{pl}.
The dependence of the outflow velocity on the XUV flux is similar to that of the massloss rates. To illustrate this, we added the position of the eight test planets to Fig. 9. The influence of the equilibrium temperatureon the velocity at the Roche radius is however unclear, possibly because variations in the equilibrium temperatureimply simultaneous changes in the atmospheric structure and in the orbital separation, thus in the Roche radius. In Kubyshkina et al. (2018), we showed that, keeping planetary mass and radius fixed, the effects of T_{eq} variations on the velocity are negligible for planets with Λ greater than 20, while for planets with lower Λ the velocity may significantly decrease with increasing temperature.
We found that the maximum value of the atmospheric temperature T_{max} behaves similarly to the outflow velocity, except for planets with small Λ values (≲15–20). For these planets, the atmospheric temperature profile is characterised by strong adiabatic cooling, which implies thatthe maximum temperature is equal to the equilibrium temperature. For planets with larger Λ values, the maximum atmospheric temperature increases almost linearly with Λ, similarly to the velocity, but it further presents a more pronounced dependence on the stellar XUV flux.
Fig. 9 Outflow velocity at the Roche radius V_{roche} as a functionof Λ. Each point is colour coded with the value of the orbital semimajor axis. The black and purple circles and squares indicate the position of the eight test planets following the same symbols as in Fig. 6. 
5.2 Atmospheric evolution of the highdensity exoplanets CoRoT7 b and HD 219134 b,c
We present here a direct application of the grid, namely, a simple scheme allowing to infer the evolution of a planetary atmosphere subject to massloss, with the massloss rates extracted from the grid. Thanks to the dense grid and possibility to interpolate across it, one can quickly derive highresolution evolutionary tracks of planets with parameters contained in the grid. The advantage isthat the tracks are obtained making use of massloss rates computed with an hydrodynamic code, rather than more approximate methods, such as the energylimited formula. As an example, we have studied the past evolution of the possible primary atmospheres of the closein highdensity planets CoRoT7 b (Léger et al. 2009; Valencia et al. 2010; Leitzinger et al. 2011; Mura et al. 2011; Barros et al. 2014) and HD 219134 b,c (Motalebi et al. 2015; Vogt et al. 2015; Gillon et al. 2017).
5.2.1 Planetary evolution modelling scheme
We first assumed that the orbital separation and stellar mass do not change with time. It follows that the massloss rates at every moment in time depend only on the planetary radius and the amount of bolometric and XUV stellar irradiation.
To infer planetary atmospheric mass fractions (f_{at}), we employed the model described by Johnstone et al. (2015a), which relates f_{at} with planetary mass and radius. However, the approximate relation between these three quantities given by Johnstone et al. (2015a) was obtained considering atmospheres having a density of 5 × 10^{12} cm^{−3} and a temperature of 250 K at the base of the simulation. These conditions significantly differ from those of our planets. Therefore, we used the code employed by Johnstone et al. (2015a) to directly compute a grid of f_{at} values (hereafter called f_{at}grid) for planets with mass and radius in the range of interest of this work and interpolate among the grid points. For each planet considered in the f_{at}grid, the core radius R_{core} was derived assuming an Earthlike density.
One of the key parameters to set to simulate the atmospheric evolution is the initial planetary radius. Various accretion models provide an estimate of the initial planetary radius, thus atmospheric mass accreted by the planet while embedded in the protoplanetary nebula (e.g. Stökl et al. 2016), but the results are rather model dependent and small variations may affect the tracks. We approached this problem in a more empirical way: we calculate tracks assuming three different initial radii and see a posteriori which is the impact of the assumption of the initial radius on the evolutionary tracks. As initial radii for each planet, we assume the values obtained by setting Λ equal to 3, 5, and 10.
For latetype stars, the stellar XUV flux depends on the stellar mass and rotation period, where the latter is timedependent. In this work, we have assumed that the rotation period varies with time as (Mamajek & Hillenbrand 2008) (26)
where P_{rot} is the rotation period (in days), τ is the stellar age (in Myr), and is the reddeningfree stellar colour. Equation (29) represents the average approximation based on a large set of latetype dwarfs, but in reality the rotation tracks of stars are nonunique, which leads to different evolutionary tracks of the XUV radiation (Johnstone et al. 2015b; Tu et al. 2015). In the most general case, this results in a variety of planetary atmosphere evolution tracks, but for the present analysis that is limited to small, closein planets we restrict ourselves to the approximation given by Eq. (29).
We set the relevant stellar parameters as follows. The stellar Xray flux can then be inferred from the rotation period as (Wright et al. 2011) (27)
where L_{X} is the Xray stellar luminosity, L_{bol} is the bolometric luminosity, C = 8.68 × 10^{−6} and β = −2.18 are empirical constants (Pizzolato et al. 2003; Wright et al. 2011), and R_{0,sat} = 0.13 is the saturation threshold corresponding to the Rossby number R_{0}. R_{0} is the ratio between the stellar rotation period and the convective turnover time (T_{conv}; Wright et al. 2011) (28)
The EUV stellar luminosity is then derived from the Xray luminosity using Eq. (22). To account for how the bolometric luminosity and equilibrium temperature change with time, we employ the MESA Isochrones and Stellar Tracks (MIST; Paxton et al. 2018).
We began the evolution at the age of 5 Myr, which is approximately the typical lifetime of protoplanetary disks (Mamajek 2009). However, since we are, in general, interested in Gyrsold planets, the exact initial time for theevolution has no significant effect on the results. Having set the initial planetary radius, f_{at}, and XUV stellar flux at the planetary orbital separation, we extracted the massloss rate from our grid, which we then used to derive how much mass is lost during the first time step. At this point, we derive the new planetary atmospheric mass that we converted into a planetary radius by interpolating over the f_{at} grid, and began the cycle again, updating at each time step the stellar XUV flux. We finally obtained atmospheric evolutionary tracks by choosing smallenough time steps, which in our case adapt to the massloss rates by ensuring that the maximum mass loss within one time step is smaller than 1% of the planetary mass, and by repeating this procedure up to the desired age (e.g. age of the given system) or till the planetary radius has reached the core radius. We ignored gravitational contraction and radioactive decay that contribute to increase the equilibrium temperature during the first phases of evolution.
5.2.2 Results for CoRoT7 b
CoRoT7 b has a mass of 5.74 M_{⊕} and a radius of 1.585 R_{⊕}, which indicates a rocky composition and a lack of a hydrogendominated envelope (Léger et al. 2009; Mura et al. 2011; Barros et al. 2014). The planet orbits an active early Ktype star (M_{*} = 0.93 M_{⊙}, R_{*} = 0.87 R_{⊙}, T_{eff} = 5275 K) at a distance of 0.0172 AU, which corresponds to a period of 0.853 days. The planet has an equilibrium temperature of 1756 K and the age of the system has been estimated to be 1.5±0.3 Gyr.
Figure 10 shows the evolution of the planetary radius as a function of time obtained for each of the three tested initial conditions corresponding to planetary radii of 8.25 R_{⊕} (f_{at} ≈ 6 × 10^{−2}), 4.95 R_{⊕} (f_{at} ≈ 2 × 10^{−2}), and 2.47 R_{⊕} (f_{at} ≈ 6 × 10^{−4}). The two tracks starting with the largest planetary radius converge quickly to the same point, implying that the result is independentof the initial condition. This is because the first part of the evolution is dominated by boiloff. The track starting with the smaller planetary radius instead does not present clear signs of a boiloff phase, but still leads to a rapid complete escape of the atmosphere. The tracks indicate that the planet is supposed to have completely lost its hydrogendominated envelope within an extremely short time of about 0.1 Myr.
The most significant changes in the size of the atmosphere occur during the first 10^{−2}−10^{−1} Myr, when the atmosphere lies in the boiloff regime, which is consistent with what was found by Owen & Wu (2016). Once the radius has reached about 2 R_{⊕}, the escape is driven by the stellar XUV flux, which for CoRoT7 b is very intense since the system is still rather young and the planet has a very short orbital period. The complete escape of the atmosphere is so fast that the initial stellar rotation rate (see, e.g. Tu et al. 2015) does not play a significant role. We therefore estimate that CoRoT7 b has lost its primary atmosphere, assuming it had accreted one to begin with, within a maximum time of about 0.1 Myr.
Atmospheric escape for CoRoT7 b has been previously studied by Jackson et al. (2010) and Leitzinger et al. (2011). Both inferred the planet’s massloss rate over time using the energylimited approximation, accounting for the Roche lobe effect (Erkaev et al. 2007), which greatly underestimates the massloss rates at the beginning of the planet’s evolution.
Jackson et al. (2010) consider the effects of the possible planetary migration through orbital tidal decay, a wide range of initial radii (up to a gas giant), an effective radius of 3 R_{pl}, various heating efficiencies up to 100%, and applied the scaling laws for the stellar XUV flux of Ribas et al. (2005). They arrive at the conclusion that CoRoT7 b could have started its evolution as a gas giant, with a mass of up to 200 M_{⊕}. In case CoRoT7 b has, instead, always been a rocky planet, they suggested that it could have lost up to half of its mass through surface melting, outgassing, and subsequent escape of the secondary atmosphere.
Leitzinger et al. (2011) do not account for planetary migration, but employed more realistic heating efficiencies of 10–25% and stellar irradiation levels consistent with those adopted in our work. Their calculations lead them to exclude that CoRoT7 b had started its evolution as a gas giant planet, with a mass similar or larger than that of Saturn, otherwise the planet would still host a significant hydrogendominated envelope, which is excluded by the Earthlike bulk density.
At the very beginning of the evolution, when the planetary atmosphere is supposed to be in boiloff, the massloss rates derived from the grid are larger than those considered by Jackson et al. (2010) and Leitzinger et al. (2011) by a factor of 10–10^{6}, depending on the initial planetary radius. This large difference is caused by the fact that the energylimited approximation is not capable of describing atmospheric escape in the boiloff phase. In the blowoff phase, instead, the massloss rates derived from the grid are about a factor of two smaller than those of Leitzinger et al. (2011) and a factor of a few smaller than those of Jackson et al. (2010).
Following the works of Jackson et al. (2010) and Leitzinger et al. (2011), we further tested the possible evolution of CoRoT7 b by increasing even more the initial planetary mass (and radius) obtaining that the planet needed to have a mass smaller than that of Uranus (about 14.5 M_{⊕}) to loose the primary atmosphere within the age of the system. An initial CoRoT7 b with a mass equal to that of Neptune (about 17 M_{⊕}) would now still hold a hydrogendominated envelope with f_{at} = 0.2. We did not explore an even heavier starting point because of the upper mass limit of 39 M_{⊕} in our grid.
Fig. 10 Evolution of the planetary radius of CoRoT7 b as a function of time. The colours indicate different initial radii, marked by the asterisks, which correspond to the values obtained by setting Λ = three (red), five (blue), and ten (green). The small dots placed along each line indicate the time steps. 
5.2.3 Results for HD 219134 b,c
The two closein, transiting superEarths HD 219134 b and c orbit a K3 mainsequence star with a radius of 0.778 R_{⊙}, a mass of 0.81 M_{⊙}, and an effective temperature of 4699 K. The estimated age of the system is 11 ± 2 Gyr. The two planets have measured masses of 5.74 and 4.74 M_{⊕}, and radii of 1.602 and 1.511 R_{⊕}, respectively.Therefore, both planets present Earthlike densities. They orbit the host star at distances of 0.039 and 0.065 AU, respectively.
The evolutionary tracks, shown in Fig. 11, indicate that the primary, hydrogendominated atmosphere escaped completely within about 12 and 80 Myr for HD 219134 b and c, respectively. The difference in time between the two planets for a complete atmospheric escape is due to their different distance from the star, thus different equilibrium temperature and stellar XUV irradiation. These times are significantly shorter than the estimated age of the system, thus allowing us to conclude that both planets have most likely completely lost their primary, hydrogendominated atmosphere through escape. This is in agreement with Dorn & Heng (2018), who arrived at the same conclusion by employing a Bayesian inference method based on the stellar properties and the energylimited approximation. However, these are extreme cases, and an approach based on the energylimited formulation would most likely lead to the wrong results for younger and/or lower density planets.
Fig. 11 Evolution of the planetary radii of HD 219134 b (top panel) and HD 219134 c (bottom panel) as a function of time. Colours, symbols, and lines are as in Fig. 10. 
6 Conclusion
We upgraded and employed an existing planetary upper atmosphere hydrodynamic code to compute a large grid of models for superEarths and miniNeptunes orbiting latetype stars. The main upgrade consists in the implementation of a scheme that automatically sets the initial parameters and profiles for each run, thus in practice allowing one to automate computations. The planets covered by the grid have masses ranging from 1 to 39 M_{⊕} and orbit early M to late Ftype stars in a wide range of orbital distances, corresponding to equilibrium temperatures between 300 and 2000 K. For each of the considered stellar masses, we have also considered three different values of the XUV flux. The wide parameter space covered by the grid allowed us to model a broad variety of planetary atmospheres, ranging from being in boiloff, to blowoff, and to very stable atmospheres.
For each planet in the grid, we computed the atmospheric temperature, number density, bulk velocity, Xray and EUV volume heating rates, and abundance of the considered species as a function of distance from the planetary centre. From these quantities, we estimated the positions of maximum dissociation and ionisation, the massloss rate, and the effective radius of the XUV absorption.
We compared the results of our grid, in particular the massloss rates, with those previously published for planets inside and outside the grid boundaries, finding excellent agreement. We also developed a tool to interpolate among the model results to infer the atmospheric properties of any planet covered by the grid. We took advantage of the large grid to explore in detail how the atmospheric characteristics vary with system parameters, finding for example that the massloss rate can be analytically described as a loglinear function of Λ and a linear function of the stellar XUV flux.
Our grid and the interpolation routine allow one to extract in a fraction of a second information that would otherwise require days or weeks to obtain. This enables one to employ the results of proper hydrodynamic computations of the massloss rates in planetary atmospheric evolution calculations. This avoids the need to use approximations, such as the energylimited formula, that have been shown (by various authors and further in this work) to significantly underestimate or overestimate in some cases the massloss rates.
We have therefore applied our grid and interpolation routine to study the evolution of the closein, highdensity planets CoRoT7 b and HD 219134 b,c. For CoRoT7 b, we found that the primary hydrogendominated atmosphere, assuming the planet has ever accreted one, was lost mostly through boiloff within about 0.1 Myr. We also concluded that the planet originally could have been as massive as Uranus, because the envelope of a heavier planet would have been too massive to completely escape within the age of the system. We arrived at a similar conclusion also for HD 219134 b,c, where for these two planets we found that the hydrogendominated atmosphere escaped completely within about 12 and 80 Myr, respectively. It is therefore likely that other similar planets, such as Kepler10b and 55 Cnc e, followed an analogous evolutionary path, which have left them with a secondary atmosphere formed by either outgassing from the magma ocean or sputtering of the stellar wind on the bare planetary surface, similar to what happens for Mercury (Mura et al. 2011; Guenther et al. 2011; Pfleger et al. 2015; Vidotto et al. 2018).
The simple evolutionary tracks we computed for CoRoT7 b and HD 219134 b,c provide just an example of what can be achieved using the grid and interpolation routine. There is, however, a large number of applications in which these tools can be used and we will explore a few of them(e.g. analytic formulation of the massloss rates as a function of system parameters) in future works currently in preparation. We are still working on increasing the size of the grid, extending it towards more massive planets and towards less massive stars. The grid and interpolation routine can be downloaded online^{1}. We are, however, planning for the near future to set up a web interface allowing users to query the grid and run the interpolation routine more easily.
Acknowledgements
We acknowledge the Austrian Forschungsförderungsgesellschaft FFG project “TAPAS4CHEOPS” P853993, the Austrian Science Fund (FWF) NFN project S11607N16, the FWF project P27256N27 and the FWF project P30949N36. N.V.E. acknowledges support by the RFBR grant No. 180500195a and 165214006 ANF_a. We thank the anonymous referee for the positive approach and the useful comments that led to a significant improvement of the manuscript.
Appendix A List of reactions and crosssections employed in the model
Reactions and crosssections employed in the model are presented in Table A.1.
Reactions and relative crosssections employed in the model.
Appendix B Normalisations employed for the computation of each model
Here we give the normalisations used for the computation of the models.
In the equations above, the subscript “0” denotes the values at the lower boundary, for example, P_{0} and C_{s0} are, respectively, the pressure and sound speed at the lower boundary.
Appendix C Planetary atmospheric massloss rates as a function of system parameters
We present here plots analogous to those in Fig. 3, but for different stellar masses (Figs. C.1–C.4).
Fig. C.4 Same as Fig. 3, but for a stellar mass of 0.4 M_{⊙}. In these plots, the temperature range was shifted to 300–1100 K, instead of 700–1500 K, to make the range of orbital separations comparable to that of other stellar masses and because of the cut on the Roche lobe for planets with an equilibrium temperature of 1500 K orbiting 0.4 M_{⊙} stars. 
References
 Barros, S. C. C., Almenara, J. M., Deleuil, M., et al. 2014, A&A, 569, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Stevenson, D. J. 2013, ApJ, 769, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Black, J. H. 1981, MNRAS, 197, 553 [NASA ADS] [Google Scholar]
 Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., & Lecavelier des Etangs, A. 2013, A&A, 557, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, Eur. Phys. J. Web Conf., 47, 03005 [CrossRef] [EDP Sciences] [Google Scholar]
 CecchiPestellini, C., Ciaravella, A., Micela, G., & Penz, T. 2009, A&A, 496, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chamberlain, J. W. 1963, Planet. Space Sci., 11, 901 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017a, MNRAS, 466, 1868 [NASA ADS] [CrossRef] [Google Scholar]
 Cubillos, P. E., Fossati, L., Erkaev, N. V., et al. 2017b, ApJ, 849, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Deming, D., Seager, S., Winn, J., et al. 2009, PASP, 121, 952 [NASA ADS] [CrossRef] [Google Scholar]
 Dorn, C., & Heng, K. 2018, ApJ, 853, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Ehrenreich, D., & Désert, J.M. 2011, A&A, 529, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Lecavelier Des Etangs, A., Hébrard, G., et al. 2008, A&A, 483, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Erkaev, N. V., Lammer, H., Odert, P., et al. 2013, Astrobiology, 13, 1011 [NASA ADS] [CrossRef] [Google Scholar]
 Erkaev, N. V., Lammer, H., ElkinsTanton, L. T., et al. 2014, Planet. Space Sci., 98, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Erkaev, N. V., Lammer, H., Odert, P., Kulikov, Y. N., & Kislyakova, K. G. 2015, MNRAS, 448, 1916 [NASA ADS] [CrossRef] [Google Scholar]
 Erkaev, N. V., Lammer, H., Odert, P., et al. 2016, MNRAS, 460, 1300 [NASA ADS] [CrossRef] [Google Scholar]
 Erkaev, N. V., Odert, P., Lammer, H., et al. 2017, MNRAS, 470, 4330 [NASA ADS] [CrossRef] [Google Scholar]
 Fleming, B. T., France, K., Nell, N., et al. 2018, J. Astron. Telesc. Instrum. Syst., 4, 014004 [NASA ADS] [CrossRef] [Google Scholar]
 Fossati, L., France, K., Koskinen, T., et al. 2015, ApJ, 815, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Fossati, L., Erkaev, N. V., Lammer, H., et al. 2017, A&A, 598, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fossati, L., Koskinen, T., France, K., et al. 2018, AJ, 155, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., & Petigura, E. A. 2018, AJ, accepted [arXiv:1805.01453] [Google Scholar]
 Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Demory, B.O., Van Grootel, V., et al. 2017, Nat. Astron., 1, 0056 [NASA ADS] [CrossRef] [Google Scholar]
 Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Guenther, E. W., Cabrera, J., Erikson, A., et al. 2011, A&A, 525, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guo, J. H., & BenJaffel, L. 2016, ApJ, 818, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Howe, A. R., Burrows, A., & Verne, W. 2014, ApJ, 787, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, B., Miller, N., Barnes, R., et al. 2010, MNRAS, 407, 910 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, A. P., Davis, T. A., & Wheatley, P. J. 2012, MNRAS, 422, 2024 [NASA ADS] [CrossRef] [Google Scholar]
 Jeans, J. 1925, The Dynamical Theory of Gases. By Sir James Jeans (Cambridge: Cambridge University Press), 1925 [Google Scholar]
 Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015a, ApJ, 815, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015b, A&A, 577, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kislyakova, K. G., Johnstone, C. P., Odert, P., et al. 2014, A&A, 562, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koskinen, T. T., Yelle, R. V., Lavvas, P., & Lewis, N. K. 2010, ApJ, 723, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013, Icarus, 226, 1678 [NASA ADS] [CrossRef] [Google Scholar]
 Kubyshkina, D., Lendl, M., Fossati, L., et al. 2018, A&A, 612, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121 [NASA ADS] [CrossRef] [Google Scholar]
 Lammer, H., Odert, P., Leitzinger, M., et al. 2009, A&A, 506, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lammer, H., Erkaev, N. V., Odert, P., et al. 2013, MNRAS, 430, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Lammer, H., Erkaev, N. V., Fossati, L., et al. 2016, MNRAS, 461, L62 [NASA ADS] [CrossRef] [Google Scholar]
 Lecavelier des Etangs, A., VidalMadjar, A., McConnell, J. C., et al., 2004, A&A, 418, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leitzinger, M., Odert, P., Kulikov, Y. N., et al. 2011, Planet. Space Sci., 59, 1472 [NASA ADS] [CrossRef] [Google Scholar]
 Léger, A., Rouan, D., Schneider, J., et al. 2009, A&A, 506, 287 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Lopez, E. D. 2017, MNRAS, 472, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, Nat. Commun., 7, 11201 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E. 2009, Am. Inst. Phys. Conf. Ser., 1158, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Masuda, K. 2014, ApJ, 783, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Menager, H., Barthélemy, M., Koskinen, T., et al. 2013, Icarus, 226, 1709 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, S., Stallard, T., Tennyson, J., & Melin, H. 2013, J. Phys. Chem. A, 117, 9770 [CrossRef] [Google Scholar]
 Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2015, ApJS, 217, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Mura, A., Wurz, P., Schneider, J., et al. 2011, Icarus, 211, 1 [NASA ADS] [CrossRef] [Google Scholar]
 MurrayClay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Öpik, E. J. 1963, Geophys. J., 7, 490 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Mohanty, S. 2016, MNRAS, 459, 4088 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Wu, Y. 2016, ApJ, 817, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfleger, M., Lichtenegger, H. I. M., Wurz, P., et al. 2015, Planet. Space Sci., 115, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quirrenbach, A., Amado, P. J., Mandel, H., et al. 2010, Proc. SPIE, 7735, 773513 [CrossRef] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Reiners, A., Schüssler, M., & Passegger, V. M. 2014, ApJ, 794, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [NASA ADS] [CrossRef] [Google Scholar]
 Rodríguez, A., Callegari, N., & Correia, A. C. M. 2016, MNRAS, 463, 3249 [NASA ADS] [CrossRef] [Google Scholar]
 Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2016, A&A, 585, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SanzForcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shematovich, V. I., Ionov, D. E., & Lammer, H. 2014, A&A, 571, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spitzer, L. 1978, Physical Processes in the Interstellar Medium (New York: WileyInterscience), 333 [Google Scholar]
 Stökl, A., Dorfi, E., & Lammer, H. 2015, A&A, 576, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stökl, A., Dorfi, E. A., Johnstone, C. P., & Lammer, H. 2016, ApJ, 825, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Tinetti, G., Drossart, P., Eccleston, P., et al. 2017, European Planetary Science Congress, 11, EPSC2017713 [Google Scholar]
 Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2018, MNRAS, 479, 4786 [NASA ADS] [CrossRef] [Google Scholar]
 Vidotto, A. A., Lichtenegger, H., Fossati, L., et al. 2018, MNRAS, 481, 5296 [NASA ADS] [CrossRef] [Google Scholar]
 Vogt, S. S., Burt, J., Meschiari, S., et al. 2015, ApJ, 814, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Yelle, R. V. 2004, Icarus, 170, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Yi, S., Demarque, P., Kim, Y.C., et al. 2001, ApJS, 136, 417 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Comparison between the massloss rates obtained from our hydrodynamic modelling (Col. 6), from the energylimited formula (Col. 7), and from the literature (Col. 8).
Stellar masses, equilibrium temperatures, planetary radii, and planetary masses considered for the computation of the grid.
All Figures
Fig. 1 Position of some of the modelled planets (black crosses) in the massradius diagram. All planets orbit a 1 M_{⊙} star and have an equilibrium temperature of 300 K (top panel) and 2000 K (bottom panel). The blue and magenta solid lines indicate the boundaries of the grid set by the cut on Λ and on bulk density, respectively. Since Λ depends on T_{eq}, the top boundary is different in the different panels. The green solid line indicates the boundary of the grid set by the cut on the Roche lobe. The position of this boundary depends on the orbital separation, thus on T_{eq}. For reference, the red dashed line indicates Earth’s density. 

In the text 
Fig. 2 From top to bottom: atmospheric profiles for density, temperature, velocity, and fraction of atomic hydrogen for planets with Λ equal to 4.8 (red) and 66.7 (blue). Both planets orbit a 1 M_{⊙} star, have T_{eq} = 1100 K and R_{pl} = 3 R_{⊕}, and are subject to an incident XUV flux of 92.6 erg cm^{−2} s^{−1}. The planet with the lower Λ has a mass of 2.1 M_{⊕}, while the other one has a mass of 29.1 M_{⊕}. The density (top panel) is normalised to its value at R_{pl}. The blue and red dashed lines show the effective radii of the XUV absorption. To allow comparing planets with significantly different Roche radii, the xaxis is in units of the planetary Roche radii, starting from the planetary surface. 

In the text 
Fig. 3 Logarithm of the massloss rates (colour coded) as a function of planetary mass and radius. The adopted T_{eq} and F_{XUV} values are given on the top of each panel. The equilibrium temperature increases from top to bottom, while F_{XUV} increases from left to right. All planets orbit a 1.0 M_{⊙} star. For reference, the dashed lines mark constant Λ values of 8, 20, and 50 (from top to bottom). The red lines indicate planets for which the Knudsen number at the upper boundary is equal to one (solid line) and ten (dashed line). 

In the text 
Fig. 4 Ratio of the massloss rates computed with the hydrodynamic model to the energylimited formula as a function of Λ (top panel) and M_{pl} (bottom panel). The green lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 300 K, F_{XUV} = 159.4 erg cm^{−2} s^{−1}. The blue lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 700 K, F_{XUV} = 4900 erg cm^{−2} s^{−1}. The red lines and circles are for systems with the following characteristics: M_{*} = 1.0 M_{⊙}, T_{eq} = 1100 K, F_{XUV} = 30784 erg cm^{−2} s^{−1}. The horizontal dotted line indicates the equality between the two values. 

In the text 
Fig. 5 Relative deviation of the interpolated massloss rates (Ṁ_{int}) from the computed ones (Ṁ) as a functionof Λ. The horizontal red lines indicate deviations of 1%, 5%, and 10%. 

In the text 
Fig. 6 Atmospheric massloss rate Ṁ as a functionof Λ for all computed planets. The colour code indicates the planetary Roche radius in units of R_{pl}. The position of the test planets listed in Table 4 is shown by black squares (Pa1 and Pa2), black circles (Pb1 and Pb2), purple squares (Pc1 and Pc2), and purple circles (Pd1 and Pd2). The lines indicate the predictions obtained by using the energylimited formula for the Pb test planets varying planetary mass only and assuming the value of R_{eff} equal to the planetary radius (dashdotted line) or equal to the value derived from the grid (dashed line). The black crosses and plus signs indicate the escape rates estimated for some of the known transiting exoplanets using the interpolation routine and direct hydrodynamic calculations, respectively. 

In the text 
Fig. 7 Effective radius (black), position of the maximum dissociation (blue), and position of the maximum ionisation (red) as a functionof Λ, for all planets in the grid. 

In the text 
Fig. 8 Effective radius R_{eff} as a function of R_{roche}. Each point is colourcoded with respect to the orbital semimajor axis. The black and purple circles and squares indicate the position of the eight test planets following the same symbols as in Fig. 6. The circle and square purple symbols largely overlap. 

In the text 
Fig. 9 Outflow velocity at the Roche radius V_{roche} as a functionof Λ. Each point is colour coded with the value of the orbital semimajor axis. The black and purple circles and squares indicate the position of the eight test planets following the same symbols as in Fig. 6. 

In the text 
Fig. 10 Evolution of the planetary radius of CoRoT7 b as a function of time. The colours indicate different initial radii, marked by the asterisks, which correspond to the values obtained by setting Λ = three (red), five (blue), and ten (green). The small dots placed along each line indicate the time steps. 

In the text 
Fig. 11 Evolution of the planetary radii of HD 219134 b (top panel) and HD 219134 c (bottom panel) as a function of time. Colours, symbols, and lines are as in Fig. 10. 

In the text 
Fig. C.1 Same as Fig. 3, but for a stellar mass of 1.3 M_{⊙}. 

In the text 
Fig. C.2 Same as Fig. 3, but for a stellar mass of 0.8 M_{⊙}. 

In the text 
Fig. C.3 Same as Fig. 3, but for a stellar mass of 0.6 M_{⊙}. 

In the text 
Fig. C.4 Same as Fig. 3, but for a stellar mass of 0.4 M_{⊙}. In these plots, the temperature range was shifted to 300–1100 K, instead of 700–1500 K, to make the range of orbital separations comparable to that of other stellar masses and because of the cut on the Roche lobe for planets with an equilibrium temperature of 1500 K orbiting 0.4 M_{⊙} stars. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.