© ESO 2020

1 Introduction

The Kepler satellite has revealed that planets are ubiquitous in the Milky Way (Coughlin et al. 2016; Thompson et al. 2018). However, most of the planets detected by Kepler are located too far away from us to make it possible to achieve an accurate estimation of their masses and characterisation of their atmospheres. To obviate this problem, the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015) and, in the near future, the PLATO mission (Rauer et al. 2014), aim to look for transiting planets orbiting bright, nearby stars. The first planet detected by TESS, π Men c, demonstrates the success of this strategy (Gandolfi et al. 2018; Huang et al. 2018). π Men c orbits a bright (V = 5.65 mag), nearby (d ≈ 18.3 pc) G0 V star previously known to host a long-period sub-stellar companion (π Men b).

Thanks to the stellar brightness, both planetary mass and radius have been successfully measured with high precision, resulting in a mass of 4.52 ± 0.81 M_⊕ and a radius of 2.06 ± 0.03 R_⊕ (Gandolfi et al. 2018). These measurements indicate that π Men c is a super-Earth with a bulk density of about 2.8 g cm⁻³, suggesting that the planet may host a significant atmosphere, possibly water-rich, which would become hydrogen dominated in the upper layers following the dissociation of water (García Muñoz et al. 2020). Preliminary one-dimensional (1D) hydrodynamic simulations (Kubyshkina et al. 2018a) of a supposedly hydrogen-dominated atmosphere showed that such an atmosphere may be subject to strong escape of the order of 1.2 × 10¹⁰ g s⁻¹, corresponding to roughly 1% of the planetary mass per Gyr (Gandolfi et al. 2018). Such a strong escape, together with a hydrogen-dominated atmosphere, would imply the presence of an extended gaseous envelope that could be directly probed by Lyα transit observations (e.g. Vidal-Madjar et al. 2003).

García Muñoz et al. (2020) presented the results of one Lyα transit observation carried out with the STIS spectrograph on board HST. The data led to a clear detection of the stellar Lyα line, but not of any planetary absorption. They obtained 1σ upper limits for the planet-to-star radius ratio at the wavelengths covered by the Lyα line (R_p,Lyα/R_star) of 0.13 and 0.12 in the [−215,−91] km s⁻¹ and [+57,+180] km s⁻¹ velocity ranges, respectively. Furthermore, the detection and reconstruction of the stellar Lyα line enabled them to constrain the stellar high-energy (X-ray + EUV; 5–912 Å; hereafter XUV) flux to 1350 erg cm⁻² s⁻¹ at the position of π Men c (i.e. about 6 erg cm⁻² s⁻¹ at 1 AU, which is close to the solar value). This estimate also takes into account the results of France et al. (2018) obtained from an HST/COS far-ultraviolet spectrum and those of King et al. (2019) obtained from archival X-ray observations. In particular, France et al. (2018) and King et al. (2019) derived stellar XUV fluxes at the distance of the planet of 1060 erg cm⁻² s⁻¹ and 1810 erg cm⁻² s⁻¹, which are within a factor of ≈1.3 from that given by García Muñoz et al. (2020).

Following their presentation of the HST observations, García Muñoz et al. (2020) showed the results of 1D hydrodynamic simulations, accounting for hydrogen and oxygen (photo)chemistry, of the planetary upper atmosphere that were aimed at reproducing the lack of planetary Lyα absorption. They showed that the mass-loss rate is somewhat sensitive to the atmospheric bulk composition, but that a water-dominated atmosphere could comply with the Lyα non-detection because the presence of a large amount of oxygen would reduce the extension of the atmosphere, which would also turn from mostly neutral to mostly ionised at a low enough altitude to make neutral hydrogen not detectable at Lyα during transit. Vidotto & Cleary (2020) presented 1D hydrodynamic simulations of a stellar-planetary wind interaction proposing that the non-detection of hydrogen absorption during transit may be due to the confinement below the sonic point of the planetary atmosphere by the stellar wind.

Here, we present fully three-dimensional (3D) hydrodynamic modelling of the atmosphere of π Men c, composed of hydrogen and helium, and its interaction with the stellar wind. Although the model does not take into account the range of (photo)chemistry that would be needed to simulate an atmosphere containing possibly large amounts of elements heavier than helium (García Muñoz et al. 2020), it overcomes the limitations of the 1D assumption imposed by García Muñoz et al. (2020) and Vidotto & Cleary (2020). The simulations cover two values of the stellar wind velocity and a wide range of stellar XUV emission and mass-loss rate values. With this multi-parameter study, we aim at identifying the physical conditions that would reproduce the non-detection of the planetary Lyα absorption. Section 2 gives a brief presentation of the modelling framework, while Sect. 3 presents the results. In Sect. 4, we discuss the results and draw our conclusions.

2 The theoretical framework

We employ the 3D multi-fluid hydrodynamic model of Shaikhislamov et al. (2018, 2020) and Khodachenko et al. (2019), which is an upgraded version of the two-dimensional model presented by Khodachenko et al. (2015, 2017) and Shaikhislamov et al. (2016). Here, we present the most relevant features of the model. The code numerically solves the hydrodynamic continuity, momentum, and energy equations for all species in the simulated multi-component flow. In this work, we consider the atmosphere to be composed by hydrogen and helium, with a helium mixing ratio of He/H = 0.1 (approximately the solar He/H abundance ratio), and account for H, H⁺, H₂, H, H, He, and He⁺. Including H is significant because it influences the planetary mass loss by up to 30% (see also García Muñoz et al. 2020). The energetic neutral atoms (ENAs) generated by charge exchange between H (of the planetary wind; hereafter PW) and H⁺ (of the stellar wind; hereafter SW) are calculated as an independent fluid, because their velocity and temperature are significantly different from those of the neutral hydrogen atoms of planetary origin. The photochemical reactions of the H–He plasma are described in Khodachenko et al. (2015) and Shaikhislamov et al. (2016).

Photoionisation of both hydrogen and helium results in strong heating by the produced photoelectrons, which is the driver of the hydrodynamic outflow of the planetary atmosphere. The model derives the corresponding heating term by integrating the stellar XUV spectrum of García Muñoz et al. (2020). As shown, for example, by Khodachenko et al. (2019), the heating term is computedas

where N_tot is the total density of all particles, including electrons, γ_a is the adiabatic specific heat ratio that we take being 5/3, n_a is the density of hydrogen atoms, hν is the photon energy, E_ion is the hydrogen ionisation energy of 13.6 eV, σ_XUV is the wavelength-dependent ionisation cross section to XUV radiation, F_XUV is the XUV stellar flux at the planetary orbital distance, n_e is the electron density, ν_Te is the terminal velocity of electrons, E₂₁ is the n = 1 level to n = 2 level hydrogen excitation energy, and σ₂₁ and σ_ion are the hydrogen excitation and ionisation cross sections by electron impact, respectively. Under these conditions, the photoionisation time of unshielded hydrogen atoms at the planetary orbit is about four hours.

The model equations are solved in a non-inertial spherical reference frame fixed at the planetary center and rotating at the same rate as the planet orbits around the star, so that the planet faces the star always with the same side. The X axis connects the star and the planet, while the Y axis is perpendicular to the X axis and lies in the ecliptic plane. The polar-axis Z is directed perpendicular to the ecliptic plane and completes the so-called tidally locked spherical reference frame. In this frame, we properly account for the non-inertial terms, namely the generalised gravity potential and Coriolis force. Since the gas beyond the planetary exobase is mostly ionised, the species in the simulation domain are collisional through the action of the Coulomb force, justifying the hydrodynamic approach (e.g. Debrecht et al. 2020; Vidotto & Cleary 2020). Furthermore, the model also takes into account radiation pressure, which we find to be mostly negligible compared to other forces (e.g. Murray-Clay et al. 2009; Khodachenko et al. 2017, 2019; Debrecht et al. 2020). We tested this by computing models with twice and ten times larger stellar Lyα emission compared to the baseline of 5.6 erg s⁻¹ cm⁻² at 1 AU. We found that only the model with the highest Lyα flux affects the Lyα planetary absorption profile, but only at velocities between −100 and 0 km s⁻¹, which are still too low to be detectable by the observations and would be contaminated by interstellar medium (ISM) absorption and geocoronal emission (García Muñoz et al. 2020).

The fluid velocity at the planetary surface is assumed to be zero. Furthermore, to keep the number of grid points in the model small enough to be manageable by the numerical code, the radial mesh is highly non-uniform, with the grid step increasing linearly from the planetary surface. This allows us to resolve the highly stratified upper atmosphere of the planet, where the required grid step is as small as Δr = R_p/400, where r is the radial distance from the center of the planet and R_p the planetary radius. As the initial state, we took a fully neutral atmosphere in barometric equilibrium composed by H₂ and He.

For all simulations, we considered the system parameters of Gandolfi et al. (2018), which we reproduce in Table 1. At the inner boundary of the simulation domain, at r = R_p, we set a temperature and a pressure of 1000 K and 0.05 bar, respectively. The former is close to the planetary equilibrium temperature obtained considering zero albedo (see Table 1), while the latter is the pressure at which the planetary atmosphere is optically thick to photons with a wavelength longer than of 10 Å (i.e., the blue edge of the XUV range).

The model also self-consistently incorporates the flow of the SW plasma, which we consider being composed by protons, in the way described by Khodachenko et al. (2019) and Shaikhislamov et al. (2020), enabling one to model the whole system. Thus, the simulation domain has one further boundary at the base of the stellar corona. At distances from the stellar center shorter than 20 stellar radii, the SW is accelerated by an empirical heating term, derived from an analytical 1D polytropic Parker-like model (Keppens & Goedbloed 1999), that we compute as (3)

In Eq. (3), γ_a is the adiabatic specific heat ratio that we take being 5/3 to correctly model shocks and r_s is the distance from the center of thestar, while γ_p, T_p, and V_p are respectively the polytropic index, temperature, and velocity obtained from the polytropic solution. The simulated stellar wind is isotropic in space, stationary in time, and in good agreement with the Parker analytical solution.

We further employ the simulations to compute the absorption at the position of the Lyα line following Shaikhislamov et al. (2018). This procedure has already been successfully employed to generate synthetic observations for the hot Jupiter HD 209458 b (Khodachenko et al. 2017; Shaikhislamov et al. 2018, 2020) and the warm Neptune GJ 346 b (Khodachenko et al. 2019), providing physically reasonable and self-consistent interpretation of Lyα transit observations, as well as of transit observations of resonance lines of minor species (e.g. C, O, Si).

3 Results

The velocity of the PW driven by XUV heating reaches moderate values of about 10 km s⁻¹, which, by far, are not enough to produce Doppler-shifted absorption beyond the wavelengths strongly contaminated by ISM absorption and geocoronal emission. Therefore, absorption in the wings of the Lyα line during transit can be caused only by ENAs generated during the interaction of the PW with the SW.

π Men is an ≈5 Gyr old solar-type star (Gandolfi et al. 2018), which is therefore likely to have solar-like wind properties. To explore how strong the Lyα absorption might be under different physical conditions, we simulate the PW-SW interaction considering two winds typical of solar-like plasma, namely a fast and a slow wind of terminal velocities V_sw,∞ = 800 km s⁻¹ and 400 km s⁻¹, respectively, further varying the stellar mass-loss rate, hence the SW density, by about an order of magnitude. Another parameter directly affecting the Lyα absorption is the stellar XUV emission, which we also vary by about an order of magnitude.

Figure 1 presents simulation results for two cases. The first case is for a very weak SW with a density ten times smaller than solar. The stellar mass-loss rate in this case is of 10¹¹ g s⁻¹. Under these conditions, the SW, which is still subsonic at the planetary orbital separation, diverts the planetary outflow only far from the planet, at an X = r/R_p value of ≈200–300. Therefore, the SW-PW interaction occurs in a low-density region of the planetary atmosphere, leading to a small ENA production. The second case is for a moderate SW, namely with a stellar mass-loss rate of 2 × 10¹² g s⁻¹ (the average solar mass-loss rate is about 2.5 × 10¹² g s⁻¹; e.g. Phillips et al. 1995). This second case is significantly different from the previous one, because the SW pressure is enough to stop the planetary outflow close to the planet, at X ≈20–30, redirecting the planetary escaping material towards the tail, further generating a large amount of ENAs. Because of its supersonic nature, this interaction generates a bow shock, as evinced by the temperature distribution (Fig. 1). The ENAs, which produce significantly high velocity absorption in the blue wing of the Lyα line, are generated inside the bow shock region.

Figure 2 presents the details of the simulated distribution of species computed for the moderate SW (i.e. top-right and bottom panels of Fig. 1) along the star-planet line. It shows that the simulation domain is split by the shock into two regions, one dominated by the escaping planetary material and one dominated by SW material. The supersonic SW remains undisturbed until the bow shock region, where it experiences a sharp deceleration, compression, and heating. The region dominated by planetary material, instead, is characterised by a relatively low temperature and velocity, though before the ionopause the PW is supersonic. At the ionopause, the normal component of the proton velocity goes to zero. However, the neutral PW particles penetrate into the shock region, where they are rapidly ionised by the hot electrons and XUV radiation. In the shocked region, adjacent to the ionopause, the energetic SW protons charge exchange with planetary atoms producing ENAs. For this simulation, we obtained a planetary mass-loss rate of about 2 × 10¹⁰ g s⁻¹, which is in good agreement with previous estimates based on 1D simulations (Gandolfi et al. 2018; García Muñoz et al. 2020); we obtain roughly the same planetary atmospheric mass-loss rate (within a few %) for all conducted simulations. We remark that, besides the difference in the geometry of the simulations, the mass-loss rates of Gandolfi et al. (2018) were obtained accounting only for H and those of García Muñoz et al. (2020) for H and O, while our simulations consider H and He.

The result that models with different assumptions and geometries lead to similar planetary mass-loss rates is found also for other planets (see, for example, Table 2 of Kubyshkina et al. 2018b). This may be because the physicalproperties (mainly density and velocity) of the gas at the sonic point are robust against the typical assumptionstaken in the different codes, though detailed comparisons of the model outputs would be necessary to confirm this. However, the most interesting result is possibly the similarity of mass-loss rates computed by 1D and 3D models, where the former are integrated across the whole planet (i.e., 1D mass-loss rates multiplied by 4π). Indeed, although it might be expected that mass loss would be strongly reduced on the night side, the large size of the upper atmosphere of close-in planets and the redistribution of heat across it lead to similar conditions throughout the majority of the upper atmosphere.

Figure 2 compiles a rough estimate of the absorption at the position of the Lyα line. The integral of the density of ENAs along the X axis (i.e. the ENA column density; hereafter NL) is equal to NL = 4.5 × 10⁴ R_p = 0.6 × 10¹⁴ cm⁻², while the resonant cross-section absorption of the Lyα line is σ = 6 × 10⁻¹⁴ cm² (Khodachenko et al. 2017). Assuming that ENAs form a shell of radius R_abs,ENAs ~ 30 R_p ~ 0.5 R_star around the planet (see Fig. 1), the Lyα absorption caused by ENAs is 0.5² NL σ ≈ 0.9, meaning that the size of the planet at Lyα wavelengths R_p,Lyα/R_star is ≈1.0.

Figure 3 shows the actual absorption profile obtained from the simulation computed considering the weak and moderate SW (i.e. top panels of Fig. 1), and the distribution of the wavelength-integrated absorption depths (i.e. 1 − exp^{−NL σ}) across the stellar disk obtained from the two simulations. The absorption profiles have been computed accounting for all neutral hydrogen particles, further considering their own velocities and temperatures. In the weak SW case, the absorption takes place mostly in the [−30,50] km s⁻¹ velocity range, where any planetary absorption signature is unobservable, because of contamination by ISM absorption and geocoronal airglow emission. In contrast, in the moderate SW case, the absorption is largest at higher negative velocities and comes mostly from the shocked region, as shown by the bottom panel of Fig. 3, which presents the distribution of the high-velocity absorption, hence from ENAs, across the stellar disk.

Finally, we performed a systematic study of the Lyα planetary absorption depth as a function of the stellar XUV flux, ranging between 3 and 20 erg s⁻¹ cm⁻² at 1 AU, and of the stellar mass-loss rate, ranging between 2 × 10¹¹ and 4 × 10¹² g s⁻¹. For this grid of models, the SW terminal velocity was kept constant at the slow SW, namely 400 km s⁻¹. We ran an additional, smaller grid considering the fast SW, namely 800 km s⁻¹, and stellar mass-loss rates, ranging between 10¹¹ and 7 × 10¹¹ g s⁻¹. We note that the estimated XUV flux of π Men at 1 AU is 6 erg s⁻¹ cm⁻², while the average solar mass-loss rate is 2.5 × 10¹² g s⁻¹, hence comparable to the strongest SW we considered (i.e., that with the highest stellar mass-loss rate). Figure 4 summarises the results of the systematic analysis. It presents the absorption in terms of R_p,Lyα/R_star integrated over the blue wing of the Lyα line in the [−215,−91] km s⁻¹ velocity range, for which García Muñoz et al. (2020) obtained 1σ and 3σ upper limits of R_p,Lyα/R_star = 0.13 and 0.24, respectively. We do not consider in the analysis the red wing of the Lyα line, because at those velocities none of the computed models led to absorption redwards of the region contaminated by ISM absorption and geocoronal emission, in agreement with the observations.

We first analyse the results obtained from the larger grid computed considering the slow SW. Figure 4 indicates that all simulations in the upper-left quadrant (i.e. lower XUV flux and larger stellar mass-loss rate) show strong Lyα absorption, in contrast to the observations. For each considered XUV flux value, there is a narrow band of stellar wind densities, following roughly the lower green full line in Fig. 4, across which the planetary absorption decreases significantly. This is because a stronger SW (i.e. larger stellar mass-loss rate) will stop the expanding PW closer to the planet (see e.g. Fig. 1), where the ENA production is efficient. In contrast, a weaker SW (i.e. smaller stellar mass-loss rate) will result in inefficient ENA production and therefore the planetary absorption in the blue wing of Lyα comes only from natural line broadening, which amounts to R_p,Lyα/R_star ≈ 0.04. Despite this small value, the planetary plasmasphere is as large as ~10 R_p because the SW is not strong enough to confine the planetary atmosphere.

Figure 4 shows that in comparison to the slow SW case (i.e. SW terminal velocity of 400 km s⁻¹), the fast SW (i.e. SW terminal velocity of 800 km s⁻¹), even with the same stellar mass-loss rate, produces significantly higher Lyα absorption. This is because the pressure of the fast SW is higher, confining the planetary atmosphere inside a bow shock closer to the planet. This implies that the fast SW interacts with a denser PW, leading to a larger ENA production that generates stronger absorption in the blue wing of the Lyα line (see bottom panel of Fig. 3).

Fig. 1 Top left: proton density distribution in the ecliptic plane across the whole simulated domain computed with a stellar XUV flux at 1 AU of 6 erg cm−2 s−1, a stellar mass-loss rate of 1011 g s−1 (i.e. weak SW), and SW terminal velocity of 400 km s−1 (i.e. velocity of 250 km s−1 and temperature of 0.65 MK at the planetary orbit), leading to a SW density at the planetary orbit of 3 × 102 cm−3. The axes arescaled in units of planetary radii and the planet sits at the center of the coordinate reference frame, while the star lies at X = 762. Top right: same as top-left, but for a stellar mass-loss rate of 2 × 1012 g s−1 (i.e. moderate SW), leading to a SW density at the planetary orbit of 6 × 103 cm−3. Bottom left: same as top-right, but for the density distribution of ENAs and closing-up to the position of the planet. Bottom right: same as bottom-left, but for the electron temperature in the X− Z plane. In all panels, the lines with arrows indicate streamlines of the corresponding fluids.

Fig. 2 Distribution profiles of major species along the planet-star line obtained in the simulation with the moderate SW (i.e. top-right and bottom panels of Fig. 1). The left axis is for the density (in log[cm−3]) of H2 (orange line), H (red line), H+ (blue line), and ENAs (green line). The right axis is for the proton temperature (in 104 K; black solid line) and velocity (in km s−1; dotted line). The vertical dashed lines indicate the approximate positions of the ionopause and of the bow shock. The planet lies at X = 0, while the star is located to the right.

Fig. 3 Top: out-of-transit (black circles) and in-transit (red circles) observed Lyα profiles (from García Muñoz et al. 2020). The dashed line shows the synthetic in-transit Lyα profile obtained from the simulation computed under the weak SW conditions (i.e. top-left panel in Fig. 1), while the solid line is for the simulation computed with the moderate SW (i.e. top-right and bottom panels of Fig. 1). Of the two simulations represented in this figure, only the weak SW case is consistent with the HST non-detection. Bottom: distribution of the Lyα absorption along the line of sight averaged over the blue wing of the line, in the [−215,−91] km s−1 velocity range, as seen by a remote observer at mid-transit and considering the moderate SW. The absorption ranges between 0 and 1, where 0 means no absorption, while 1 means full absorption. The black circle at the outer boundary indicates the star and the white horizontal line shows the planetary orbital path accounting for the impact parameter (as the planet moves from left to right). At mid-transit, the planet is located at the intersection of the horizontal white solid line and of the vertical dashed black line.

4 Discussion and conclusions

We ran 3D hydrodynamic simulations of the interaction between the expanding upper atmosphere of π Men c and of the wind of its host star, considering a planetary atmosphere composed of H and He. We ran simulations assuming two distinct values of the SW terminal velocity (400 and 800km s⁻¹) and a range of stellar XUV fluxes and mass-loss rates. We find, assuming a slow SW and for stellar XUV fluxes close to those estimated for π Men, that the non-detection of Lyα absorption during transit can be reproduced by considering stellar winds that are significantly weaker thanthe average solar wind, with a density more than six times smaller. Reproducing the Lyα non-detection employing a faster SW would require an even lower SW density. We find that with a solar-like SW, fitting the Lyα non-detection would require an improbably higher stellar XUV flux, namely about 4 times of that estimated for π Men, while the highest XUV estimate is of just ≈1.3 times larger(King et al. 2019) and uncertainties on the reconstructed XUV fluxes based on Lyα measurements are typically of the order of 30% (Linsky et al. 2014). This is because a higher XUV flux more rapidly ionises hydrogen, increasing the upper atmospheric heating and expansion, pushing the interaction region with the SW farther away from the planet. Furthermore, at such high stellar XUV fluxes the planet would have lost almost half of its current mass within the estimated age of the system, without accounting for the fact that the star had been more active in the past. On the basis of the estimated stellar XUV flux, accounting for its uncertainty and of evolutionary considerations, we conclude that a high stellar XUV emission is unlikely to be the cause of the non-detection of planetary Lyα absorption.

This result is driven by the fact that Lyα planetary atmospheric absorption at the velocity probed by the observations can be caused solely by ENAs, which become more abundant with increasing the velocity or the density of the SW. Therefore, similarly to Vidotto & Cleary (2020), we find that a stronger (i.e., faster or denser) SW compresses the planetary atmosphere on the side facing the star, reducing its size. However, a stronger SW penetrates deeper into the neutral part of the expanding planetary atmosphere, increasing the density of ENAs, thus the Lyα absorption at the velocities probed by the observations.

However, the set of simulations described in Sect. 3 do not reproduce the case assumed by Vidotto & Cleary (2020) in which the SW is so strong that it compresses the planetary atmosphere below the sonic point. Therefore, we ran a further simulation considering the slow SW (i.e. 400 km s⁻¹), a stellar XUV flux at 1 AU of 6 erg cm⁻² s⁻¹, and a stellar mass-loss rate of 1.7 × 10¹³ g s⁻¹ and density of 5 × 10⁴ cm⁻³, closely corresponding to Model A in Vidotto & Cleary (2020). In agreement with Vidotto & Cleary (2020), we find that the bow shock and ionopause lie much closer to the planet, namely at 16.5 and 10.5 R_p, respectively. We also find that the planetary flow on the day side reaches its maximum speed of 4.9 km s⁻¹ at a distance of 6 R_p, while the sound speed is 14.4 km s⁻¹, implying a subsonic flow.

However, we also find that the planetary mass loss rate is just 2.5% smaller than what we obtained from the other simulations, therefore confirming that a subsonic PW does not necessarily entail a smaller mass-loss rate (García Muñoz 2007). In one-dimensional models, this is the result of partial compensation in density and velocity changes when the prescribed downstream pressure is increased. Indeed, higher pressures slow down the flow, which responds by increasing the temperature and in turn the density (through a larger scale height). Because the mass loss rate scales with the product of density and velocity, their individual changes tend to cancel out. There is a limit to this, as for instance expected if the planetary wind temperatures become very high and the plasma loses a significant amount of energy through radiation. This picture is qualitatively consistent with the parametric study conducted by Christie et al. (2016), who report a decrease in the mass-loss rate when the planetary wind becomes weaker and easier to confine by the stellar wind (see their Fig. 7). The quantitative differences with Christie et al. (2016) probably arise from the different treatment of the equation of state of the gas. It is unclear if their prescription of a polytropic equation of state, which results in roughly isothermal temperature of the planetary wind, can capture the compensation effects described above.

Furthermore, we find that absorption in the blue wing of Lyα is still very large, namely R_p,Lyα/R_star ≈ 0.64, which is in disagreement with the observations. Vidotto & Cleary (2020) could not have reached this result, because of the assumption of 1D and, more importantly, the lack of ENAs. Therefore, while a very strong SW is indeed capable of slightly reducing the planetary mass-loss rate, it would also lead to a stronger, rather than weaker, absorption signature in the blue wing of the Lyα line profile.

Since estimates indicate that the wind strength of π Men is close to solar (or even stronger; Vidotto & Cleary 2020), our simulations clearly suggest that it is very unlikely that π Men c hosts an atmosphere dominated by hydrogen and helium, a finding that is in agreement with the considerations of Gandolfi et al. (2018) and García Muñoz et al. (2020). Therefore, we argue that, despite the rather low bulk planetary density, π Men c’s atmosphere cannot be strongly hydrogen-dominated and that it should, thus contain a non-negligible amount of heavier elements, as suggestedby García Muñoz et al. (2020). Future observations should focus on looking for elements, such as He, C, and O, to shed more light on the planetary atmospheric composition. Furthermore, additional observations are needed to check that the available measurements have not been made at non-typical and rare stellar conditions of either very low SW density or high XUV flux.

Fig. 4 Planetary Lyα absorption in Rp,Lyα/Rstar integrated inthe [−215,−91] km s−1 velocity range as a function of stellar XUV flux (at 1 AU; x-axis) and input stellar mass-loss rate (left y-axis). Red circles indicate the results obtained with a SW temperature and terminal velocity of 0.64 MK and 400 km s−1, respectively. Blue circles are for a SW temperature and terminal velocity of 1.2 MK and 800 km s−1, respectively.The size of the red circles is proportional to the planetary absorption, whose value is given beside each circle. For a comparison with the observations, the green solid lines show the approximate position of the 0.13 and 0.24 absorption levels with respect to the results given by the red circles. The gray curve gives the planetary mass-loss rate (right y-axis) obtained from the simulations as a function of stellar XUV flux in the case of the slow SW. The planetary mass-loss rate is very weakly dependent of the stellar mass-loss rate (see Sect. 4).

Acknowledgements

I.S., M.Kh., and M.R. received support by the RSF project 18-12-00080 in the frame of which the numerical modeling, key for this study, has been developed. Parallel computing has been performed at the Computation Center of Novosibirsk State University, the SB RAS Siberian Supercomputer Center, the Joint Supercomputer Center of RAS, and the Supercomputing Center of the Lomonosov Moscow State University. I.S. and M.Kh. also acknowledge the RFBR project 20-02-00520. M.Kh. acknowledges support from the projects I2939-N27 and S11606-N16 of the Austrian Science Fund (FWF). M.Kh. acknowledges grant number 075-15-2019-1875 from the government of the Russian Federation under the project called “Study of stars with exoplanets”. We thank the anonymous referee for the useful comments that led to improve the manuscript.

References

All Tables

All Figures

Fig. 1 Top left: proton density distribution in the ecliptic plane across the whole simulated domain computed with a stellar XUV flux at 1 AU of 6 erg cm−2 s−1, a stellar mass-loss rate of 1011 g s−1 (i.e. weak SW), and SW terminal velocity of 400 km s−1 (i.e. velocity of 250 km s−1 and temperature of 0.65 MK at the planetary orbit), leading to a SW density at the planetary orbit of 3 × 102 cm−3. The axes arescaled in units of planetary radii and the planet sits at the center of the coordinate reference frame, while the star lies at X = 762. Top right: same as top-left, but for a stellar mass-loss rate of 2 × 1012 g s−1 (i.e. moderate SW), leading to a SW density at the planetary orbit of 6 × 103 cm−3. Bottom left: same as top-right, but for the density distribution of ENAs and closing-up to the position of the planet. Bottom right: same as bottom-left, but for the electron temperature in the X− Z plane. In all panels, the lines with arrows indicate streamlines of the corresponding fluids.

In the text

Fig. 2 Distribution profiles of major species along the planet-star line obtained in the simulation with the moderate SW (i.e. top-right and bottom panels of Fig. 1). The left axis is for the density (in log[cm−3]) of H2 (orange line), H (red line), H+ (blue line), and ENAs (green line). The right axis is for the proton temperature (in 104 K; black solid line) and velocity (in km s−1; dotted line). The vertical dashed lines indicate the approximate positions of the ionopause and of the bow shock. The planet lies at X = 0, while the star is located to the right.

In the text

Fig. 3 Top: out-of-transit (black circles) and in-transit (red circles) observed Lyα profiles (from García Muñoz et al. 2020). The dashed line shows the synthetic in-transit Lyα profile obtained from the simulation computed under the weak SW conditions (i.e. top-left panel in Fig. 1), while the solid line is for the simulation computed with the moderate SW (i.e. top-right and bottom panels of Fig. 1). Of the two simulations represented in this figure, only the weak SW case is consistent with the HST non-detection. Bottom: distribution of the Lyα absorption along the line of sight averaged over the blue wing of the line, in the [−215,−91] km s−1 velocity range, as seen by a remote observer at mid-transit and considering the moderate SW. The absorption ranges between 0 and 1, where 0 means no absorption, while 1 means full absorption. The black circle at the outer boundary indicates the star and the white horizontal line shows the planetary orbital path accounting for the impact parameter (as the planet moves from left to right). At mid-transit, the planet is located at the intersection of the horizontal white solid line and of the vertical dashed black line.

In the text

Fig. 4 Planetary Lyα absorption in Rp,Lyα/Rstar integrated inthe [−215,−91] km s−1 velocity range as a function of stellar XUV flux (at 1 AU; x-axis) and input stellar mass-loss rate (left y-axis). Red circles indicate the results obtained with a SW temperature and terminal velocity of 0.64 MK and 400 km s−1, respectively. Blue circles are for a SW temperature and terminal velocity of 1.2 MK and 800 km s−1, respectively.The size of the red circles is proportional to the planetary absorption, whose value is given beside each circle. For a comparison with the observations, the green solid lines show the approximate position of the 0.13 and 0.24 absorption levels with respect to the results given by the red circles. The gray curve gives the planetary mass-loss rate (right y-axis) obtained from the simulations as a function of stellar XUV flux in the case of the slow SW. The planetary mass-loss rate is very weakly dependent of the stellar mass-loss rate (see Sect. 4).

In the text