© The Authors 2025

1 Introduction

Close-in exoplanets can be subject to significant atmospheric escape due to the strong irradiation they receive (Lammer et al. 2003). The ionisation of upper planetary atmospheres and the heat deposition are thought to mostly be due to the absorption of the XUV (X-ray and extreme ultraviolet radiation, from 15 to 912 Å) from the host star (Yelle 2004; García Muñoz 2007). Depending on the ionisation conditions at the location where this flux is absorbed, the input energy can ionise, heat, or excite the atmosphere (Gillet et al. 2023). When the heating is sufficiently strong, it triggers a thermal atmospheric escape that is generally partially ionised (for a review of escape mechanisms, see Gronoff et al. 2020). The escaping planetary material then interacts with the ambient, fully ionised stellar wind, leading to further ionisation and neutralisation due to, for example, chargeexchange phenomena (Holmström et al. 2008). The extended tail of neutrals, whether escaping neutrals from the planet or neutralised stellar wind (Tremblin & Chiang 2013), leaves a detectable in-transit signal of the Lyα line around hot Jupiters (Vidal-Madjar et al. 2003) and hot Neptunes (Kulow et al. 2014). In addition, strong escape can rapidly deplete the atmosphere of hot planets. It could reduce the loss time of their atmospheres to a few hundred thousand of years for low-mass planets (Kubyshkina et al. 2018), which can have a strong impact on their migration path as the stellar system evolves (Lazovik 2023). A precise, quantitative estimate of the atmospheric escape rate and of its detectable signals is therefore very much needed.

Numerous modelling efforts have been undertaken over the last two decades to quantify the amount of planetary material susceptible to escape in this manner. Starting from a 1D representation of the star-planet line, works such as Yelle (2004) and García Muñoz (2007) characterised the physical and chemical aeronomy of hot Jupiters, using the example of the hot Jupiter HD 209458 b. The modelling efforts continued, going from 1D (Chadney et al. 2017; Koskinen et al. 2022) to 3D (Rumenskikh et al. 2022), and considering various chemical compositions (García Muñoz et al. 2021; Khodachenko et al. 2021), multi-fluid approaches (Shaikhislamov et al. 2021), the effect of magnetic fields (Erkaev et al. 2017; Daley-Yates & Stevens 2019), interaction with a magnetised stellar wind (Carolan et al. 2021), and non-local thermodynamic effects (García Muñoz & Schneider 2019).

Despite these efforts, even in the simple case of an atmosphere composed purely of atomic hydrogen, the energy repartition between heating, excitation and ionisation was only recently quantified taking photoelectrons into electrons by Guo & Ben-Jaffel (2016) and Gillet et al. (2023) using 1D models of planetary atmospheres. They report strong effects on the predicted mass-loss rate for hot Jupiters (reduction in the predicted mass-loss rate by a factor of 2 when photoelectrons are taken into account), reinforcing the importance of precisely characterising the details of the energy deposition in the upper atmosphere. Furthermore, the development of realistic planetary escape models should allow the community to directly compare model results with observed atmospheric escape tracers to put planetary escape theory on firm grounds, and use available observational information to characterise the aeronomy of planets.

The repeated detection of Lyα transits in multiple systems also revealed some temporal variability. For instance, Lecavelier Des Etangs et al. (2012) showed that the Lyα transit was visible for one epoch following a strong stellar flare of HD 189733 and not visible at another epoch, questioning whether stellar variability can strongly influence the atmospheric escape. On one hand, Chadney et al. (2017) show with a 1D model of HD 189733 b that the flare timescale is too short for an individual flare to have a significant effect on the atmospheric escape (but it could strongly influence the ionosphere of the planet). On the other hand, Hazra et al. (2022) show using a 3D model that a constant XUV flux corresponding to a strong flare should strongly impact the rate of atmospheric escape. As of today, the effect of successive, repeated strong flares on the atmospheric escape remains an open question.

In this work we aim to extend the investigation carried out in Gillet et al. (2023) by considering a more realistic 2D geometry, embedding a day and a night side for the atmosphere. A 2D geometry allows us to assess the effect of photoelectrons on the detectable absorption in Lyα while keeping the computational cost reasonable with respect to a 3D geometry. In our self-consistent treatment, we model the heating, ionisation and excitation caused by XUV stellar photons in the atmosphere. We present the 2D self-consistent model that we developed to describe the thermal escape of the planetary wind and its interaction with the stellar wind in Section 2. We then compare the results obtained in 1D (Gillet et al. 2023) and those obtained in 2D in Section 3. We describe our methodology for producing synthetic Lyα transits from our simulations in Section 4, and show how photoelectrons and turbulent variability in our model affect such profiles in Section 5. We study the case of a strong flare in Section 6, and discuss and summarise our main results in Section 7.

2 Model description

2.1 Extension of the 1D physical model

Our 2D Cartesian simulations are run with the hydrodynamics code PLUTO (Mignone et al. 2007), which solves the Euler equations in a rotating reference frame. Similarly to our 1D cases in Gillet et al. (2023), the 2D equations for the conservation of mass, momentum, and energy solved in PLUTO are $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)=0,$$(1) $$\frac{\partial (\rho u)}{\partial t}+\nabla \cdot \left[\rho uu+P_{T}I\right]=-\rho \nabla \varphi +\rho \left(F_{\text{cent}}+F_{\text{cor}}\right),$$(2) $$\frac{\partial (E+\rho \varphi )}{\partial t}+\nabla \cdot \left[\left(E+P_T+\rho \varphi \right)u\right]=\rho u\left(F_{\text{cent}}+F_{\text{cor}}\right)+H-C,$$(3)

where t is the time, u is the fluid velocity, and ρ is the total density, P_T is the thermal pressure, E = P_T /(γ − 1) + ρu²/2 is the total energy with γ = 5/3, the adiabatic index for a mono-atomic gas. C and H are the cooling and heating terms defined in Gillet et al. (2023). The joint gravitational potential of the planet plus the star is defined as $$\varphi =-\frac{G{M}_p}{r}-\frac{G{M}_{\star }}{R_{\text{orbit }}-r},$$(4)

where G is the gravitational constant, and M_★ and M_p are the masses of the star and the planet, respectively. The Coriolis and centrifugal force are defined as F_cor=-2(Ω × u) and F_cent= Ω × (Ω × R) where R is the position vector relative to the star and $\text{Ω}=\sqrt{\left(G{M}_{\star }/R_{orb}^3\right)}{e}_{ɀ}$ is the orbital rotation rate with e_ɀ is the unit vector along the coordinate ɀ. The model presented here solves these equations on the (x,y) orbital plane.

We consider a planetary atmosphere composed of atomic hydrogen in neutral H I and ionised H⁺ forms plus thermal electrons. When the stellar radiation hits the planetary atmosphere, it ionises the gas and contributes to its heating. The photoionisation rate coefficient J [s⁻¹] is given by: $$J={\displaystyle {\int }_{{\lambda }_{\min }}^{{\lambda }_0}{\sigma }_{\lambda }}{F}_{\star }\left(\frac{\lambda }{hc}\right)\left(1+{\text{Φ}}_{\lambda ,xe}\right)\text{d}\lambda ,$$(5)

where F_★ [erg cm⁻²s⁻¹Å⁻¹] is the attenuated stellar flux. Φ_λ,xe is the number of secondary ions created per photoionisation. The heating deposition rate H [erg s⁻¹ cm⁻³] in the energy equation can be expressed as: $$H=n_{\text{HI}}{\displaystyle \int {\sigma }_{\lambda }}{F}_{\star }\left(1-\frac{\lambda }{{\lambda }_0}\right){\eta }_{\lambda ,x_e} \text{d}\lambda$$(6)

with ${\eta }_{\lambda ,x_e}$ being the heating efficiency or rate of conversion from kinetic energy of the photoelectron into actual heating. This equation subtracts from the photoelectron energy the fraction of energy that goes into excitation and ionisation of the gas and, therefore, does not contribute to heating. The full description of the chemistry and radiative transfer can be found in Gillet et al. (2023) along with the implementation of the photoelectrons. As there, the photoelectron physics is introduced through a parameterisation inferred from calculations in García Muñoz (2023).

2.2 Stellar wind

We considered the stellar wind to be a fluid in which all species (electrons, neutrals and ions) have the same velocity and temperature. The stellar wind is expected to have a significant impact on the structure and velocity of the planet’s atmospheric escape. When the planetary wind meets the stellar wind, its structure changes in a way similar to that of a comet. Structures of different sizes are then formed: bow shock, turbulence as well as instabilities (such as shearing due to the velocity differential between the two winds). The overall structure of the planetary wind can vary depending on the ram and thermal pressure ratios of each wind, and have been classified in particular by Matsakos et al. (2015).

If the stellar wind is supersonic, it will form a shock as it slows down to adapt to the velocity of the obstacle (the planet in this case). The velocity of the post-shock stellar wind seen by the planet will therefore be lower than the pre-shock velocity (via the Rankine-Hugoniot shock conditions). The interaction of the two winds results in an equilibrium point of thermal pressure and ram pressure (see Figure 2) with $P_{ram }=\frac{1}{2}\rho v^2$. If the equilibrium point is achieved between the ram pressure of the planetary wind and the thermal pressure of the stellar wind, then the mass-loss rate of the planetary escape remains independent of the stellar wind conditions (Christie et al. 2016). If the equilibrium point is achieved by the thermal pressures, then we have a modification of the boundary conditions of the planetary escape. Strong stellar winds can prevent the planetary wind from reaching supersonic speeds in the vicinity of the planet (Vidotto & Cleary 2020; Carolan et al. 2021). It is therefore crucial to properly model the stellar wind in order to study these interactions with the planetary wind.

In our model, the stellar wind is injected into the 2D Cartesian grid at the left edge of the domain along the y axis. We based our stellar wind prescription on the wind model of Parker (1958). We imposed a stellar wind mass-loss rate as a fixed parameter in our model and set it to 50 Ṁ_⊙ with Ṁ_⊙ = 2 × 10⁻¹⁴ M_⊙ yr⁻¹ (see also Carolan et al. 2021). Such a large mass-loss is representative of young Sun-like stars (Wood et al. 2021) and leads to a stellar-planetary wind interaction within 20 planetary radii, which appropriately limits the size of the box required for our modelling. A weaker stellar wind does not allow us to contain the planetary wind in a reasonable numerical box size. For close-in planets, a global modelling embedding also the star itself is then needed (see Matsakos et al. 2015; Strugarek et al. 2015). The stellar wind temperature is set to T_w = 1.5 × 10⁶ K. The density ρ is prescribed by the stellar mass-loss as ρ_w = M_⊙/4πr²v_w with v_w = 400 km/s and the pressure can be derived from density and temperature. We corrected v_x and v_y to take the rotating frame into account. The neutral fraction x_HI is set to zero when the stellar flow enters the box.

Fig. 1 Illustration of the dynamic characteristics of the interaction of the planetary wind (with the planet at the centre delimited by a circle) and stellar wind for the most extreme cases M0.02 (left) and M0.69 (right). Top panel: instantaneous temperature (in Kelvins) in logarithmic scale. Bottom panel: instantaneous density of neutrals in logarithmic scale.

2.3 Numerical model

We solved the Euler equations on a 2D Cartesian grid,x, y = [−20, 20]R_p parallel to the orbital plane with the planet located at centre of the box. We carried out high resolution 2D simulations (5000 cells × 5000 cells) cut into a uniform grid of 1000 × 1000 cells covering the square [−1,1] around the planet and stretched by 2000 points in each direction between ±1 and ±20 R_p . In all our cases, we initialised our model by interpolating the corresponding 1D spherical converged solution (obtained with the model described in Gillet et al. 2023) on the 2D Cartesian grid. We aimed to have the highest possible resolution between the surface of the planet and the interface of interaction with the stellar wind in order to study the mechanisms of photoionisation, energy deposition, but also resolve the small structures of the planetary wind.

For the numerical method, we used an advanced option in PLUTO named “shock flattening” in its “MULTID” version (Colella & Woodward 1984), which provides more dissipation in the star-planet interaction region where the shock occurs. To solve the Euler equations, as in 1D, we used the Riemann solver in the Harten-Lax-Van approximation,exactly resolving the contact discontinuities between the cells (Harten-Lax-van Leer, see Toro 2009). We used a linear reconstruction for the spatial order of integration and a third-order Runge-Kutta (RK3) scheme for the temporal evolution. Typically, a run takes about two weeks to reach a relaxed state on 400 cores.

In the same way as in the 1D case (Gillet et al. 2023), at the interior boundary of each planet, we fixed the density to 1.326 × 10⁻¹⁰ g/cm³, the pressure to 12 µbar and the velocity to zero, and we assumed that the gas is entirely neutral. This pressure is typical of 1D models and ensures that at the base of the model atmosphere the gas remains in hydrostatic equilibrium and most of the stellar XUV energy is deposited at higher altitudes. The local temperature at the level of the atmosphere is 1100 K, which is broadly representative of conditions in close-in planets. We prescribed the same boundary conditions at the interior limit of each planet. The planet is now an internal boundary condition in our grid. The circular limit of this boundary condition is therefore discretised on the Cartesian grid. We have described how the edge through which the stellar wind enters is treated. Next, we describe the other three edges: the condition at the outer boundaries of the box on the other three edges is defined by a zero gradient for ρ, P and x_HI in all circumstances. In the case of velocity, we prescribed a zero gradient if the flow leaves the box and set it to a zero value if it enters the box at any of the three boundaries where the stellar wind is not imposed.

In this study we covered the same planetary masses as in Gillet et al. (2023): 0.02 0.05, 0.1, and 0.69 M_J (see Table 1). They range from a sub-Neptune planet to a Jovian-like planet similar to HD 209458b.

For the stellar XUV irradiation, we adopted a solar spectrum downloaded from the SOLID data exploitation project¹, as observed on December 13, 2021. We used R_★ = 7 × 10¹⁰ cm and M_★ = 1.98 × 10³⁰ g as the stellar radius and mass, respectively, and an orbital distance of 0.045 AU. The XUV-integrated flux of the star at 1 AU is 4.4 erg cm⁻² s⁻¹, and it is about 500 times larger at the orbit of the planet. We assumed the planet is tidally locked and that its spins and orbits have the same period.

Fig. 2 Thermal and ram pressure profiles for the M0.05 case. We compare between 1D (dashed lines) and 2D (solid lines) simulations. In 2D, the profiles are taken on the star-planet line on the day side and averaged over a time period of 42 h.

3 Self-consistent hydrodynamic escape in 2D

3.1 Dynamics of the planetary outflow

The planetary wind forms due to the interaction of incident XUV radiation from the host star with all layers of the planet’s atmosphere. It produces an extended envelope of partially ionised gas that escapes and is directed against the comparatively hot, fast, fully ionised stellar wind.

In our simulations, we captured the formation of multi-scale vortices for all our cases. These simulated instabilities only appear at high resolution (Tremblin & Chiang 2013). In all cases, the bow shock is created in front of the planet upward to the interaction of the two winds (stellar + planetary) regardless of the resolution. These dynamic characteristics are illustrated in Figure 1, for the M0.02 case on the left panels and for the M0.69 case on the right panels.

We observe that the temperature (top panels) can rise up to several million kelvins in the zone between the bow shock and the planet’s atmosphere, creating an extreme temperature gradient. The atmosphere of the least massive case extends up to several planetary radii in the line of sight of the star on the day side, with the creation of vortices lasting beyond 15 R_p in the cometary tail of the planet. Case M0.69 has a highly compressed atmosphere with a smoother planetary tail that seemingly has fewer vortices. In case M0.02, neutrals can be found up to 5 R_p in all directions (see the lower panels of Figure 1) whereas they are limited to a region much closer to the planet (in relative scales) in case M0.69. In addition, case M0.02 displays a more complex structure with a denser region of neutrals inside the tail.

The ram pressure at the sonic point of the stellar wind is one of the essential parameters that determine whether or not the exoplanetary atmosphere is confined by the stellar wind (Vidotto & Cleary 2020). If the ram pressure of the stellar wind is greater than the total pressure of the planetary wind (their ratio is equal to M² where M is the Mach number), the planterary wind remains subsonic and is therefore sensitive to the confinement by the stellar wind. In the case of our simulations, all our atmospheres are confined. We illustrate this in Figure 2, where we can see that below the interface between the subsonic and supersonic stellar wind, the total pressure (P_ram + P_T) of the prescribed stellar wind is greater (above 4 R_p) than that of the planetary wind (below 4 R_p). For the M0.05 case, deep in the atmosphere, the ram and thermal pressures of the planetary wind in 1D and 2D overlap up to 1.5 R_p (solid and dashed lines), these then differentiate, until completely diverging at the interface between the subsonic planetary wind and the subsonic stellar wind at 2.8 R_p.

Deep in the atmosphere, the pressure and temperature in the 2D simulations correspond to the equivalent 1D case (Gillet et al. 2023) as shown by the temperature-pressure profile in Figure 3. However, in the upper atmosphere, the 2D case deviates considerably from the 1D case, as expected, due to the deviation from the 1D steady state (represented by the black curve) and modulated by the interaction of the two winds leading to a compression of the planetary atmosphere. We note the small temporal variability on the T/P profiles vs r/R_p of all our cases, calculated over a time period around 42 hours corresponding to about half of the 3.5 days orbital period of the planet (pink area, in figure 3 for the case M0.05).

Fig. 3 Pressure and temperature profile for the M0.05 case. We compare the 1D profile (in black) with the 2D profile (in red) averaged over 42 hours. Temporal variability (in pink) is averaged over half of the orbital period. 1 Rp corresponds to the surface at which we start studying the upper atmosphere.

3.2 Secondary ionisation by photoelectrons in 2D

We next focused on the effect of photoelectrons, and compare the photoionisation and heating rates from our 2D simulations with those from the 1D geometry (Gillet et al. 2023). In this section, we illustrate it with a mini-Neptune type planet (case M0.05 in Table 1).

We first describe the case where secondary ionisation is not included (red curves) in the top panel of Figure 4. In this approach, all the excess energy above the ionisation threshold goes into heating. We observe that, in both the 1D (dotted lines) and 2D (solid lines) cases, the heating rates peak at similar altitudes. Both the 1D and 2D solutions have two maxima. The two maxima are more easily discernible in the 1D case. The 2D case decreases more slowly in the upper atmosphere up to 3.5 R_p, this is because the planetary wind is comparatively denser, as seen in Figure 3. Above this altitude, the ionisation fraction is close to 1, because the fully ionised stellar wind dominates, and the heating deposition rate naturally drops sharply. At low altitude, the photoionisation rates of the 1D and 2D cases increase similarly, up to 1.3 R_p.

At higher altitudes, the photoionisation rates of the 2D case (bottom panel) increase more rapidly, as a result of the compression exercised by the stellar wind. When secondary ionisation is taken into account (shown by the solid and dashed blue lines in Figure 4), we observe a shape similar to the thermal deposition rate curve in the two cases. The ionisation rate nevertheless increases more quickly in 2D between 1.5 R_p and 2.4 R_p before dropping due to the interaction with the stellar wind.

However, we observe significant differences in 2D when secondary ionisation is included (blue) or not (red). The collision point of the two subsonic flows is located at 2.6 R_p for the case of photoelectrons, while it is located much further up at 3.5 R_p when these are not included. The planetary wind is shaped by the balance between its thermal pressure and the ram pressure of the stellar wind. Since the total pressure is lower in the case that includes photoelectrons, due to lower loss rates and therefore lower velocities, the stellar wind compresses the planetary wind more. The 2D profile for the photoionisation rate deviates between 1.5 and 2.6 R_p (peak in the bottom panel of Figure 4) with a stronger rate in the 2D geometry due to a greater quantity of neutrals at this altitude. Finally, at the lowest altitudes, the geometry has very little impact on the heating and ionisation rates on the star-planet axis.

The ram (in red) and thermal (in blue) pressure profiles are also affected, as shown in Figure 5. The planetary ram pressure shows different structures when photoelectrons are included (dashed lines), while the thermal pressure decreases. The stellar wind shock moves closer to the planet when photoelectrons are considered, influencing the ram pressure structure deep in the atmosphere.

Fig. 4 Illustration of the heating rate along with the photoionisation rate and its dependence on the ionisation fraction, XHI , as a function of the radial distance on the day side and averaged over 42 hours. Top panel: comparison of the heating rate for the M0.05 case in 1D and 2D, with (in blue) and without (in red) the photoelectrons. Bottom panel: comparison of the photoionisation rate J×xHI for the M0.05 case in 1D and 2D, with (in blue) and without (in red) the photoelectrons.

Fig. 5 Ram (solid line) and thermal (dashed line) pressure profiles for case M0.05 with photoelectrons included or not, averaged over 42 hours. We show of the effect of secondary ionisation without (in red) and with (in blue) the pressures of planetary and stellar winds.

3.3 Atmospheric mass-loss rates

In the 2D geometry we are concerned with, the closed area over which the flux of mass must be calculated to determine the mass loss rate is the area of a cylinder. In a first step, the flux of mass is calculated through an annulus around the planet of radius, depending on the case, between 2 and 6 R_p . Using this range helps stabilise the mass-loss rate, meaning that it minimises the temporal variations in and out of the integration volume. One example is shown in Figure 6. In a second step, we multiplied the above determination by L=2R_p , where L represents the height of the cylinder. This choice ensures the match between the expressions for the mass loss rate between the 1D (spherical shell symmetry) and our current 2D geometry. We note that because the 1D and 2D geometries differ in their radial extension, the mass loss rates are not expected to match exactly between the two types of models even when using L = 2R_p . We report in Table 2 the mass-loss rate in 2D for all cases.

We observe a reduction in mass-loss rate from 22% (M0.69) to 48% (M0.02) when considering the effect of photoelectrons in the 2D cases. The same dependence on planetary mass as in 1D is shown: secondary ionisation impacts low-mass planets more easily. The variation in mass-loss is slightly less pronounced in 2D than in 1D, for example a reduction of 48% in 2D, compared to 54% in 1D for the M0.02 case.

We also observe a difference between the mass-loss rate in 1D and 2D, as reported in Table 3. The relative 1D-2D difference is reported as a percentage. In all cases, the mass-loss rate is slightly lower in 2D than in 1D. For example, the 2D M0.69 case presents a mass-loss rate of 91% of the analogous 1D case without photoelectrons. This is probably due that, in 2D, only the day part of the atmosphere receives the incident flow. In addition, we have shown that the planetary wind is confined and compressed in our 2D models, which a priori leads to a moderate reduction in the mass loss rate, as predicted in Vidotto & Cleary (2020).

We ran one more simulation with the 2D geometry, without stellar wind, and including the specifics of photoelectrons in order to assess the 1D versus 2D effects. We obtained a massloss rate of 5.0 ×10¹⁰ g/s, which should be compared with the 1.78 ×10¹⁰ g/s for the case with stellar wind. This shows that, in the 2D case, stellar wind confinement reduces the mass-loss rate by a factor of 3.

Fig. 6 Illustration of the calculation of the mass-loss rate in 2D for the M0.1 case. The mass flux integration ring is represented by the red ring.

4 Lyα absorption of planetary atmospheric escape

4.1 Theory of Lyα

The Lyα absorption profile can provide insight into the magnitude of the atmospheric escape process (Owen et al. 2023). The relevant information is encoded in the signal strength and its wavelength dependence. We used our simulation results to calculate synthetic transmission spectra. The absorption cross-section of a line is: $$\sigma (\lambda )=\frac{\pi e^2}{m_e{c}^2}{f}_{lu}{\lambda }_0^2\psi (\lambda ,u,T),$$(7)

where e is the electronic charge, m_e the electronic mass, c the speed of light, λ₀ is the wavelength of the line at rest equal to 1215.67 Å and f_lu is the strength of its oscillator 0.41641 for the Lyα line (Kramida et al. 2018).

The so-called Voigt profile ψ(λ, u, T) (normalised so that ∫ ψ(λ, u, T)dλ = 1) is a convolution of a Lorentzian profile resulting from the natural widening (which dominates in the wings) and a Gaussian profile resulting from the Doppler shift, dominating in the centre of the line. It is calculated by Doppler effect (λ-λ₀)(1 − u/c) where u is the velocity component along a single ray in the star-planet line of sight.

The Voigt distribution is calculated with width at half maximum for the Doppler component $\text{Δ}{\lambda }_D={\lambda }_0\times 7.16\times {10}^{-7}\sqrt{T}$, with T the local temperature in Kelvins. The width at half maximum for the Lorentzian component is $\text{Δ}{\lambda }_L={\lambda }_0^2{A}_{ul}/(2\pi c)$, where A_ul=6.2648 × 10⁸ s⁻¹ is the transition probability of the line (Kramida et al. 2018).

The differential optical thickness at wavelength λ for a ray passing through a grid element along the line of sight from the star to the planet (for all x component in the box) is calculated as $$\text{Δ}\tau (\lambda ,x,y)=n_{\text{HI}}(x,y)\sigma (\lambda )dx,$$(8)

where n_H I is the neutral particle density of hydrogen and σ(λ) is the absorption cross-section. The optical thickness τ_λ(y) at a given position in the y direction in our 2D spatial grid is calculated as follows: $${\tau }_{\lambda }(y)={\displaystyle \sum _x\text{Δ}}\tau (\lambda ,x,y).$$(9)

Finally, we calculated the effective size of the planet for a given wavelength, λ, in our 2D geometry as follows: $$D_{\lambda }={\displaystyle {\int }_{-\infty }^{+\infty }\left[1-\exp \left(-{\tau }_{\lambda }(y)\right)\right]}\text{d}y.$$(10)

To construct the synthetic Lyα profile, we chose 100 wavelength intervals over the wavelength range [1214.17–1217.17 Å]. We verified that this choice did not impact our results by increasing the spectral resolution by a factor of 2 in a few selected tests. We then calculated the wavelength-dependent effective radius as R_eff = D_λ/2 and compared it to the solid optical radius, R_opt, corresponding to a planet without atmosphere. We then calculated the wavelength-dependent effective radius as R_eff = D_λ/2 and compared it to the solid optical radius, R_opt, corresponding to a planet without atmosphere (what is measured with optical photometry).

Fig. 7 Illustration of the Lyα signal. Top panel: 2D representation of the instantaneous number density of gas escaping from a Jupiter-like planet. The coloured hatching shows the different regions scanned in our simulation grid to obtain the Lyα profile. Lower panel: results of instantaneous Lyα profiles associated with the different regions scanned around the planet.

4.2 Origin of the signal for a hot Jupiter

To better understand the origin of the Lyα signal, we illustrated the atmosphere of the M0.69 case. As shown in Equation (8), the signal is sensitive to the particle density of neutral hydrogen n_HI , the temperature and the velocity of the gas (both contained in ψ(λ, u, T)). We expect a strong signal contribution near the planet, where neutral densities are high, and no contribution in regions where the fully ionised stellar wind dominates. As for the shape of the absorption line, it will be partly dictated by the structure of the neutral gas flow velocity.

We split the signal R_eff/R_opt into several contributions in the y direction near the planet (see the top panel of Fig. 7). As expected, we observe (bottom panel) that two regions contribute the most to the Lyα signal, those in magenta (b) and green (c) colours. In the same panel, the total signal is represented by the black line. The two regions scan the planetary atmosphere where the particle density of neutral hydrogen is the highest as well as the planetary tail where we can still find neutrals up to 20 R_p in the positive x direction.

The contribution of zone (b) is saturated at the core of the line, contributing to the total signal by the geometric size of the scanned region. As expected, there is a very small contribution from region (a) and region (d), with less than 0.01% of the total signal, both regions combined. Ionised stellar wind dominates in these regions, preventing the build-up of a significant Lyα signal. On the other hand, the wings of the distribution are dominated by the purple part (b) of the domain due to a greater density of neutrals and the temperature of the gas.

Figure 8 shows the velocity v_x of the flow. The material escapes with relatively modest velocities of 20–40 km/s. However, at the interface between the planetary and stellar wind in the planet’s wake, the velocities reach 50–100 km/s (dotted black contour in Figure 8). This acceleration is caused by the interaction with the stellar wind. To assess why such velocities are reached, we went back to the extra model that we ran without stellar wind. This model (not shown here) presents absorption at high velocities of up to |v| ≃ 100 km/s, but without a clearly preferred blue or red shift. The large blueshift seen in the figures therefore originates from the interaction of the planetary wind with the stellar wind, and this interaction helps increase speeds closer to the planet due to entrainment.

The peaks of the two dominant contributions (regions b and c in the upper panel of Figure 7 are all blueshifted due to the stellar wind accelerating the gas to high velocities towards the planetary tail and hence towards the observer. The total signal is blueshifted, resulting in a particularly noticeable contribution at a speed of 100 km/s on the left wing of the line. Unlike the 1D case where we observe lower speeds, here, the speeds observed in 2D are caused by the entrainment of neutrals by the stellar wind. The asymmetry of the black line is due to the presence of patches of gas moving towards the star, creating a redshifted component at 20 km/s in the zone (b) as seen in Figure 8.

To understand how well the gas absorbs, we calculated the absorbance as follows: ${\kappa }_{{\lambda }_k}=n_H(x,y)\sigma \left(x,y,{\lambda }_k\right)$ (also equal to dτ(x, y, λ_k)/dx). It is a property of the gas, independent of the numerical grid adopted. Figure 9 illustrates the absorbance ${\kappa }_{{\lambda }_k}$ at different wavelengths represented by the pink, orange, and cyan dashed curves in Figure 7. We observe that the absorbance is higher on the short-wavelength wing of the Lyα at λ = 1214.47 Å compared to the long-wavelength wing at λ = 1216.89 Å. Near the centre of the line at λ = 1215.68 Å, the gas absorbs on a larger spatial scale, with an additional signal coming from the interface between the planetary and stellar wind.

Fig. 8 Profile of the 42h-averaged velocity , vx, computed at the inertial frame for the case M0.05 along with a velocity line of 100 km/s, shown by the dotted black contour.

4.3 Impact of the planetary mass on the Lyα profile

The structure of the planetary atmosphere depends on the stellar parameters and on how the deposited energy is used by the planetary wind to overcome the gravitational potential of the planet. Low-mass planets have atmospheres that are more extended than more massive planets than higher planet mass planets. This can impact the shape and amplitude of the Lyα profile. This can impact the shape and amplitude of the Lyα profile. Figure 10 shows the Lyα profile (now normalised to R_★) for our four different cases: M0.69, M0.1, M0.05 and M0.02 described in Table 1.

As expected, the effective Lyα radius of the planet normalised to the stellar radius (R_★ =7 × 10¹⁰ cm) is larger for the more massive planet in the core of the line, and goes down to about 0.5 on average for M0.02. Interestingly, all profiles are blueshifted with similar speeds ranging from 50 km/s in the core of the line and up to 100 km/s on the left wing of the line where it is observable. The lifetime of hydrogen in the planetary tail depends on the rate of photoionisation but also on the dynamics of the winds (planetary + stellar). To observe neutral hydrogen at such high speeds, understanding the role of the stellar wind which accelerates the planetary wind, is essential.

We observe asymmetries at the heart of the line that differ from one case to another. The ‘irregularities’ observed in case M0.02 are caused by the presence of patches of neutral gas around the planet. The instantaneous shape of the absorption line is dictated by the transient bursts of the planetary wind rather than by the characteristics of the steady-state wind. If the planetary wind is more turbulent, the asymmetry at the heart of the line will increase.

5 Effect of secondary ionisation on Lyα

We next studied the impact of secondary ionisation on the Lyα profile by analysing the two extreme cases of our sample: the cases M0.69 (in blue) and M0.02 (in red). We focused only on the qualitative aspect: a comparison with observations will be made in a future study.

We observe that when the secondary ionisation is included (dashed line in the Figure 11), the amplitude of the Lyα is generally reduced at positive velocities (those more closely connected with the dayside of the planet). The reason for this reduction is that the size of the neutral cloud is reduced as the ionisation proceeds faster. The same effect seems to occur on the night side of the M0.02 simulation. The situation for the M0.69 case is somewhat more complex. Looking into the details of the flowfield, we find that in the simulation with the detailed photoelectron, there is a reduced number density of neutrals on the nightside but higher velocities. The combination of these two features appears as the broadening of the Lyα absorption signal at negative velocities in Figure 11). This serves as a reminder that multiple physical processes can combine differently in ways that are difficult to anticipate.

Both cases M0.02 and M0.69 show a decrease of 0.1 R_eff/R_* around the core of the line. However, on the left wing, the signal is increased by 0.1 R_eff/R_* for M0.69 while decreased by 0.2 for the case M0.02. The effect of photoelectrons is therefore even more important for low-mass planets. We observe an increase in the signal at shorter wavelengths for the case M0.69 when photoelectrons are included, which is the spectral region that is observable. The origin of this interesting feature is not yet fully clear, and a more detailed investigation will be carried out in a future work to assess its significance at these and other planets for which Lyα observations have been made. Faster ionisation when photoelectrons are included tends to shrink the cloud of neutrals that absorb the Lyα photons of the star. The change in the mass-loss rate, however, has additional implications on the dynamics of the atmosphere, which may affect the core and the wings of the absorption line differently.

Fig. 9 Instantaneous absorbance coefficient at different wavelengths. From left to right: absorbance for λ = 1214.47 Å, λ = 1215.68 Å and λ = 1216.89 Å, corresponding to different parts of the Lyα profile shown in Figure 7 (from left to right, the panels correspond to the magenta, orange and cyan vertical lines of Figure 7).

Fig. 10 Illustration of the impact of the planetary mass on the Lyα signal for the different cases normalised with the Sun radius, averaged over 42 hours.

6 Impact of stellar flares on Lyα

Stellar flares are frequent on young low-mass stars, notably for spectral types M or K around which many detected exoplanets orbit, but also on more massive stars like our Sun of spectral type G (Hawley et al. 2014). These events are very energetic, impulsive and difficult to predict, and it is potentially important to take into account their influence on exoplanetary atmospheres when interpreting transit observations.

Very few studies or modelling of flares and their impact on exoplanetary atmospheres have been carried out. Hazra et al. (2022) studied the impact of these flares on HD 189733b’s atmosphere using 3D radiation hydrodynamic simulations of atmospheric exhaust that include heating by stellar photons using a constant flare model.

To determine if a dynamic event originating from the star influences the loss of atmospheric mass, we modelled and parameterised the flare as time-varying multiplicative factor ƒ_λ(t) applied to the X-ray part of the stellar spectrum and defined as follows: $$f_{\lambda }(t)=\{\begin{array}{ll}1\hfill & t<t_0\hfill \\ 20\times e^{\left(t-t_0\right)}\hfill & t_0<t<t_1\hfill \\ 10\times e^{-2.06\times {10}^{-4}\left(t-t_1\right)}\hfill & t_1<t<t_2\hfill \\ 1\hfill & t>t_2\hfill \end{array},$$(11)

with t in seconds. It increases rapidly in order to reach at t₁ = 300 s an X-ray spectrum 10 times greater, then a slow decrease until returning at t₂ = 3 h to its initial state. We only multiplied the incoming stellar flux by f_λ on the most energetic segments of the solar spectrum λ < 100 Å corresponding to the X-ray range (see Figure 12).

We next looked at the impact of the flare on the Lyα profile shown in Figure 13 for cases M0.02 (upper panel), and M0.69 (lower panel). The Lyα profile at the peak of the flare is modelled by the dashed lines. We observe, in all cases, that the effect of the flare is located within the natural temporal variability of the Lyα profile (pink and grey zones, temporal average of the profiles over 42 hours) whether secondary ionisation is included or not.

This means that we would not be able to directly observe an flare on these profiles. If we calculate the shortest heating timescale (by comparing the total energy E = P_T/(γ − 1) + ρu²/2 with the energy allocated to heating H, t_heat = min(E/H)) for the M0.02 case, we find t_heat = 51.08 h and for the M0.69 case t_heat = 9.89 h, during the peak of the flare. The time of the flare is therefore too short compared to the characteristic heating times of the planet to disturb the atmosphere in all cases. Therefore, we do not observe a change in the mass-loss rate. Nevertheless, it is possible that flares of greater intensity, or that a fast repetition of closely flares could temporarily modify the structure of the atmosphere, especially for the most massive planets due to their shorter heating timescales. Our results differ from Hazra et al. (2022), where they maintain the conditions for a longer period of time in order to let their model reach steady-state. We, however, considered a realistic temporal evolution of the XUV flux of a flare.

Fig. 11 Illustration of the impact of secondary ionisation on the Lyα profile for cases M0.02 (in red) and M0.69 (in blue). Solid lines represent cases without photoelectrons, while dashed lines represent cases with photoelectrons normalised with the Sun radius. These profiles are also averaged over a period of 42 hours.

Fig. 12 Modelling of the flare on the solar spectrum used in our simulations.

7 Conclusion

After a first study in 1D (Gillet et al. 2023) in which we determined its limitations due to the simplified geometry, we extended our model to a more realistic 2D geometry that integrates both a day part and a night part for the atmosphere. Thanks to this 2D geometry, we were able to evaluate the impact of photoelectrons on the Lyα absorption, which is one of the main observables of atmospheric escape. We simulated planets with the same masses as the 1D cases, and studied their effect on heating, ionisation, and excitation while comparing them to the 1D results. We also studied the interaction between the stellar wind and the planetary wind.

Planetary wind is formed by the interaction of the host star’s XUV radiation with the planet’s atmosphere, creating an envelope of partially ionised gas directed towards the ionised stellar wind. The interaction of the two winds strikes a balance between thermal pressure and ram dynamic pressure. We are able to see the classic features of multi-scale vortex formation. This interaction forms a shock in front of the planet, whatever the resolution of the computational grid, with temperatures reaching several million kelvins at the level of this interaction zone due to compression of the gas.

In the absence of secondary ionisation, the heating rate peaks at similar altitudes in the 1D and 2D models, but the 2D behaviour is more complex. Photoionisation is primarily affected at higher altitudes in 2D compared to 1D, with significant differences when secondary ionisation is included. Overall, the ionisation rate is higher in 2D due to a greater amount of neutrals at these altitudes. On the other hand, the lowest altitudes are only slightly impacted by the geometry of our simulations. The mass-loss rates are lower in 2D models than in the equivalent 1D models, regardless of whether photoelectrons are considered. The inclusion of photoelectrons decreased the mass-loss rate by 22–48% in the studied sample, similar to what was observed in 1D (Gillet et al. 2023), showing the impact secondary ionisation has on low-mass planets.

To link the simulations with observations, we examined the origin of the Lyα signal in the atmosphere of case M0.69. We notice that the signal depends on the density of neutral hydrogen particles, the temperature, and the velocity of the gas. The signal contribution is significant near the planet, where neutral densities are high, but zero in regions dominated by the ionised stellar wind. The shape of the absorption line is related to the structure of the neutral gas flow. Our analysis shows that the largest contributions to the Lyα signal come from close to the planet. However, the high velocities of the gas towards the observer induce a blueshift of the signal, notably at 100 km/s, which is greater than the velocities observed in 1D and is linked to the interaction with the stellar wind. We do not take energetic neutral atoms into consideration, but a future work will study the charge exchange between the two winds.

The structure of an exoplanet’s atmosphere depends on the stellar radiative output and on the use of that energy by the atmospheric gas to escape the gravitational potential. In our simulations, low-mass planets have a more extended atmosphere than those with higher masses relative to their radii. This directly affects the Lyα profile. All profiles are blueshifted at similar speeds, but asymmetries at the core of the line vary from one case to another: they are linked to the dynamics of planetary and stellar winds.

We have discussed the qualitative aspects of how secondary ionisation affects the Lyα profile. We generally observe a decrease in the amplitude of the Lyα signal, which is more significant for the low-mass (M0.02) case. As photoelectrons increase the ionisation rate of the atmosphere, thereby reducing the presence of neutrals in the upper atmosphere, this directly affects the Lyα profile.

Stellar flares are common on young, low-mass stars of spectral type M or K, as well as on more massive G-type stars such as our Sun. These events, which are energetic and difficult to predict, can affect the atmospheres of transiting exoplanets. We examined the impact of a solar flare on two simulated Jovian- and mini-Neptune-mass planets orbiting a solar-like star. We find that the effect of the flare lies within the natural variability of the Lyα profile, even when considering secondary ionisation. The profiles do not show a significant change in mass-loss rates. However, we believe that more intense or fast and repeated flares could temporarily change the structure of the atmosphere.

Fig. 13 Illustration of the impact of a flare (dashed line) on the Lyα profiles of cases M0.02 (upper panel) and M0.69 (lower panel) with (in red) and without (in black) photoelectrons, averaged over 42 hours.

Acknowledgements

We acknowledge funding from the Programme National de Planétologie (INSU/PNP). A.S. acknowledges funding from the European Union’s Horizon-2020 research and innovation programme (grant agreement no. 776403 ExoplANETS-A) and the PLATO/CNES grant at CEA/IRFU/DAp, the French Agence Nationale de la Recherche (ANR) project STORMGENESIS #ANR-22-CE31-0013-01, and the European Research Council project ExoMagnets (grant agreement no. 101125367).

References

All Tables

All Figures

	Fig. 1 Illustration of the dynamic characteristics of the interaction of the planetary wind (with the planet at the centre delimited by a circle) and stellar wind for the most extreme cases M0.02 (left) and M0.69 (right). Top panel: instantaneous temperature (in Kelvins) in logarithmic scale. Bottom panel: instantaneous density of neutrals in logarithmic scale.
In the text

	Fig. 2 Thermal and ram pressure profiles for the M0.05 case. We compare between 1D (dashed lines) and 2D (solid lines) simulations. In 2D, the profiles are taken on the star-planet line on the day side and averaged over a time period of 42 h.
In the text

	Fig. 3 Pressure and temperature profile for the M0.05 case. We compare the 1D profile (in black) with the 2D profile (in red) averaged over 42 hours. Temporal variability (in pink) is averaged over half of the orbital period. 1 Rp corresponds to the surface at which we start studying the upper atmosphere.
In the text

Fig. 4 Illustration of the heating rate along with the photoionisation rate and its dependence on the ionisation fraction, XHI , as a function of the radial distance on the day side and averaged over 42 hours. Top panel: comparison of the heating rate for the M0.05 case in 1D and 2D, with (in blue) and without (in red) the photoelectrons. Bottom panel: comparison of the photoionisation rate J×xHI for the M0.05 case in 1D and 2D, with (in blue) and without (in red) the photoelectrons.

In the text

	Fig. 5 Ram (solid line) and thermal (dashed line) pressure profiles for case M0.05 with photoelectrons included or not, averaged over 42 hours. We show of the effect of secondary ionisation without (in red) and with (in blue) the pressures of planetary and stellar winds.
In the text

	Fig. 6 Illustration of the calculation of the mass-loss rate in 2D for the M0.1 case. The mass flux integration ring is represented by the red ring.
In the text

	Fig. 7 Illustration of the Lyα signal. Top panel: 2D representation of the instantaneous number density of gas escaping from a Jupiter-like planet. The coloured hatching shows the different regions scanned in our simulation grid to obtain the Lyα profile. Lower panel: results of instantaneous Lyα profiles associated with the different regions scanned around the planet.
In the text

	Fig. 8 Profile of the 42h-averaged velocity , vx, computed at the inertial frame for the case M0.05 along with a velocity line of 100 km/s, shown by the dotted black contour.
In the text

	Fig. 9 Instantaneous absorbance coefficient at different wavelengths. From left to right: absorbance for λ = 1214.47 Å, λ = 1215.68 Å and λ = 1216.89 Å, corresponding to different parts of the Lyα profile shown in Figure 7 (from left to right, the panels correspond to the magenta, orange and cyan vertical lines of Figure 7).
In the text

	Fig. 10 Illustration of the impact of the planetary mass on the Lyα signal for the different cases normalised with the Sun radius, averaged over 42 hours.
In the text

	Fig. 11 Illustration of the impact of secondary ionisation on the Lyα profile for cases M0.02 (in red) and M0.69 (in blue). Solid lines represent cases without photoelectrons, while dashed lines represent cases with photoelectrons normalised with the Sun radius. These profiles are also averaged over a period of 42 hours.
In the text

	Fig. 12 Modelling of the flare on the solar spectrum used in our simulations.
In the text

	Fig. 13 Illustration of the impact of a flare (dashed line) on the Lyα profiles of cases M0.02 (upper panel) and M0.69 (lower panel) with (in red) and without (in black) photoelectrons, averaged over 42 hours.
In the text