Free Access
Issue
A&A
Volume 612, April 2018
Article Number A25
Number of page(s) 10
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201731816
Published online 16 April 2018

© ESO 2018

1 Introduction

With the number of extrasolar planets that we have discovered rapidly increasing and the development of characterization techniques, the comparative analysis of exoplanet properties is becoming possible. To date, this effort has focused largely on the comparison of planets of different masses and stellar irradiation levels. The temporal component of this complex question, that is, the evolution of planetary atmospheres with time, has been touched upon mostly in the context of planets orbiting Gyrs old stars, evolved away from the zero-age main sequence, and the reaction of planetary atmospheres to slowly increasing irradiation caused by a host star evolving off the main sequence (e.g. Owen & Wu 2016; Lopez & Fortney 2016). This is mostly because the early evolution of planetary atmospheres, which is the most relevant in terms of atmospheric evolution (e.g. Owen & Wu 2013; Lopez & Fortney 2013; Jin et al. 2014; Chen & Rogers 2016), has been inaccessible due to various observational difficulties related to the activity and variability of young stars that have precluded the detection of planets in their early phases of evolution.

After the dispersal of the protoplanetary disk, extrasolar planets cool and contract to their currently observed sizes (e.g. Baraffe et al. 2006; Stökl et al. 2016). Recently, Stökl et al. (2015) and Owen & Wu (2016) suggested that planets emerging from their disks rapidly find themselves in a state of strong mass loss, called “boil-off”, due to Parker winds in the planetary atmosphere driven by the stellar thermal irradiation. During this short lived phase, the planet rapidly cools and contracts, losing a substantial fraction of its gaseous envelope before the planetary atmosphere settles into a more stable configuration and continues to evolve on slower timescales (e.g. Jin et al. 2014; Fossati et al. 2017; Jin & Mordasini 2017; Owen & Wu 2017). For planets close enough to the host star, the atmosphere may then go through a long phase of efficient hydrodynamic escape, or even “blow-off”, which, depending on the planet’s atmosphere mass, may or may not significantly affect the planet (e.g. Lammer et al. 2003; Yelle 2004; Lecavelier des Etangs et al. 2004).

In this work, we present simulations of the upper atmospheric structure of the youngest confirmed exoplanet, K2-33b. Identified from K2 photometry by Vanderburg et al. (2016) and determined to be a pre-main-sequence low-mass member of the Upper Scorpius OB association by David et al. (2016) and Mann et al. (2016), K2-33b is a planet with a radius of about 5 R (i.e. 1.3 times Neptune’s radius). The semi-major axis is 0.0409 au and the rotational period is 5.43 days. The planetary mass is undetermined, yet a conservative upper limit of 5.4 MJ (Mann et al. 2016) indicates that the object is clearly a planet. Owing to the planet’s size, the true mass is likely to be somewhat smaller than that of Jupiter (i.e. ≈318 M). The K2-33 system still possesses a protoplanetary disk; however, the lack of emission from hot material indicates that the disk appears to be truncated at 2 AU (David et al. 2016). Age estimates for the K2-33 system are 5–10 Myr (David et al. 2016) and <20 Myr (Mann et al. 2016).

The planet K2-33b is thus an object that has only very recently emerged from its birthplace inside the protoplanetary disk in which it was formed, and its atmosphere may still be in the process of cooling and contracting to a more stable configuration. The planet may either have formed further out in the system and travelled to its present location via disk-driven migration (e.g. Lin et al. 1996; Ida & Lin 2008; Rafikov 2006; Schlaufman et al. 2009; Lubow & Ida 2010; Raymond & Cossou 2014), or have formed at (or very near) its present location (Chiang & Laughlin 2013; Ogihara et al. 2015). Owing to the young age of the system, the planet is likely experiencing extremely high levels of high-energy X-ray and extreme ultraviolet (XUV) stellar irradiation (e.g. Wright et al. 2011). Its youth makes this object a key target for atmospheric studies, providing a snapshot on the earliest stages of planetary evolution.

We present here a grid of upper atmosphere models for K2-33b, aiming to determine the possible state of this planet’s atmosphere and to look for possible constraints on the planet’s mass from upper atmosphere modelling. In Sect. 2 we present the physical model used in our calculations, the results of which are presented in Sect. 3. In Sect. 4 we offer a first interpretation of the possible atmospheric properties of K2-33b and summarize our conclusions in Sect. 5.

2 Modelling setup and formalism

2.1 System parameters

David et al. (2016) and Mann et al. (2016) estimate compatible values of 1.1 ± 0.1 and 1.05 ± 0.07 R for the radius of the host star K2-33. From the rather well determined stellar radius follows a planetary radius of 5.76 R (David et al. 2016), or 5.04 R (Mann et al. 2016). However, stellar mass estimates diverge depending on whether or not the stellar density, inferred from the transit light curve (following Seager & Mallén-Ornelas 2003), is included in the analysis. By considering stellar densities inferred from the transit light curve (0.49 ρ), Mann et al. (2016) find a stellar mass of 0.56 ± 0.09 M, while David et al. (2016) report a mass of 0.31 ± 0.05 M and a corresponding stellar density of 0.34 ± 0.12 ρ, based on stellar evolutionary models. For our calculations, we base ourselves on Mann et al. (2016) and use therefore a value of 0.56 ± 0.09 M. We also consider an orbital separation of d = 0.04 au.

We estimated the stellar X-ray luminosity from the basic stellar parameters: mass, radius, and rotation rate. The X-ray luminosity, LX, can be calculated from (e.g. Wright et al. 2011)

where R0,sat is the saturation threshold corresponding to the Rossby number R0, which is the ratio between the stellar rotation period and the convective turnover time (R0,sat = 0.13; Wright et al. 2011). The X-ray emission reaches the saturation value RX,sat for the fast rotating stars, with R0 smaller than R0,sat. Mann et al. (2016) estimated the stellar rotation period to be 6.27 ± 0.17 days and the theoretical prediction of the convection turnover time for a stellar mass of 0.5 M is 35 days (Wright et al. 2011). This implies that the star is already out of the saturated regime. By using values of β = −2.18 and C = 8.68 × 10−6 (Pizzolato et al. 2003; Wright et al. 2011), we obtain LX = 2.16 × 1029 erg s−1 and consequently an X-ray flux at the planetary orbit of 48 257 erg s−1 cm−2, which is about 105 times stronger than that currently experienced at Earth.

The X-ray and extreme ultraviolet (EUV) luminosities are related as (Sanz-Forcada et al. 2011) (1)

which leads to a stellar EUV luminosity of LEUV = 1.066 ×1030 erg s−1 and an EUV flux at the planet’s orbit of 2.382 ×105 erg s−1 cm−2.

Since the planetary mass is essentially unknown, we perform our simulations over a grid in planetary mass, assuming values of 2, 4, 7, 10, 20, and 40 M. We also consider different values of the equilibrium temperature (Teq) of 850, 1000, 1150, and 1300 K, where the lowest value of 850 K is roughly the Teq value calculated on the basis of the system parameters (Mann et al. 2016), assuming a Bond albedo of zero and zero energy redistribution. By considering higher temperatures, we account for possible additional internal heating originating from gravitational contraction and radioactive decay.

2.2 Physical model

To model the planetary upper atmosphere and infer the mass-loss rates, we employ an updated version of the 1D hydrodynamical code described in detail in the appendix of Erkaev et al. (2016). This model treats the planetary upper atmosphere as a clear hydrogen gas envelope, including the effects of recombination, dissociation, ionization of the hydrogen atoms and molecules, and Ly-α cooling. In this work, we additionally include cooling from H3+ molecules (see details below).

The original code does not account for the spectral dependence of stellar high-energy flux and considers just the integrated EUV flux, excluding X-rays. To make the model more appropriate for treating young stars with intense X-ray fluxes, such as K2-33, we implemented also X-ray heating (see details below).

The lower boundary of the calculation domain, Rpl, is set to be at the planetary photospheric radius, where the optical and infrared stellar photons are absorbed. The upper boundary is set at the Roche radius , where , and Mpl and M* are the planetary and stellar masses, respectively. We also assume that the atmosphere at the lower boundary consists exclusively of molecular hydrogen, while the rest of the atmosphere is composed of molecular and atomic hydrogen, their ions, and H3+.

For each modelled planet, the atmospheric pressure at the planetary photospheric radius (i.e. the lower boundary) was estimated employing the same technique described in Fossati et al. (2017) and Cubillos et al. (2017), which is based on the radiative transfer code BART (Bayesian Atmospheric Radiative Transfer; Cubillos 2016), and assuming solar metallicity. The obtained pressure values, which range between about 40 and 400 mbar, are listed in Table 1. We further set the temperature at the lower boundary to be equal to the equilibrium temperature (see Fossati et al. 2017, for a discussion on this point).

The height-averaged heating efficiency (η) of the upper atmosphere exposed to high-energy stellar radiation is usually considered to be between 10% and 60% (e.g. Watson et al. 1981; Yelle 2004; Murray-Clay et al. 2009; Cecchi-Pestellini et al. 2009; Owen & Jackson 2012; Shematovich et al. 2014; Salz et al. 2016). In this work, we follow the considerations of Erkaev et al. (2016) and adopt a value of 15%, based on the results of the direct simulation Monte Carlo model calculations of Shematovich et al. (2014), which model photolytic and electron impact processes in the thermosphere by solving the kinetic Boltzmann equation.

With the inclusion of Ly-α and H3+−cooling, the equations for mass, momentum, and energy conservation are respectively

Here ρ, v, T, and P are mass density, velocity, temperature, and pressure of the gas as a function of the radial distance r from Rpl, respectively1. The thermal energy E is equal to . The thermal conductivity is χ, which is given by (Watson et al. 1981), where the exponent 0.7 was chosen for the neutral gas (Banks & Kockarts 1973). The terms QXUV and QLy − α are the stellar XUV volume heating rate and Lyα cooling, respectively, while is the gravitational potential.

In the absence of a measured stellar XUV spectrum and to speed up computations for the XUV stellar flux absorption, we follow the approach of Erkaev et al. (2016), who considered that the whole stellar EUV flux is emitted at a single wavelength of 60 nm, which corresponds to the position of the peak of the stellar EUV emission. To account for X-ray heating, we additionally consider absorption of stellar X-ray photons by the planetary atmosphere. Here, we assume that they are all emitted at a wavelength of 5 nm, which lies roughly in the middle of the X-ray wavelength range. To summarize, we set the integrated X-ray and EUV stellar flux to be emitted at 5 and 60 nm, respectively. The XUV heating function QXUV is therefore composed of two terms, QEUV and QX, which describe the heating by the EUV and X-ray stellar flux, respectively. These two functions are constructed in the same way, but the absorption cross-sections and absorption functions of the stellar flux inside the planetary atmosphere are defined at the two wavelengths given above. This results in a redistribution of the energy absorption inside the atmosphere, which, however, does not significantly affect the large-scale characteristics of the atmosphere.

The total heating function is QXUV = QEUV + QX. Each heating function is given by (5)

where m stands for either EUV or X, σm is the absorption cross-section for the specific wavelength, and ϕm is the flux absorption function, given by (6)

Here, Jm(r, θ) is a function in spherical coordinates describing the spatial variation of the EUV, or X-ray, flux due to atmospheric absorption (Erkaev et al. 2015) and r in this case is radial distance from the planetary centre, in units of Rpl.

Following Murray-Clay et al. (2009), the absorption cross-sections are given by , where Ei is the hydrogen ionization energy and Eλ is the photon energy at a specific wavelength range, that is Eλ = 20 eV in the EUV and 248 eV in the X-ray domain. This gives EUV flux absorption cross-sections of 2 × 10−18 cm−3 and 1.2 × 10−18 cm−3 for atomic and molecular hydrogen, respectively. The X-ray absorption cross-section is σx ~ 0.5 × 10−3σeuv.

We calculate Ly-α cooling using (Watson et al. 1981) (7)

We also include H3+ cooling, as given by Miller et al. (2013), in the form of (8)

where Cn are the temperature-dependent coefficients listed in Table 5 of Miller et al. (2013). To include H3+ in the model, we add the following additional chemical reactions (with rate coefficients taken from Yelle 2004)

both having a cross-section of γH2 = 2 × 10−9,

with a cross-section of , and

with a cross-section of . By including these reactions, the continuity equations for the atmospheric species become

Assuming quasi-neutrality, the electron density is given by (15)

while the total hydrogen number density is the sum of the number densities of all species. The ionization and dissociation reactions used in the above equations and their cross-sections are listed in Table 2. We note that our photo-electric cross-sections do not account for X-ray ionization (see for example Lorenzani & Palla 2001; Locci et al. 2018, for a more thorough approach), which is, however, orders of magnitude smaller than that caused by EUV irradiation.

The atmospheric mass density ρ and pressure P are given by (16)

and (17)

while normalization and numerical realization are the same as in Erkaev et al. (2016).

Table 1

Dependence of the lower pressure boundary (in mbar) as a function of planetary Teq and mass (Mpl).

Table 2

Reactions and relative cross-sections employed in the model.

3 Results

3.1 Escape rates

For each cell in the model grid, we calculate the properties of the upper atmosphere including temperature, number density, and bulk flow velocity as a function of radial distance. From these quantities, we estimate the hydrodynamic hydrogen escape rate Lhy as the flow of particles through the upper boundary (18)

where n1 and V 1 are the number density and velocity at the Roche radius Rroche. Figure 1 shows the values of Lhy obtained from the simulations as a function of planetary mass and restricted Jeans escape parameter (19) This last quantity corresponds to the value of the Jeans escape parameter calculated for the observed planetary radius and mass for the planet’s equilibrium temperature, and considering atomic hydrogen (Fossati et al. 2017). This parameter is particularly useful because it gives a rough measure of the planetary gravitational potential versus the intrinsic planetary thermal energy (i.e. excluding XUV heating), and hence it allows us to estimate the stability of the atmosphere against its own thermal energy and planetary gravity. The physical properties of the modelled planets and the derived escape rates are listed in Table 3.

As expected, the escape rates increase with decreasing planetary mass. At high mass (i.e. Mpl > 10 M), the escape is driven by the heating caused by the absorption of the stellar XUV flux and the atmosphere is in the blow-off regime. Heating caused by absorption of the stellar X-ray flux contributes just a few percent of the total escape rates. For lower masses, where Lhy increases steeply, the escape is driven mostly by the planetary high thermal energy and low gravity. For these planets, the atmosphere lies in the boil-off regime (see Sect. 4.1).

thumbnail Fig. 1

Hydrodynamic escape rates Lhy as a function of planetary mass (top) and restricted Jeans escape parameter Λ (bottom). Asterisks and circles mark models for which Lhy is higher and lower, respectively, than that estimated from the energy-limited formulation (see Sect. 4.1, for details). The dashed line in the bottom panel denotes Eq. (22). The vertical dotted line indicates the minimum Λ value for which Eq. (22) is valid.

Open with DEXTER
Table 3

Physical properties of the modelled planets and escape and mass-loss rates derived from the simulations.

3.2 Temperature and velocity profiles

The atmospheric temperature profiles are mostly dependent on the planetary mass, while the adopted Teq value has only a small influence (Fig. 2). For Mpl > 10M the models present efficient heating, with pronounced maxima, caused by the absorption of the stellar XUV radiation. For the lowest mass planets, the temperature profile is dominated by adiabatic cooling, driven by the strong hydrodynamic expansion of the atmosphere in the boil-off regime. For planets with intermediate masses (i.e. 7–15 M), the temperature profile presents a minimum close to Rpl, the position and amplitude of which decrease with increasing mass. We discuss the origin of this temperature minimum in Sect. 4.

For atmospheres dominated by adiabatic cooling, the temperature decreases at approximately the same rate throughout the atmosphere, independently of Teq. For the profiles presenting XUV heating in the atmosphere, the value of the maximum temperature increases slightly with decreasing temperature, with a difference of about 250 K between the coolest and hottest temperature maxima.

As stated in Sect. 2, the planetary atmosphere at the lower boundary consists of pure molecular hydrogen. Due to the XUV irradiation, hydrogen dissociates not far from the planetary surface and the upper part of the atmosphere is thus dominated by atomic hydrogen. The position of the dissociation barrier, the thickness of the hydrostatic H2 layers, and the atmospheric temperature mainly depend on the planetary mass. The dissociation appears to be more effective for the more massive planets, where the stellar flux penetrates deeper into the planetary atmosphere. Furthermore, the thickness of the pure molecular hydrogen layers (i.e. where the other species contribute less than 1%) varies with planetary mass and equilibrium temperature. For Teq = 850 K, the thickness of these layers ranges between 1.3 Rpl for the lowest mass planets to 0.04 Rpl for the highest mass planets. By increasing the equilibrium temperature from the lowest to the highest value, the thickness of the molecular hydrogen layer increases by about 30%.

Figure 3 shows the velocity profiles obtained for the modelled planets. For the lowest mass planets, the bulk of the atmosphere moves at a subsonic speed, with the local sound speed being . For the other planets, the upper part of the atmosphere becomes supersonic above the position of the maximum temperature. Similar to the behaviour observed for the temperature profiles, the position at which the atmosphere becomes supersonic gets closer to Rpl as the mass increases. The sensitivity of both temperature and velocity profiles to variations in Teq decreases significantly with increasing mass: Teq only seems to play a role in shaping the atmospheric profiles for the lower mass planets, for which the atmosphere is in the boil-off regime.

thumbnail Fig. 2

Temperature profiles of the modelled planets as a function of planetary mass and equilibrium temperature. For planets with a mass between 7 and 15 M, the temperature profile presents a minimum close to Rpl.

Open with DEXTER
thumbnail Fig. 3

Gas velocity profiles for the planets with masses of 4 (top-left), 7 (top-right), 20 (bottom-left), and 40 (bottom-right) M. The velocity profiles for planets with 10 and 15 M are very similar to those of the 20 M planet. The colours correspond to different values of the equilibrium temperature (blue: 850 K; red: 1000 K; yellow: 1150 K; magenta: 1300 K). The line style indicates whether a certain velocity is below (solid) or above (dash-dotted) the local sound speed.

Open with DEXTER

4 Discussion

4.1 Atmospheric escape regime

We compare here the escape rates derived from the models with those estimated through the energy-limited formulation (Watson et al. 1981; Erkaev et al. 2007) (20)

where FXUV is the stellar XUV flux incident on the planet, Reff is the radius at which most of the XUV flux is absorbed, and K is a term accounting for Roche lobe effects (Table 3). The value of Reff is calculated as (21)

In Eq. (21), JXUV corresponds to the stellar XUV flux inside the atmosphere, which depends on r and on the spherical angle (see Erkaev et al. 2007, 2015, for more details). The Reff values correspond roughly to the positions of the temperature maxima shown in Fig. 2 and hence Reff increases with decreasing Λ. For the lowest mass planets, where the temperature decreases adiabatically, Reff lies instead above the Roche radius, which is the upper boundary of our calculations. For such planets, we therefore consider Reff to be equal to Rroche, hence slightly underestimating the Len values.

Figure 4 shows the Lhy /Len ratio as a function of Λ. This plot indicates that for Λ values ranging between 10 and 40, the two escape rates are roughly similar, while for lower Λ the escape rates calculated using the hydrodynamic model are significantly higher than those derived using Eq. (20). This agrees with the findings of Fossati et al. (2017), who carried out a similar analysis, but for 5 and 10 M planets subject to much smaller amounts (i.e. factor of 10 to 1000) of XUV irradiation. The atmospheres of planets with the lowest Λ values lie therefore in the boil-off regime. This is also shown, for example, by the strong dependence of the escape rates on Teq for these planets (see Fig. 1).

Given the monotonic increase of the escape rates with decreasing Λ , we fit the Lhy values with a logarithmic function of Λ, obtaining (22)

where the fit was done excluding the planets for which Reff lies above Rroche. It is important to remark that this fit is valid just for planets with Λ greater than about five and for the stellar parameters, including the XUV fluxes, adopted for these calculations. Higher or lower XUV fluxes would lead to higher or lower escape rates, respectively.

thumbnail Fig. 4

Ratio between the escape rates calculated with the hydrodynamic model Lhy and with the energy limited formulation Len as a function of the restricted Jeans escape parameter Λ. The dashed horizontal line indicates the equality between Lhy and Len. The symbols are colour-coded as in the top panel and Fig. 1.

Open with DEXTER

4.2 Impact of the uncertainties of the stellar parameters

As mentioned in Sect. 2.1, there is a significant difference in the stellar mass estimated by David et al. (2016) and Mann et al. (2016). This study is based on the stellar mass of 0.56 M derived by Mann et al. (2016) and we discuss here the impact of this choice on the results. We re-calculated the grid of planetary upper atmosphere models considering a stellar mass of 0.3 M (David et al. 2016). We found that, because of the decreased gravitational influence of the star on the planet, hence of the larger planetary Roche lobe, the escape rates for the lower mass planets are about half of those obtained considering the higher mass star. For the higher mass planets, instead, the modification of the gravitational potential has no effect on the escape rates.

The lower stellar mass leads also to lower stellar XUV fluxes. This is due to the dependence of the X-ray luminosity on the stellar mass, such that LX,0.56 ≈ 1.8 LX, 0.3, where the subscript indicates the stellar mass. Thus, when considering the lower stellar mass, the escape rates decrease by about 36% compared to those derived considering the higher stellar mass.

A further, closely related, uncertainty is that of the stellar bolometric luminosity, which enters also in the calculation of the XUV fluxes. By considering for example the estimate of Lbol based on the stellar age (Preibisch & Feigelson 2005), the XUV fluxes would increase by a factor of two. This would in turn lead to an increase of 40% in the escape rates for the most massive planets. These considerations indicate that the major uncertainties in the stellar parameters lead to escape rates that are within a factor of two of those presented in Sect. 3.

4.3 Possible impact of the modelling formalism

We discuss here how the results depend on the major modelling assumptions. One dimensional models are widely used for exoplanet studies, particularly with respect to upper atmosphere modelling (see e.g. Murray-Clay et al. 2009; Shaikhislamov et al. 2014; Salz et al. 2016; Erkaev et al. 2016, 2017; Fossati et al. 2017; Lammer et al. 2016). An upgrade to 2D or 3D simulations is instead required to model the interaction between the expanding planetary atmosphere with the stellar wind, particularly when considering ion pick-up processes related to planetary atmospheric escape. The mass loss expected to occur because of charge exchange processes is, however, about an order of magnitude smaller than the escape caused by XUV heating (see e.g. Kislyakova et al. 2014), and hence does not significantly affect our results. Furthermore, the stellar wind does not play a significant role if the planetary outflow is supersonic before the stellar-planetary interaction occurs, because the stellar wind particles do not penetrate into the deeper atmospheric layers. This is the case of K2-33b, except possibly for the lowest considered masses. However, for these cases, the planetary atmosphere fully escapes within 5 Myr.

The presence of strong transversal planetary winds would redistribute the thermal energy in the lower atmosphere from the day to the night side, thus lowering the temperature at the lower boundary of our simulations. Given the small dependence of our results on the planetary equilibrium temperature, we expect that a lower temperature by a few hundred degrees would not significantly affect our results.

As was stated in Sect. 2.2, we assume solar metallicity when estimating the lower pressure boundary P0 . We can therefore consider the effect of a different metallicity varying the lower pressure boundary. We thus conducted additional runs assuming a metallicity of 10× and 0.1× solar. We found that, at a temperature of 850 K, the escape rates change by about 5% in the case of the Mpl = 20 M planet and by about 15% in the case of 7 M planet. Other, more local, parameters, such as temperature profiles, ionization rate, and position of the dissociation barrier, vary within 1%.

Finally, we tested the effect of the choice of chemical reaction rates. We ran further tests on the two simulated planets mentioned above (Teq = 850 K;Mpl = 7 and 20 M) substituting the reaction rates given by Yelle (2004; used in this work) with those provided by the UMIST Database for Astrochemistry (UDfA)2 . We found that the escape rates vary by less than 0.5%, with the largest variations (still below 5%) being for the ionization rate.

4.4 Evolutionary status

We discuss here whether the derived escape rates, together with the age of the system, can provide any constraint on the planet’s minimum mass. We do not aim to realistically trace the evolution of the planet’s atmosphere, which would require more sophisticated models. We aim instead at providing an order of magnitude estimation for the timescale expected for the complete escape of a hydrogen-dominated atmosphere. By comparing it with the system’s age, we derive an estimate of the possible minimum planetary mass. To this end, one needs to assume a value for the current envelope mass fraction f, which is unknown. In the literature, the maximum f value that is commonly considered is 50% (e.g. Rogers et al. 2011; Mordasini et al. 2012; Owen & Wu 2016). We therefore consider possible envelope mass fractions of 50%, 30%, 10%, and 1%. Given also the weak dependence of the escape rates on the equilibrium temperature, we consider the results obtained for a Teq value of 850 K.

Figure 5 shows the timescale estimated for the complete escape of a hydrogen-dominated atmosphere (i.e. ratio between atmospheric mass and mass-loss rate) as a function of planetary mass, for the four considered f values. This plot shows that, if the planet had amass smaller than 7 M, the planetary hydrogen-dominated envelope would have completely escaped by now. However, such a low planetary mass would imply a low bulk density, and therefore the presence of a thick hydrogen-dominated atmosphere, which is in contradiction with the previous statement. The result is that the planet cannot be less massive than 7 M, regardless of the considered envelope mass fraction. In the same way, if the planet has an f value smaller than 30% (or 10%), one could also safely exclude planetary masses smaller than 10 M (or 15 M). For masses larger than 15–20 M, the planet would keep its hydrogen-dominated atmosphere, even considering small envelope mass fractions, hence ending up as a hot Neptune.

These estimates do not consider that the stellar XUV flux, hence the atmospheric escape rates, will decrease with time. This will lead to a slight decrease of the possible minimum mass. To conclude, we can safely exclude the possibility that K2-33b is less massive than 7 M, or 10 M if the envelope mass fraction is below 30%.

thumbnail Fig. 5

Timescale, in Myr, estimated for the complete escape of a hydrogen-dominated atmosphere as a function of planetary mass, assuming four different values of the atmospheric mass fraction f of 50% (black dot), 30% (red asterisk), 10% (blue diamond), and 1% (green triangle). The horizontal dashed line indicates the maximum age (i.e. 20 Myr) derived for the K2-33 system (Mann et al. 2016).

Open with DEXTER

4.5 Temperature minima

As mentioned in Sect. 3, the structure of the temperature profiles shown in Fig. 2 for planets with masses between 7 and 15 M is particularly unusual. To understand the origin of the temperature minima, we look at the effects of H3+ cooling, X-ray heating, and EUV heating.

4.5.1 H cooling

Interestingly, most of the cooling caused by H3+ molecules is located close to the position at which the temperature minimum occurs. We therefore ran simulations for the 10 M planet without considering cooling, obtaining results that were not appreciably different from those of the original runs with the inclusion of H3+ cooling. This result is in agreement with those of Shaikhislamov et al. (2014) and Chadney et al. (2015) who did similar tests for close-in Jupiter-mass planets.

4.5.2 X-ray heating

The X-ray heating rate is comparable to that produced by EUV, but, because of the much smaller absorption cross-section (σX ~ 0.5 × 10−3σEUV), the X-rays penetrate significantly deeper into the atmosphere and deposit their energy in a thin layer where the atmosphere is dominated by H2 molecules (Fig. 6). Already at 1.5 Rpl, X-ray heating becomes negligible compared to that produced by the EUV. Additional models run without the inclusion of the stellar X-ray fluxes show that X-ray heating is not responsible for the temperature minima either.

thumbnail Fig. 6

Volume heating rates of EUV (solid) and X-ray (dashed) stellar radiation for the 10 M planet with an equilibrium temperature of 850 K.

Open with DEXTER

4.5.3 EUV heating

We further explored the effects of the intense EUV flux on the temperature profiles by carrying out additional simulations with the EUV fluxes reduced by a factor of 100 and 1000, compared to what we estimated for K2-33. We found that the amplitude of the minimum decreases in the first case and vanishes completely in the second.

Figure 7 shows the atmospheric profiles for temperature, velocity, and EUV volume heating rate for the planets with a mass of 10 M and an equilibrium temperature of 850 and 1300 K. For each planet, we ran two simulations: one considering the EUV flux estimated for K2-33, and one considering the same EUV flux reduced by a factor of 103 . In all cases, the atmospheres are characterized by adiabatic cooling close to Rpl, but only the planets subject to the stronger EUV irradiation present a temperature minimum close to the position of the peak of the volume heating rate.

The observed temperature minima are therefore caused by the extreme EUV heating. For the more irradiated planets, the EUV heating rate is about five orders of magnitude higher than that of the less irradiated planets, which also leads to the formation of a hydrodynamic flow much closer to Rpl . This results in a strong adiabatic cooling caused by the fast gas flow, which leads to a local decrease in the temperature. However, when the heating reaches values of about 1.5–3 × 10−6 erg cm−3 s−1, it dominates the adiabatic cooling, producing the observed temperature minimum.

thumbnail Fig. 7

Top three panels: EUV volume heating rate (top), bulk velocity (middle), and temperature (bottom) profiles for the 10 M planet with an equilibrium temperature of 850 K. The solid line is for the planet irradiated by the XUV flux estimated for K2-33 and the dash-dotted line is for the same planet, but with a 1000 times weaker level of XUV irradiation. The vertical dashed lines mark the position at which the velocity begins to rise (left) and that of the temperature minimum for the most irradiated planet (right). Bottom three panels: same as the top three panels, but for a planet with an equilibrium temperature of 1300 K. For the lower and higher temperature planets, Reff lies at about 2.3 and 3.0 Rpl, respectively.

Open with DEXTER

5 Conclusions

The young K2-33b planet offers the unique opportunity to study the very early phases of planetary evolution at a time when the planetary atmosphere is subject to extremely intense high-energy stellar fluxes. This is therefore the time when atmospheric escape is strongest and shapes the planetary envelope, setting it onto its final evolutionary path. We employed an upper atmosphere hydrodynamic model to study the physical characteristics of the planetary envelope and infer the escape rates. In order to more adequately model such a young system, we included in the model X-ray heating and H3+ cooling, in addition to the processes included in the model in previous studies (Erkaev et al. 2016). In the absence of a measured planetary mass, we consider masses ranging between 2 and 40 M, while the planetary radius was kept fixed to the measured value.

We found that for values of the restricted Jean escape parameter Λ smaller than 10, the escape rates are considerably larger than those derived from the energy-limited formula. For these planets, the escape is driven by the high intrinsic thermal energy of the atmosphere and low planetary gravity, while XUV heating plays a negligible role. For planets with a Λ value between 10 and 40, the hydrodynamic escape rates are comparable to those provided by the energy limited formula, while for larger Λ values they are lower.

For the lowest mass planets, the temperature profile is characterized by adiabatic cooling, caused by the strong hydrodynamic outflow of the atmosphere. For the highest mass planets, the temperature profiles present a maximum close to the position of maximum absorption of the stellar EUV flux (we find that X-ray heating plays a negligible role compared to EUV heating). For the intermediate mass planets, we find instead the presence of a local minimum close to the planetary radius. Additional dedicated simulations showed that this feature is caused by the combination of adiabatic cooling and the absorption of the strong stellar EUV fluxes.

Finally, we roughly estimate whether the derived escape rates, together with the age of the system, can provide any constraint on the minimum mass possible for the planet. We found that for masses smaller than 7 M, the planet would have completely lost its hydrogen-dominated envelope, which is in contradiction with the low average density that would result from a low planetary mass. We therefore place a lower limit of 7 M on the planetary mass.

Acknowledgements

We acknowledge the Austrian Forschungsförderungsgesellschaft FFG project “TAPAS4CHEOPS” P853993, the Austrian Science Fund (FWF) NFN project S11607-N16, and the FWF project P27256-N27. NVE acknowledges support by the RFBR grant No. 15-05-00879-a and 16-52-14006 ANF_a.

References

  1. Banks, P. M., & Kockarts, G. 1973, Aerospace (New York: Academic Press) [Google Scholar]
  2. Baraffe, I., Alibert, Y., Chabrier, G., & Benz, W. 2006, A&A, 450, 1221 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Black, J. H. 1981, MNRAS, 197, 553 [NASA ADS] [Google Scholar]
  4. Cecchi-Pestellini, C., Ciaravella, A., Micela, G., & Penz, T. 2009, A&A, 496, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Chadney, J. M., Galand, M., Unruh, Y. C., Koskinen, T. T., & Sanz-Forcada, J. 2015, Icarus, 250, 357 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180 [NASA ADS] [CrossRef] [Google Scholar]
  7. Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cubillos, P. E. 2016, ArXiv e-prints [arXiv:1604.01320] [Google Scholar]
  9. Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017, MNRAS, 466, 1868 [NASA ADS] [CrossRef] [Google Scholar]
  10. David, T. J., Hillenbrand, L. A., Petigura, E. A., et al. 2016, Nature, 534, 658 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  11. de Wit, J., & Seager, S. 2013, Science, 342, 1473 [NASA ADS] [CrossRef] [Google Scholar]
  12. Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Erkaev, N. V., Lammer, H., Odert, P., et al. 2015, MNRAS, 448, 1916 [NASA ADS] [CrossRef] [Google Scholar]
  14. Erkaev, N. V., Lammer, H., Odert, P., et al. 2016, MNRAS, 460, 1300 [NASA ADS] [CrossRef] [Google Scholar]
  15. Erkaev, N. V., Odert, P., Lammer, H., et al. 2017, MNRAS, 470, 4330 [NASA ADS] [CrossRef] [Google Scholar]
  16. Fossati, L., Erkaev, N. V., Lammer, H., et al. 2017, A&A, 598, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Ida, S., & Lin,D. N. C. 2008, ApJ, 673, 487 [NASA ADS] [CrossRef] [Google Scholar]
  18. Jin, S., & Mordasini, C. 2017, ApJ, 853, 163 [NASA ADS] [CrossRef] [Google Scholar]
  19. Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [NASA ADS] [CrossRef] [Google Scholar]
  20. Johnstone, C. P., Güdel, M., Lüftinger, T., Toth, G., & Brott, I. 2015, A&A, 577, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Kislyakova, K. G., Johnstone, C. P., Odert, P., et al. 2014, A&A, 562, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Kislyakova, K. G., Pilat-Lohinger, E., Funk, B., et al. 2016, MNRAS, 461, 988 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121 [NASA ADS] [CrossRef] [Google Scholar]
  24. Lammer, H., Erkaev, N. V., Fossati, L., et al. 2016, MNRAS, 461, L62 [NASA ADS] [CrossRef] [Google Scholar]
  25. Lecavelier des Etangs, A., Vidal-Madjar, A., McConnell, J. C., & Hébrard, G. 2004, A&A, 418, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606 [NASA ADS] [CrossRef] [Google Scholar]
  27. Locci, D., Cecchi-Pestellini, C., Micela, G., Ciaravella, A., & Aresu, G. 2018, MNRAS, 473, 447 [NASA ADS] [CrossRef] [Google Scholar]
  28. Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lopez, E. D., & Fortney, J. J. 2016, ApJ, 818, 4 [NASA ADS] [CrossRef] [Google Scholar]
  30. Lorenzani, A., & Palla, F. 2001, in From Darkness to Light: Origin and Evolution of Young Stellar Clusters, ASP Conf. Proc., 243, 745 [NASA ADS] [Google Scholar]
  31. Lubow, S. H., & Ida, S. 2010, in Exoplanets, ed. S. Seager (Tucson, AZ: Univ. of Arizona Press), 347 [Google Scholar]
  32. Mann, A. W., Newton, E. R., Rizzuto, A. C., et al. 2016, AJ, 152, 61 [NASA ADS] [CrossRef] [Google Scholar]
  33. Miller, S., Stallard, T., Tennyson, J., & Melin, H. 2013, J. Phys. Chem. A, 117, 9770 [CrossRef] [Google Scholar]
  34. Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [NASA ADS] [CrossRef] [Google Scholar]
  36. Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 578, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [NASA ADS] [CrossRef] [Google Scholar]
  38. Owen, J. E., & Morton, T. D. 2016, ApJ, 819, L10 [NASA ADS] [CrossRef] [Google Scholar]
  39. Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [NASA ADS] [CrossRef] [Google Scholar]
  40. Owen, J. E., & Wu, Y. 2016, ApJ, 817, 107 [NASA ADS] [CrossRef] [Google Scholar]
  41. Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [NASA ADS] [CrossRef] [Google Scholar]
  42. Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Preibisch, T., & Feigelson, E. D. 2005, ApJS, 160, 390 [NASA ADS] [CrossRef] [Google Scholar]
  44. Rafikov, R. R. 2006, ApJ, 648, 666 [NASA ADS] [CrossRef] [Google Scholar]
  45. Raymond, S. N., & Cossou, C. 2014, MNRAS, 440, L11 [NASA ADS] [CrossRef] [Google Scholar]
  46. Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
  47. Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2016, A&A, 585, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Schlaufman, K. C., Lin, D. N. C., & Ida, S. 2009, ApJ, 691, 1322 [NASA ADS] [CrossRef] [Google Scholar]
  50. Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
  51. Shaikhislamov, I. F., Khodachenko, M. L., Sasunov Yu. L., et al. 2014, ApJ, 795, 132 [NASA ADS] [CrossRef] [Google Scholar]
  52. Shematovich, V. I., Ionov, D. E., & Lammer, H. 2014, A&A, 571, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Stökl, A., Dorfi, E., & Lammer, H. 2015, A&A, 576, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Stökl, A., Dorfi, E. A., Johnstone, C. P., & Lammer, H. 2016, ApJ, 825, 86 [NASA ADS] [CrossRef] [Google Scholar]
  55. Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41 [NASA ADS] [CrossRef] [Google Scholar]
  56. Vanderburg, A., Latham, D. W., Buchhave, L. A., et al. 2016, ApJS, 222, 14 [NASA ADS] [CrossRef] [Google Scholar]
  57. Vidal-Madjar, A., Désert, J.-M., Lecavelier des Etangs, A., et al. 2004, ApJ, 604, L69 [NASA ADS] [CrossRef] [Google Scholar]
  58. Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150 [NASA ADS] [CrossRef] [Google Scholar]
  59. Wright, J. T. 2004, AJ, 128, 1273 [NASA ADS] [CrossRef] [Google Scholar]
  60. Wright, N. J., Drake, J. J., & Civano, F. 2011, in 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ASP Conf. Ser., 448, 1333 [NASA ADS] [Google Scholar]
  61. Yelle, R. V. 2004, Icarus, 170, 167 [NASA ADS] [CrossRef] [Google Scholar]

1

We note that here r and v are called R and V in Erkaev et al. (2016).

All Tables

Table 1

Dependence of the lower pressure boundary (in mbar) as a function of planetary Teq and mass (Mpl).

Table 2

Reactions and relative cross-sections employed in the model.

Table 3

Physical properties of the modelled planets and escape and mass-loss rates derived from the simulations.

All Figures

thumbnail Fig. 1

Hydrodynamic escape rates Lhy as a function of planetary mass (top) and restricted Jeans escape parameter Λ (bottom). Asterisks and circles mark models for which Lhy is higher and lower, respectively, than that estimated from the energy-limited formulation (see Sect. 4.1, for details). The dashed line in the bottom panel denotes Eq. (22). The vertical dotted line indicates the minimum Λ value for which Eq. (22) is valid.

Open with DEXTER
In the text
thumbnail Fig. 2

Temperature profiles of the modelled planets as a function of planetary mass and equilibrium temperature. For planets with a mass between 7 and 15 M, the temperature profile presents a minimum close to Rpl.

Open with DEXTER
In the text
thumbnail Fig. 3

Gas velocity profiles for the planets with masses of 4 (top-left), 7 (top-right), 20 (bottom-left), and 40 (bottom-right) M. The velocity profiles for planets with 10 and 15 M are very similar to those of the 20 M planet. The colours correspond to different values of the equilibrium temperature (blue: 850 K; red: 1000 K; yellow: 1150 K; magenta: 1300 K). The line style indicates whether a certain velocity is below (solid) or above (dash-dotted) the local sound speed.

Open with DEXTER
In the text
thumbnail Fig. 4

Ratio between the escape rates calculated with the hydrodynamic model Lhy and with the energy limited formulation Len as a function of the restricted Jeans escape parameter Λ. The dashed horizontal line indicates the equality between Lhy and Len. The symbols are colour-coded as in the top panel and Fig. 1.

Open with DEXTER
In the text
thumbnail Fig. 5

Timescale, in Myr, estimated for the complete escape of a hydrogen-dominated atmosphere as a function of planetary mass, assuming four different values of the atmospheric mass fraction f of 50% (black dot), 30% (red asterisk), 10% (blue diamond), and 1% (green triangle). The horizontal dashed line indicates the maximum age (i.e. 20 Myr) derived for the K2-33 system (Mann et al. 2016).

Open with DEXTER
In the text
thumbnail Fig. 6

Volume heating rates of EUV (solid) and X-ray (dashed) stellar radiation for the 10 M planet with an equilibrium temperature of 850 K.

Open with DEXTER
In the text
thumbnail Fig. 7

Top three panels: EUV volume heating rate (top), bulk velocity (middle), and temperature (bottom) profiles for the 10 M planet with an equilibrium temperature of 850 K. The solid line is for the planet irradiated by the XUV flux estimated for K2-33 and the dash-dotted line is for the same planet, but with a 1000 times weaker level of XUV irradiation. The vertical dashed lines mark the position at which the velocity begins to rise (left) and that of the temperature minimum for the most irradiated planet (right). Bottom three panels: same as the top three panels, but for a planet with an equilibrium temperature of 1300 K. For the lower and higher temperature planets, Reff lies at about 2.3 and 3.0 Rpl, respectively.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.