<figure><img src="image_1.png" /></figure>

HAL
open science

Secondary ionisation in hot atmospheres and
interactions between planetary and stellar winds

A. Gillet, A. Strugarek, A. García Muñoz

► To cite this version:

A. Gillet, A. Strugarek, A. García Muñoz. Secondary ionisation in hot atmospheres and interactions between planetary and stellar winds. Astronomy & Astrophysics - A&A, 2025, 696, pp.A64. 10.1051/0004-6361/202452559. cea-05021854

HAL Id: cea-05021854

https://cea.hal.science/cea-05021854v1

Submitted on 5 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Secondary ionisation in hot atmospheres and interactions between planetary and stellar winds

A. Gillet*, A. Strugarek, and A. García Muñoz

Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

Received 10 October 2024 / Accepted 4 February 2025

## ABSTRACT

**Context.** The loss of close-in planetary atmospheres is influenced by various physical processes, such as photoionisation, which could potentially affect the atmosphere survivability on a secular timescale. The amount of stellar radiation converted into heat depends on the energy of the primary electrons produced by photoionisation and the local ionisation fraction. The HI Lyman-alpha (Lyα) line is an excellent probe for atmospheric escape, but the reason for the origin of the high velocities detected in this line is not yet fully understood.

**Aims.** We characterise the effect of secondary ionisation by photoelectrons on the ionisation and heating of the gas for different planet-star systems in a 2D geometry. We study the interaction between the planetary and stellar winds, the difference in the predicted mass-loss rates between 1D and 2D models, the Lyα signal, and the impact of stellar flares.

**Methods.** Using the PLUTO code, we performed 2D hydrodynamic simulations for four planets in the Neptune to Jupiter size range. We produced synthetic Lyα profiles to determine the origin of the signal and in particular its high-velocity Doppler shift.

**Results.** We find a similar trend to that shown in the 1D models, with a decrease in the planetary mass-loss rate for all systems when secondary ionisation is taken into account. The mass-loss rates decrease by 48% for the least massive planet when secondary ionisation is accounted for. However, the decrease is less pronounced in 2D than in 1D. We observe differences in the Lyα profile between the different cases and significant asymmetries in all of them, especially for the lower-mass planets. Finally, we observe that stellar flares do not affect the mass-loss rate because they act, in general, on a timescale that is too short.

**Conclusions.** We find and confirm in our extended 2D model that photoelectrons affect mass-loss rates by factors that are potentially important for planetary evolution theories, and that they also affect the Lyα profile. We find velocities in the escaping atmosphere of up to 100 km/s, with the gas moving away from the star, which could be the result of the interaction with the stellar wind. We find that stellar flares generally occur on a timescale that is too short to have a visible impact on the mass-loss rate of the atmosphere.

**Key words.** atomic processes – hydrodynamics – planets and satellites: atmospheres – planet-star interactions – stars: winds, outflows

## 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, charge-exchange 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

<footnote>* Corresponding author; alexandre.gillet@cea.fr</footnote>
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 \mathbf{u}) = 0, \qquad (1) $$

$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla \cdot [\rho\mathbf{u}\mathbf{u} + P_T \mathbf{I}] = -\rho\nabla\phi + \rho(\mathbf{F}_{\text{cent}} + \mathbf{F}_{\text{cor}}), \qquad (2) $$

$$ \frac{\partial(E + \rho\phi)}{\partial t} + \nabla \cdot [(E + P_T + \rho\phi)\mathbf{u}] = \rho\mathbf{u}(\mathbf{F}_{\text{cent}} + \mathbf{F}_{\text{cor}}) + H - C, \qquad (3) $$

where $t$ is the time, $\mathbf{u}$ is the fluid velocity, and $\rho$ is the total density, $P_T$ is the thermal pressure, $E = P_T/(\gamma - 1) + \rho\mathbf{u}^2/2$ is the

total energy with $\gamma = 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

$$ \phi = -\frac{GM_p}{r} - \frac{GM_*}{R_{\text{orbit}} - r}, \qquad (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 $\mathbf{F}_{\text{cor}} = -2(\mathbf{\Omega} \times \mathbf{u})$ and $\mathbf{F}_{\text{cent}} = \mathbf{\Omega} \times (\mathbf{\Omega} \times \mathbf{R})$ where $\mathbf{R}$ is the position vector relative to the star and $\mathbf{\Omega} = \sqrt{(GM_*/R_{\text{orb}}^3)} \mathbf{e}_z$ is the orbital rotation rate with $\mathbf{e}_z$ is the unit vector along the coordinate $z$. The model presented here solves these equations on the $(x,y)$ orbital plane.

We consider a planetary atmosphere composed of atomic hydrogen in neutral HI 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^{-1}$] is given by:

$$ J = \int_{\lambda_{\text{min}}}^{\lambda_0} \sigma_{\lambda} F_{\star} \left(\frac{\lambda}{hc}\right) (1 + \Phi_{\lambda,xe}) d\lambda, \qquad (5) $$

where $F_*$ [$\text{erg cm}^{-2} \text{s}^{-1} \text{Å}^{-1}$] is the attenuated stellar flux. $\Phi_{\lambda,xe}$ is the number of secondary ions created per photoionisation. The heating deposition rate $H$ [$\text{erg s}^{-1} \text{cm}^{-3}$] in the energy equation can be expressed as:

$$ H = n_{\text{HI}} \int \sigma_{\lambda} F_{\star} \left(1 - \frac{\lambda}{\lambda_0}\right) \eta_{\lambda,xe} d\lambda \qquad (6) $$

with $\eta_{\lambda,xe}$ 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_{\text{ram}} = \frac{1}{2}\rho v^2$. If the
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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 $\dot{M}_{\odot}$ with $\dot{M}_{\odot} = 2 \times 10^{-14} \dot{M}_{\odot} \text{yr}^{-1}$ (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 \times 10^6 \text{ K}$. The density $\rho$ is prescribed by the stellar mass-loss as $\rho_w = M_{\odot}/4\pi^2 v_w$ with $v_w = 400 \text{ 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.

## 2.3. Numerical model

We solved the Euler equations on a 2D Cartesian grid, $x, y = [-20, 20] \text{ R}_p$ parallel to the orbital plane with the planet located at centre of the box. We carried out high resolution 2D simulations ($5000 \text{ cells} \times 5000 \text{ cells}$) cut into a uniform grid of $1000 \times 1000 \text{ cells}$ covering the square $[-1,1]$ around the planet and
<table><caption>**Table 1.** Planet parameters for planetary systems M0.69-M0.02.</caption><thead><tr><th>Parameters</th><th>M0.69</th><th>M0.1</th><th>M0.05</th><th>M0.02</th></tr></thead><tbody><tr><td>M<sub>p</sub> (M<sub>1</sub>)</td><td>0.69</td><td>0.10</td><td>0.05</td><td>0.02</td></tr><tr><td>R<sub>p</sub> (R<sub>J</sub>)</td><td>1.32</td><td>0.69</td><td>0.55</td><td>0.40</td></tr></tbody></table>

**Notes.** $M_p$ and $R_p$ are the planetary mass and radius respectively in Jovian units.

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 \times 10^{-10}$ g/cm$^3$, the pressure to $12 \mu\text{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 $\rho$, P and $x_{H1}$ 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<fnref n="1" />, as observed on December 13, 2021. We used $R_* = 7 \times 10^{10}$ cm and $M_* = 1.98 \times 10^{30}$ 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$^{-2}$ s$^{-1}$, and it is about 500 times

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

larger at the orbit of the planet. We assumed the planet is tidally locked and that its spins and orbits have the same period.

## 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^2$ where M is the Mach number), the planetary wind remains subsonic and is therefore sensitive to the confinement by the stellar wind. In the case of our simulations, all our

<footnote n="1">¹ European comprehensive solar irradiance data exploitation; https://www.mps.mpg.de/solar-variability/solid</footnote>
<figure><img src="image_5.png" /><figcaption>**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 $R_p$ corresponds to the surface at which we start studying the upper atmosphere.</figcaption></figure>

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).

## 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,

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Illustration of the heating rate along with the photoionisation rate and its dependence on the ionisation fraction, $X_{HI}$, 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 \times X_{HI}$ for the M0.05 case in 1D and 2D, with (in blue) and without (in red) the photoelectrons.</figcaption></figure>

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
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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.

### 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

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Mass-loss rates.</caption><thead><tr><th>Planet</th><th>Photoelectrons not included</th><th>Photoelectrons included</th><th>Diminution</th></tr></thead><tbody><tr><td>M0.69</td><td>2.31</td><td>1.78</td><td>22%</td></tr><tr><td>M0.1</td><td>5.16</td><td>3.43</td><td>33%</td></tr><tr><td>M0.05</td><td>4.96</td><td>2.79</td><td>43%</td></tr><tr><td>M0.02</td><td>5.36</td><td>2.80</td><td>48%</td></tr></tbody></table>

**Notes.** Atmospheric mass-loss rates [$10^{10}$ g/s] calculated in 2D for each of the exoplanetary systems in this work, with photoelectrons included or not.

<table><caption>**Table 3.** Differences between 1D and 2D.</caption><thead><tr><th rowspan="2">Planet</th><th colspan="2">2D/1D</th></tr><tr><th>Photoelectrons not included</th><th>Photoelectrons included</th></tr></thead><tbody><tr><td>M0.69</td><td>91%</td><td>88%</td></tr><tr><td>M0.1</td><td>81%</td><td>76%</td></tr><tr><td>M0.05</td><td>82%</td><td>80%</td></tr><tr><td>M0.02</td><td>83%</td><td>81%</td></tr></tbody></table>

**Notes.** Relative difference in atmospheric mass-loss between 2D and 1D with photoelectrons included or not. 100% corresponds to the same value between 2D and 1D, a lower percentage indicates a decrease in 2D.

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 mass-loss 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.

# 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), \qquad (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 $\psi(\lambda, u, T)$ (normalised so that $\int \psi(\lambda, u, T)d\lambda = 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 ($\lambda-\lambda_0$)(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 $\Delta\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 $\Delta\lambda_L = \lambda_0^2 A_{ul} / (2\pi c)$, where $A_{ul}=6.2648 \times 10^8 \text{ s}^{-1}$ is the transition probability of the line (Kramida et al. 2018).

The differential optical thickness at wavelength $\lambda$ 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

$$ \Delta\tau(\lambda, x, y) = n_{HI}(x, y) \sigma(\lambda) dx, \qquad (8) $$

where $n_{HI}$ is the neutral particle density of hydrogen and $\sigma(\lambda)$ is the absorption cross-section. The optical thickness $\tau_{\lambda}(y)$ at a given position in the y direction in our 2D spatial grid is calculated as follows:

$$ \tau_{\lambda}(y) = \sum_{x} \Delta\tau(\lambda, x, y). \qquad (9) $$

Finally, we calculated the effective size of the planet for a given wavelength, $\lambda$, in our 2D geometry as follows:

$$ D_{\lambda} = \int_{-\infty}^{+\infty} [1 - \exp(-\tau_{\lambda}(y))] dy. \qquad (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

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

$R_{eff} = D_{\lambda}/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_{\lambda}/2$ and compared it to the solid optical radius, $R_{opt}$, corresponding to a planet without atmosphere (what is measured with optical photometry).

## 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 $\psi(\lambda, 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
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Profile of the 42h-averaged velocity, $v_x$, 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.</figcaption></figure>

the most to the $Ly\alpha$ 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\alpha$ 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| \approx 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(x, y, \lambda_k)$ (also equal to

$d\tau(x, y, \lambda_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\alpha$ at $\lambda = 1214.47$ Å compared to the long-wavelength wing at $\lambda = 1216.89$ Å. Near the centre of the line at $\lambda = 1215.68$ Å, the gas absorbs on a larger spatial scale, with an additional signal coming from the interface between the planetary and stellar wind.

### 4.3. Impact of the planetary mass on the $Ly\alpha$ 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\alpha$ profile. This can impact the shape and amplitude of the $Ly\alpha$ profile. Figure 10 shows the $Ly\alpha$ 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\alpha$ radius of the planet normalised to the stellar radius ($R_* = 7 \times 10^{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\alpha$

We next studied the impact of secondary ionisation on the $Ly\alpha$ 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\alpha$ 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\alpha$ absorption signal at negative velocities in Figure 11). This serves as a reminder that
<figure><img src="image_13.png" /><figcaption>Fig. 9. Instantaneous absorbance coefficient at different wavelengths. From left to right: absorbance for $\lambda = 1214.47$ Å, $\lambda = 1215.68$ Å and $\lambda = 1216.89$ Å, corresponding to different parts of the Ly$\alpha$ profile shown in Figure 7 (from left to right, the panels correspond to the magenta, orange and cyan vertical lines of Figure 7).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10. Illustration of the impact of the planetary mass on the Ly$\alpha$ signal for the different cases normalised with the Sun radius, averaged over 42 hours.</figcaption></figure>

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_{\text{eff}}/R_*$ around the core of the line. However, on the left wing, the signal is increased by 0.1 $R_{\text{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$\alpha$ observations have been made. Faster ionisation when photoelectrons are included tends to shrink the cloud of neutrals that absorb the Ly$\alpha$ 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.

## 6. Impact of stellar flares on Ly$\alpha$

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

<figure><img src="image_15.png" /><figcaption>Fig. 11. Illustration of the impact of secondary ionisation on the Ly$\alpha$ 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.</figcaption></figure>

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 $f_{\lambda}(t)$ applied to the X-ray part of the stellar spectrum and defined as follows:

$$f_{\lambda}(t) = \begin{cases} 1 & t < t_0 \\ 20 \times e^{(t-t_0)} & t_0 < t < t_1 \\ 10 \times e^{-2.06 \times 10^{-4}(t-t_1)} & t_1 < t < t_2 \\ 1 & t > t_2 \end{cases} \quad (11)$$

with $t$ in seconds. It increases rapidly in order to reach at $t_1 = 300$ s an X-ray spectrum 10 times greater, then a slow decrease
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Modelling of the flare on the solar spectrum used in our simulations.</figcaption></figure>

until returning at $t_2 = 3$ h to its initial state. We only multiplied the incoming stellar flux by $f_\lambda$ on the most energetic segments of the solar spectrum $\lambda < 100$ Å corresponding to the X-ray range (see Figure 12).

We next looked at the impact of the flare on the Ly$\alpha$ profile shown in Figure 13 for cases M0.02 (upper panel), and M0.69 (lower panel). The Ly$\alpha$ 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$\alpha$ 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/(\gamma-1)+\rho u^2/2$ with the energy allocated to heating H, $t_{\text{heat}} = \min(E/H)$) for the M0.02 case, we find $t_{\text{heat}} = 51.08$ h and for the M0.69 case $t_{\text{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.

## 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$\alpha$ 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

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 13.** Illustration of the impact of a flare (dashed line) on the Ly$\alpha$ profiles of cases M0.02 (upper panel) and M0.69 (lower panel) with (in red) and without (in black) photoelectrons, averaged over 42 hours.</figcaption></figure>

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$\alpha$ 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 densi-
ties are high, but zero in regions dominated by the ionised stellar
wind. The shape of the absorption line is related to the struc-
ture 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 neu-
tral 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 sim-
ulations, 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 sec-
ondary 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 spec-
tral 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.

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 STORMGE-
NESIS #ANR-22-CE31-0013-01, and the European Research Council project
ExoMagnets (grant agreement no. 101125367).

References

Carolan, S., Vidotto, A., Hazra, G., Villarreal D'Angelo, C., & Kubyshkina, D. 2021, MNRAS, 508, 6001

Chadney, J., Koskinen, T., Galand, M., Unruh, Y., & Sanz-Forcada, J. 2017, A&A, 608, A75

Christie, D., Arras, P., & Li, Z.-Y. 2016, ApJ, 820, 3

Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174

Daley-Yates, S., & Stevens, I. R. 2019, MNRAS, 483, 2600

Erkaev, N., Odert, P., Lammer, H., et al. 2017, MNRAS, 470, 4330

García Muñoz, A. 2023, Icarus, 392, 115373

García Muñoz, A. 2007, Planet. Space Sci., 55, 1426

García Muñoz, A., & Schneider, P. C. 2019, ApJ, 884, L43

García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36

Gillet, A., García Muñoz, A., & Strugarek, A. 2023, A&A, 680, A33

Gronoff, G., Arras, P., Baraka, S., et al. 2020, J. Geophys. Res. Space Phys., 125, e2019JA027639

Guo, J., & Ben-Jaffel, L. 2016, ApJ, 818, 107

Hawley, S. L., Davenport, J. R., Kowalski, A. F., et al. 2014, ApJ, 797, 121

Hazra, G., Vidotto, A. A., Carolan, S., Villarreal D'Angelo, C., & Manchester, W. 2022, MNRAS, 509, 5858

Holmström, M., Ekenbäck, A., Selsis, F., et al. 2008, Nature, 451, 970

Khodachenko, M., Shaikhislamov, I., Lammer, H., et al. 2021, MNRAS, 507, 3626

Koskinen, T. T., Lavvas, P., Huang, C., et al. 2022, ApJ, 929, 52

Kramida, A., Ralchenko, Yu., Reader, J. & NIST ASD Team 2018, NIST Atomic Spectra Database (ver. 5.5.6), (Gaithersburg, MD: NIST)

Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, ApJ, 866, L18

Kulow, J. R., France, K., Linsky, J., & Loyd, R. P. 2014, ApJ, 786, 132

Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121

Lazovik, Y. A. 2023, MNRAS, 520, 3749

Lecavelier Des Etangs, A., Bourrier, V., Wheatley, P., et al. 2012, A&A, 543, L4

Matsakos, T., Uribe, A., & Königt, A. 2015, A&A, 578, A6

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Owen, J. E., Murray-Clay, R. A., Schreyer, E., et al. 2023, MNRAS, 518, 4357

Parker, E. N. 1958, ApJ, 128, 664

Rumenskikh, M., Shaikhislamov, I., Khodachenko, M., et al. 2022, ApJ, 927, 238

Shaikhislamov, I., Khodachenko, M., Lammer, H., et al. 2021, MNRAS, 500, 1404

Strugarek, A., Brun, A. S., Matt, S. P., & Réville, V. 2015, ApJ, 815, 111

Toro, E. F. 2009, in Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin: Springer), 315

Tremblin, P., & Chiang, E. 2013, MNRAS, 428, 2565

Vidal-Madjar, A., Des Etangs, A. L., Désert, J.-M., et al. 2003, Nature, 422, 143

Vidotto, A., & Cleary, A. 2020, MNRAS, 494, 2417

Wood, B. E., Müller, H.-R., Redfield, S., et al. 2021, ApJ, 915, 37

Yelle, R. V. 2004, Icarus, 170, 167