Self-consistent simulation of photoelectrons in exoplanet winds: Faster ionisation and weaker mass-loss rates

Context. Close-in exoplanets undergo extreme irradiation levels leading to hydrodynamic atmospheric escape and the formation of planetary winds. The planetary mass-loss is governed by several physical mechanisms, including photoionisation that may impact the evolution of the atmosphere. The stellar radiation energy deposited as heat strongly depends on the energy of the primary electrons following photoionisation and on the local fractional ionisation. All these factors affect the model-estimated atmospheric mass-loss rates and other characteristics of the outflow in ways that have not been clearly elucidated. Moreover, the shape of the XUV stellar spectra strongly influences the photoionisation and heating deposition in the atmosphere. Substantial changes are to be expected in the planetary mass-loss rate. Aims. We study the effect of secondary ionisation by photoelectrons on the ionisation and heating of the gas for different planet-star systems. We focus on the local and planet-wise effects, to clearly demonstrate the significance of these interactions. Methods. Using the PLUTO code, we performed 1D hydrodynamics simulations for a variety of planets and stellar types. We included planets in the size range from Neptune to Jupiter, and stars from M dwarfs to Sun-like. Results. Our results indicate a significant decrease in the planetary mass-loss rate for all planetary systems when secondary ionisation is taken into account. The mass-loss rate is found to decrease by 43% for the more massive exoplanets to 54% for the less massive exoplanets orbiting Sun-like stars, and up to 52% for a Jovian-like planet orbiting an M-type star. Our results also indicate much faster ionisation of the atmosphere due to photoelectrons. Conclusions. We built a self-consistent model including secondary ionisation by photoelectrons to evaluate its impact on mass-loss rates. We find that photoelectrons affect the mass-loss rates by factors that are potentially important for planetary evolution theories. We also find that enhanced ionisation occurs at altitudes that are often probed with specific atomic lines in transmission spectroscopy. Future modelling of these processes should include the role of photoelectrons. For this purpose, we make a simple but accurate parameterisation for atomic hydrogen atmospheres available.


Introduction
It has long been proposed that the atmospheres of exoplanets orbiting close-in to their host stars must be rapidly escaping (Lammer et al. 2003;Baraffe et al. 2004).Accurately predicting the mass loss rates remains however a difficult task due to the complexity of the star-planet interactions that participate in the mass loss process.
Planets can lose their atmospheres through multiple mechanisms, broadly separated into thermal and non-thermal (see Gronoff et al. 2020, for a complete review on escape processes).Non-thermal processes such as ion sputtering, polar outflow or charge-exchange with the stellar wind dominate for planets orbiting at large distances from their star.In the context of planets on short-period orbits, thermal processes, which include Jeans and hydrodynamic escape, generally dominate.We will focus here on hydrodynamic escape, which occurs when intense stellar radiation deposits its energy in the planetary atmosphere (at various altitudes depending on the radiation wavelength), leading to the escape of the atmosphere into space.
The planetary wind interacts with the stellar wind and forms a number of hydrodynamic features such as a bow shock, a comet-like tail and Kelvin-Helmholtz instabilities (Tremblin & Chiang 2013).Insight into these processes can be gained by measuring with in-transit spectroscopy the excess absorption that occurs at a number of atomic lines such as Hi Lyα and Hα, or the He i triplet at 10830 Å.These features have been measured for hot Jupiters (e.g.Vidal-Madjar et al. 2003;Lecavelier Des Etangs et al. 2010;Jensen et al. 2012) and hot Neptunes (e.g.Kulow et al. 2014;Ehrenreich et al. 2015;Ben-Jaffel et al. 2022).A number of studies have been conducted to understand the global problem of planetary escape with hydrodynamical 1D models (e.g.Yelle 2004;García Muñoz 2007;Murray-Clay et al. 2009;Koskinen et al. 2013).Such models predict velocities for the planetary wind of a few km/s in the planet's proximity, which are consistent with the line broadening of Hα and the He i triplet measured with high-resolution spectroscopy (e.g.Salz et al. 2018;García Muñoz & Schneider 2019).However, they fail to predict the velocities of ∼100 km/s that are found with Lyα spectroscopy (e.g.Vidal-Madjar et al. 2003;Ben-Jaffel 2007).Whether the velocities measured with Lyα spectroscopy are representative of the planetary wind or are instead indicative of stellar wind protons that become neutralised by charge-exchange with the planetary wind remains an open question (Tremblin & Chiang 2013).
Many authors have looked into this problem with 2D and 3D models that include self-consistent radiative transfer, with the goal of constraining the mass loss rates from the planets (e.g.Tripathi et al. 2015;Debrecht et al. 2019).Recently, Shaikhislamov et al. (2021) developed a global 3D multi-fluid model to investigate the He 10830 Å line; Daley-Yates & Stevens (2019) explored the star-planet interaction of the wind in the presence of a magnetic field; and Carolan et al. (2021) looked at the effects of the stellar wind strength on a magnetised planetary wind, showing that the increased polar loss compensates for the decreased mass loss at the dead zones of the planet.
Accurately predicting the atmospheric mass loss rate is essential to answering fundamental questions about the evolution of planets.Indeed, some atmospheres, in particular those of lowmass planets, may be escaping so efficiently that they can be completely lost to space on timescales shorter than the planet's lifetime.In the case of hydrodynamic escape of strongly irradiated exoplanets, the aim of this work, calculating the mass loss rate involves determining how much of the stellar radiation energy is converted into heat, and how much goes otherwise to ionisation and excitation of the atmospheric gas.
We focus here on atmospheric gas made of hydrogen atoms H because although H 2 should be prevalent in primary atmospheres, the molecule will typically dissociate rapidly in the planet's upper atmosphere.Also, the physics of photoelectron interactions with H 2 is more complex than for H atoms because, being a molecule, H 2 offers additional channels for excitation, dissociation and ionisation that must be tracked individually (Hallett et al. 2005).We postpone such a study to future work and focus here on the essentials for an atomic atmosphere.Upon absorption by the atmospheric atomic gas, the stellar Xray and Extreme Ultraviolet photons (jointly referred to as XUV, and covering wavelengths from a few Å to the Lyman continuum threshold at 912 Å) release high velocity electrons.These so-called photoelectrons have energies E 0 =hc/λ−13.6eV, where h and c are Planck's constant and the speed of light, respectively, and λ is the photon wavelength.The photoelectrons can excite the gas and produce secondary electrons while slowing down.The gas heating corresponds to the fraction of the photoelectron's initial energy E 0 that is not expended in excitation or ionisation but that goes instead into kinetic energy.
We emphasize the importance of two parameters that dictate the fraction of the photoelectron's initial energy that is deposited as heat, namely: the local fractional ionisation x e and the energy E 0 of the primary photoelectrons.In particular, if the local fractional ionisation is high, elastic collisions between the fast and thermal electrons ensure that most of E 0 is transferred to the thermal electrons, thereby heating the gas.Also, if E 0 is less than 10.2 eV, i.e. below the lowermost threshold for inelastic collisions in the H atom, all of E 0 is transferred as heat to the background gas in elastic collisions.In these two limits, the prediction of the gas heating is relatively straightforward.Importantly, when the fractional ionisation is low or moderate and E 0 is sufficiently large, additional electrons and ions are created during the secondary ionisation process, when the primary photoelectron interacts with other hydrogen atoms.For the most energetic primary electrons created after photoionisation, it follows a cascade of multiple ionisation events.
The fundamental information for the treatment of photoelectrons, and for the assessment of how much of their energy goes into excitation, ionisation and heating of the gas has been available for decades.For example, Habing & Goldsmith (1971) and Shull (1979) have looked at the production of secondary elec-trons induced by soft X-rays and how they interact with a primordial gas using a Monte Carlo method.The outcome of such modelling efforts can be readily parameterised as a function of x e and E 0 to take into account the effect of photoelectrons in the net ionisation rate and heating rate induced by XUV photons impacting the upper atmosphere of short-period orbit exoplanets.
In this work, we revisit the problem by exploring selfconsistently the simultaneous heating and ionisation that occurs as the photoelectrons slow down, taking advantage of recent advances in the characterisation of the XUV spectra of stars.This effect is often neglected in the studies dedicated to modelling hydrodynamic escape of planetary atmospheres (e.g.Lammer et al. 2003;Wu & Lithwick 2013;García Muñoz et al. 2020, 2021).With the ultimate goal of better predicting the loss rate and the main features of planetary winds, in some cases an arbitrary pre-fixed fraction (typically, 2-30%) of the photoelectron energy E 0 is assumed to transfer into heating (Linssen et al. 2022).Our work shows how to partly overcome such simplified treatments and elaborates on the importance of photoelectrons in the strongly-irradiated atmospheres of some exoplanets.Guo & Ben-Jaffel (2016) have covered related ideas, with an emphasis on how the spectral energy distribution of the star affects the mass loss rate and ionisation of strongly-irradiated atmospheres.Following their initial study, in this work we put significant effort to elucidate how the photoelectrons affect the atmosphere both locally and globally, and provide simple yet accurate descriptions for the self-consistent implementation of such processes in the continuity and energy conservation equations of the gas.These prescriptions should be useful to other modelling efforts.In addition, we leverage recent developments in the reconstruction of the XUV spectra for cool stars by observations at X-ray and far ultraviolet wavelengths (France et al. 2016).
To examine the quantitative effect of photoelectrons and assess whether the planet's gravity plays a role, we performed 1D spherical hydrodynamic simulations with the PLUTO code (Mignone et al. 2007).We focused on the escaping atmospheres of four hypothesised planetary systems with masses ranging from 0.02 M J to 0.69 M J (similar to HD209458b) and assessed the effect of different XUV spectra taken from the MUSCLES survey (France et al. 2016).We defined the planet size R p so that the planet density remains constant in all cases.This is a somewhat arbitrary choice, but useful because of its simplicity.In the energy-limited limit, the mass loss rate (Erkaev et al. 2007) is proportional to the inverse of the planet density, in which case assuming a constant bulk density serves as a valid basis to compare against that limit.Our simulations incorporate the relevant physics of photoelectrons in the continuity and energy equations.
The plan of the paper is the following: the model, including the PLUTO setup and our treatment of radiative transfer with photoelectrons, is described in section §2.In section §3, we describe the relevant physics of photoelectrons and its implementation in PLUTO.Section §4 is dedicated to the description of the results of the atmospheric escape of a Neptune-like planet.Section §5 describes the dependence of our results on the planet mass, and finally in section §6 we study the impact of the shape of the stellar spectra of M and K type stars on secondary ionisation.

Physical model
The model is constructed with the hydrodynamics code PLUTO (Mignone et al. 2007), which solves the Euler equations in a rotating reference frame.The 1D equations for conservation of mass, momentum and energy solved in PLUTO are: where r is distance as measured from the planet's centre and t is time, and where u is the fluid velocity, ρ the density, P T the thermal pressure, E = P T /(γ − 1) + ρu 2 /2 the total energy with γ = 5/3, the adiabatic index for a mono-atomic gas, and C and H are the cooling and heating terms to be defined in §2.3.The joint gravitational potential of the planet plus the star is defined as: where G is the gravitational constant, and M ⋆ and M p are the mass of the star and the planet, respectively.Finally, the centrifugal force is written as: where R orbit is the distance from the planet centre to the stellar centre.e r is the unit vector in the outwards radial direction.
We consider here an atmosphere composed of atomic hydrogen in neutral Hi and ionised H + forms, plus thermal electrons.The density in the equations above corresponds to the total mass density ρ = ρ Hi + ρ H + , which neglects the contribution of thermal electrons.We define the neutral fraction x Hi = ρ Hi /ρ and the fractional ionisation as x e = 1 − x Hi , both taking values between 0 and 1.On the basis of charge neutrality, the number density of ions and electrons is the same.To determine the partitioning between Hi, H + and electrons, we consider the processes for: collisional ionisation: H + e - H + + 2 e -, radiative recombination: H + + e - H + hν , photoionisation: H + hν H + + e -.In our hydrodynamic simulations, we have adapted the Simplified Non-Equilibrium Cooling (SNeq) module from PLUTO (Teşileanu et al. 2008) in the optically thin limit to track the fraction of neutrals and ions of atomic hydrogen in the chemical reaction network equation.In that module, PLUTO solves the equation for the neutral fraction evolution: where c r = 2.6 • 10 −11 × T −0.5 and c i = 5.83 • 10 −11 √ T exp(−157890/T ) are the recombination and ionisation rate coefficients in cm 3 /s, both dependent on the temperature T [K] of the gas, and on the electron number density n e = (ρ/m p )x e , where m p is the proton mass [g/particle].In this work, we added the extra term in Eq. 6 to account for the photoionisation of the gas by the stellar XUV photons impacting the planetary atmosphere.The formulation of the photoionisation rate coefficient J [s −1 ] is given in §2.3.
Our model considers a single thermal temperature to describe the kinetic energy of the neutral, ion and electron components of the gas.This tacitly assumes that the transfer of kinetic energy between them occurs very rapidly, which should be a valid approximation for the relevant fractional ionisations.

Radiative transfer
Photoionisation is the only heating source of the hydrogen gas included in our model.It is caused by XUV photons coming from the star (see fig 1).In principle, the XUV stellar spectrum should consider all possible wavelengths below 912 Å, which is the threshold at which ground-state H atoms can absorb.In practice, at sufficiently short wavelengths the atmosphere of a planet becomes thin, and these wavelengths can safely be neglected in the radiative transfer.In this work, we consider a solar spectrum which starts at 15 Å (see below), and we take that limit also as the shortest-wavelength in our implementation of other stellar fluxes and the H atom cross-sections.In our simulations, we compute the local stellar flux F ⋆ (r, λ) from a reference stellar flux spectrum at 1 AU, which is scaled by (1AU/R orbit ) 2 to produce the top-of-the-atmosphere spectrum F ⋆ (r = T OA, λ), and is attenuated with an optical depth τ λ using Beer-Lambert's law.Namely: The wavelength-dependent optical depth τ λ (r) in the direction towards the star is calculated through: where the integral is performed along the planet-star line of sight.Here, n Hi = (ρ/m p )x Hi is the neutral number density [cm −3 ] and σ λ [cm 2 ] the wavelength-dependent photoionisation cross-section of hydrogen, given by: where σ 0 = 6.3×10 −18 cm 2 is the cross-section at the threshold wavelength λ 0 = 912 Å.For the stellar irradiation, we adopted a solar spectrum downloaded from the SOLID 1 project, as observed on December 13th 2021.We degraded the spectrum in 20 bins of equal size in the range from λ min =15 Å to λ 0 =912 Å.We verified that the choice of 20 bins does not impact the results presented in what follows by performing a few simulations with 40 equal-size bins, which did not show any significant difference.Figure 2 displays the original solar spectrum in black along with the degraded spectrum, the latter represented by the blue bars.The XUV-integrated flux at 0.045 AU is 2174 erg/cm 2 /s.

Cooling, heating and ionisation rates
The cooling term C [erg s −1 cm −3 ] in Eq.3 represents the sum of the thermal energy of the electrons lost by collisional ionisation and radiative recombination: Here, n H = n Hi + n H + is the total number of hydrogen nuclei, and c r and c i are the rate coefficients from Eq 6.Our model does not include molecular chemistry and therefore it neglects for example cooling by H + 3 emission in the infra-red.In this expression, 2.17 • 10 −11 erg is equivalent to the energy of 13.6 eV required to ionise the neutral atom.Correspondingly, 1.07 • 10 −12 (T/11590) erg represents the ∼0.67 kT that are extracted from the kinetic energy of the thermal electrons and lost to radiation if the environment is optically thin.In future work, we intend to move away from the optically thin approximation and treat the recombination into the different bound levels separately, as well as keeping track of the radiated energy.
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: 1 European comprehensive solar irradiance data exploitation; https://projects.pmodwrc.ch/solid where F ⋆ [erg cm −2 s −1 Å −1 ] is the attenuated stellar flux described in Eq. 7. Φ λ,xe is the number of secondary ions created per photoionisation, to be described in section 3. The heating deposition rate H [erg s −1 cm −3 ] in the energy equation can be expressed as: with η λ,x e being the heating efficiency to be described in section 3. Omitting the negative term in the parenthesis of the above equation, and assuming η λ,x e = 1, would provide the incident energy that is deposited in the gas.The negative term in the parenthesis subtracts from it the energy that is expended to produce the primary photoelectrons during photoionisation.Models that adopt a number of secondary ions Φ λ,x e = 0 and a heating efficiency η λ,xe = 1 assume that the surplus energy of the photoelectrons after photoionisation goes entirely into heating the gas.However, this is not generally true, since the photoelectrons also cause excitation and ionisation of the surrounding neutral atoms.Ideally, one has to determine precisely the heating efficiency by estimating which fraction of the energy goes into heating and which fraction is being used for the other processes.
Here, we aim to properly calculate the effect of photoelectrons on the chemistry and heating of the atmospheric gas, taking into account the dependence of these effects on the photoelectron energy and fractional ionisation.In section 3 we detail the physical processes in play and how to determine and parameterise η λ,xe and Φ λ,xe .

Numerical Method and initial set-up
To solve Eqs.1-3 we use the Harten-Lax-Van Leer approximate Riemann Solver that solves exactly stationary contact discontinuities between the cells (HLLC, see Toro 2009).We use linear reconstruction for the spatial order of integration and a third order Runge-Kutta scheme (RK3) for the time evolution.According to the stiffness of equations 3 (including cooling 10 and heating 12) and 6, a dynamically-adaptive integration strategy is adopted (see Teşileanu et al. 2008, for more details).
PLUTO uses dimensionless units rather than for example cgs units so that flow quantities can be properly scaled to "reasonable" numbers (Mignone et al. 2007).We work with three fundamental units : ρ 0 , L 0 and V 0 .The reference density at the base of the atmosphere is fixed for all planetary systems at ρ 0 = 1.326 × 10 −10 g/cm 3 corresponding to a pressure of 12 µbar (value taken from most models to define their boundary condition for a temperature of 1100 K.The reference length L 0 is equal to the planetary radius R p and the reference velocity V 0 is calculated as the following V 0 = (GM p /L 0 ) 0.5 .Their values are given in Table 1 for the models considered in this study.The computational domain is composed of a stretched grid of 500 cells between 1 and 30 R p with a stretch factor of 1.017.The grid is extended below 1 R p by adding 10 cells of uniform size between 0.999 and 1 R p .This extension is used by PLUTO to establish the boundary conditions at the bottom of our model.
In our simulations, the atmosphere is initialised with the density profile:  Notes.R orbit is the orbital distance of the planet calculated as R orbit = (0.045AU/L 0 ), M ⋆ /M p , L 0 and V 0 are the orbital distance, star-to-planet mass ratio, reference length and velocity, respectively.
where α ρ sets the initial density scale height in the atmosphere.
For the simulations presented in the work, α ρ is set to 20.The atmosphere is initialised everywhere at the equilibrium temperature T eq , which is equal to 1100 K in all our models in which the planet orbits a Sun-like star.The ideal gas law is used to initialise the pressure, and the neutral hydrogen fraction profile is initially set as: and the velocity is set to zero.
In sections 4 and 5, all planetary systems are considered to orbit a solar-type star with a separation of R orbit = 0.045 AU.Table 1 lists the physical parameters used for all simulations of each planetary systems.We choose four different planetary masses : 0.02 0.05, 0.1 and 0.69 M J .This ranges from a sub-Neptune planet to a Jupiter-like planet similar to HD209458b.In section 6, we consider different orbital distances for the planets around different stars to ensure having the same integrated flux.
At the inner boundary, we set the density and the pressure (12 µbar) to their initial value, the velocity to zero, and assume that the gas is entirely neutral.At the outer boundary, we use a free outflow boundary condition, with the gradients of all variables (P, ρ, v r and x HI ) equal to 0.

Theoretical baseline
Before thermalising and depositing their kinetic energy into heat, photoelectrons interact with the neutral hydrogen atoms in the gas multiple times, losing at each collision a part of their energy.This interaction has been investigated by a number of authors for a variety of astrophysical applications (Habing & Goldsmith 1971;Shull 1979;Furlanetto & Stoever 2010).In the investigation of exoplanetary atmospheres, Cecchi-Pestellini et al. ( 2006) and Shematovich et al. (2014) studied the effect of stellar X-ray irradiation on the heating by including the effect of photoelectrons and found that the heating efficiency is notably less than 1 if the atmosphere is mostly neutral when photoelectrons are included.Their simulations are not self-consistent, though, as those studies do not let the modified hydrodynamic solution alter the conditions that the photoelectrons experience.Guo & Ben-Jaffel (2016) investigated the production of secondary ions in H 2 /H atmospheres, and found an increase of the total ionisation rate by a factor of 10 in the region r < 1.05 R p that is reached only by the shortest-wavelength photons.They also report a drop in the mass loss rate by less than a factor of 2 when photoelectrons are included, although the connection with the local heating efficiency is not established.García Muñoz (2023) has looked into the photoelectron-driven processes that affect the population of the first excited level of hydrogen that is sensed in transmission spectroscopy of the Hα line.
Based on the physics of the Hi atom excitation and ionisation, several ranges of energies may be considered.For E 0 <10.2 eV, the threshold for the first excitation level, all the energy of photoelectrons is deposited as heat (Dalgarno & McCray 1972;Cravens et al. 1975).For energies 10.2 eV <E 0 <13.6 eV, the surplus of energy goes either into heating and into excitation.Finally, for energies E 0 > 13.6eV, which defines the ionisation threshold, it can additionally be used to ionise further the gas (see Fig. 3).The energies 10.2 eV and 13.6 eV are the thresholds for E 0 that enable excitation and ionisation of the atom after photoionisation, corresponding to the wavelengths 520 and 455 Å respectively.For even shorter-wavelength photons, the ejected photoelectron will be able to excite/ionise the hydrogen atoms multiple times.It is reasonable to assume that the excited Hi (resulting from either collisions of ground state H with photoelectrons or with thermal electrons or from the recombination of protons) atom will radiate away its excitation energy.Indeed, al-though the Lyman-α line (and other lines in the Lyman series) is not optically thin, it occurs that even after accounting for opacity its effective radiative timescale is shorter than the collisional deexcitation timescale at the relevant atmospheric pressures (see García Muñoz 2023).
We calculated the heating efficiency and ionisation yield with the model described in García Muñoz (2023).The calculations were conducted over a grid in energy E 0 that resolves well the details in the heating efficiency and the ionisation yield at five values of the fractional ionisation, namely, x e =10 −4 , 10 −3 , 10 −2 , 10 −1 and 1. Figure 4 shows the heating E h as a function of the primary electron energy E 0 and x e (different colours).The black solid line for x e = 1 coincides with E h =E 0 .We define the heating efficiency, η λ,x e = E h /E 0 which corresponds to the ratio between the energy that goes into heating E h and the initial energy of the photoelectron E 0 .As noted earlier, for any prescribed x e , the heating efficiency depends on E 0 or, equivalently, on the energy of the incident photon hc/λ, and takes values from 0 to 1.
When η λ,xe = 1 all energy is deposited as heat while when η λ,xe < 1 a fraction of the energy is diverted to excite or ionise other atoms.To facilitate its use in our hydrodynamical model and in other similar models, we fitted the MC calculations of η λ,xe to 4 th order polynomials: The polynomial coefficients a n,λ can be found in Appendix A. When the fractional ionisation falls below 10 −4 , η λ,x e becomes essentially independent of x e .In those conditions, we fix x e to 10 −4 to estimate η λ,x e .x e =1 x e =1.e-4 x e =1.e-3 x e =1.e-2 x e =1.e-1 Fig. 4. Energy E h deposited as heat by primary electrons of energy E 0 .Solid curves represent different fractional ionisations.The heating efficiency is given by η λ,xe = E h /E 0 .Adapted from García Muñoz (2023).
During the secondary ionisation process, a number Φ λ,xe of secondary ions and electrons are produced from the primary electrons of energy E 0 .Figure 5 reports Monte Carlo model calculations of the number of secondary ions created for a primary electron of energy E 0 , and for a range of conditions of fractional ionisation x e .For large photoelectron energies, numerous secondary ions can be created.For instance, for E 0 = 300 eV (λ = 39.5 Å) up to 5 secondary ions can be created if x e = 10 −1 (green curve) and up to 10 if x e = 10 −2 (orange curve).Simon Wedlund et al. (2011) have shown that about 8.5 secondary ions are expected to be created at E 0 = 300 eV in a non-ionised atmosphere.Here, we consider equally-spaced bins to sample the wavelength of the impacting photons, which translates into a first wavelength bin covering 15 Å (826.5 eV) to 59.85 Å (207.15 eV) (see Fig. 2).In this bin, we predict the average production of about 14 secondary ions, which is significantly larger because the considered bin encompasses a range of E 0 that is on average larger than 300 eV.When considering smaller wavelength bins, we recover the production of 8.5 secondary ions at E 0 = 300 eV predicted by Simon Wedlund et al. (2011).We reiterate nonetheless that our results are not significantly affected when considering smaller wavelength bins.
Conversely, for low photoelectron energies, no additional ions are created (as expected, see Fig. 3).
Similarly, we performed a 4th order polynomial fit: with b n,λ the polynomial coefficients found in Appendix A. x e =1 Fig. 5. Secondary ions Φ λ,xe produced by primary electron of energy E 0 .Solid curves represent different fractional ionisations.We note that these yields refer to energy bins of finite size rather than to the specific energies quoted there.This distinction makes some difference at the higher energies because the energy bins are larger for them.We have nevertheless confirmed by running a few tests with a refined spectral grid of 40 bins that this has no significant impact on the overall gas properties and mass loss rates.Adapted from García Muñoz (2023).
Because for sufficiently large energies E 0 , η λ,xe is small and Φ λ,xe is large when x e is small, significant effects on the heating and ionisation of the atmosphere are expected in the highpressure region where the gas remains neutral and that only very energetic photons can reach.Conversely, we expect that further out where the local fractional ionisation is high and the newly created photoelectrons have small energies, the effect of photoelectrons will be weak.
A note is due on the terminology used throughout this paper.We will say that a calculation incorporates secondary ionisation when it is done with the parameterised forms of η λ,x e and Φ λ,x e described above, which obviously includes the effect of both excitation and ionisation by the primary and secondary electrons.
Otherwise, we will say that the calculation does not incorporate secondary ionisation when it is assumed that η λ,x e =1 and Φ λ,x e =0.

Photoionisation and heating rates
A main goal of this paper is to investigate the feedbacks between photoelectron-driven processes and the hydrodynamic outflow, which we detail in section 4. Before that, and to demonstrate the effects of photoelectrons, we will illustrate how they affect the energy deposition and ionisation for a prescribed atmospheric profile.For this demonstration, we assume the atmospheric profile to be fixed, meaning that we do not consider the feedbacks driven by the photoelectron deposition in the gas.Through this exercise, we want in particular to assess the intricate link between the shape of the solar spectrum (fig. 1) and the deposition of the associated energy into the atmosphere.The shape of the stellar XUV spectrum partly dictates where in the atmosphere and at which wavelengths the photoionisation and heat deposition are more prominent.
For this exercise, we consider a reference profile taken from one of our numerical simulations: the steady-state atmosphere of a Jupiter-like planet with a mass of 0.69 M J (M0.69 in Table 1).The top panel of Figure 6 shows the fractional ionisation x e as a function of the neutral density n Hi .In addition, we show in the bottom panels the wavelength contributions of photoionisation (left) and heating (right) at the five representative altitudes depicted in the top panel.
Each left panel of Fig. 6 shows two curves.The solid curve corresponds to J ′ x Hi = σ λ F ⋆ λ hc x Hi and shows the contribution to ionisation of each wavelength bin, excluding the effect of secondary ionisation by the photoelectrons.In turn, the dotted curve corresponds to J ′ 1 + Φ λ,xe x Hi and shows the total ionisation, including the production of secondary ions.The integral of the latter over wavelengths corresponds to the net ionisation rate J x Hi of Eq. 6.
When the unattenuated stellar flux hits the upper atmosphere of the planet (see panel for 3.71 R p , in magenta), only the longest wavelengths contribute to photoionisation, producing a large peak observed at the threshold of 912 Å.This is largely because of the strong modulating effect of the photoionisation cross section, that prioritises the longer wavelengths.As we go deeper into the atmosphere (1.58 R p , green), the number density increases by two orders of magnitude, the atmosphere becomes optically thicker, especially at the longer wavelengths and, a broader range of wavelengths contribute to photoionisation with the presence of multiple peaks (at 300 Å, 600 Å and 912 Å).At 1.29 R p (in red), photoionisation is dominated by a peak around 300 Å.The atmosphere is optically thick at the longest wavelengths, and only the ones lower than 400 Å will reach this altitude, as the photoionisation cross-sections are smaller at short wavelengths.For the same reason, only the most energetic photons penetrate below 1.05 R p (in blue).Photoionisation occurs at a slower rate at this altitude than in the upper atmosphere, yet the number of electrons produced per volume unit is much larger, because it partly follows the hydrogen density.
Interestingly, much faster ionisation occurs when photoelectrons are taken into account (dotted lines).The enhancement of ionisation rate by photoelectrons strengthens at the lowest altitudes.This is clearly seen at 1.05 R p (in blue) and 1.02 R p (in cyan) due to the multiplicative term (1 +Φ λ,x e ) in the integral of J x Hi .In these deep layers, the photoelectrons enhance the double peak structure between 15 and 200 Å, and ionisation is increased by a factor of two.The existence of multiple peaks results from the shape of the stellar spectrum at short wavelengths, and the extra modulating effect of Φ λ,x e at very short wavelengths.
Similarly, the right panels of Fig. 6 elaborate on the energy term of Eq. 12, and we can define three contributions.H ′ = n Hi σ λ F ⋆ depicts the energy available before photonionisation (solid lines); H ′ 1 − λ λ 0 represents how much energy might potentially be effectively deposited as heat after photoionisation (dashed lines); lastly H ′ 1 − λ λ 0 η λ,x e shows the total heating rate per wavelength bin, taking into account the effect of photoelectrons (dotted lines).
The energy deposition at different altitudes in the atmosphere is depicted in the right panels of fig 6.In the upper atmosphere, for the highest altitudes, 3.71 and 1.58 R p , the energy available before (solid line) and after (dashed line) photoionisation originates from a wide range of wavelengths.However, the heat deposition rate is 100 times larger at 1.58 R p than at 3.70 R p , as a result of a higher density.For low-energy photons, the term 1 − λ λ 0 is close to 0. Those photons spend most of their energy in the primary ionisation event, and very limited energy remains for heating (or for further ionisation or excitation).As altitude decreases, the peak of energy deposition shifts towards shorter wavelengths, and most of the net energy extracted (dashed line) from the available energy (solid lines) originates from λ < 400 Å.The peak observed at the photoionisation threshold of 912 Å at the two highest altitudes disappears because most of the energy available at these wavelengths has been used for primary ionisation at higher altitudes.For the region the closest to the surface (in cyan), we observe that the energy available before H ′ and after H ′ 1 − λ λ 0 photoionisation overlap.Indeed, for wavelengths lower than 200 Å, the term 1 − λ λ 0 is close to 1 and E 0 is marginally impacted after removing 13.6 eV.
An additional modulation occurs when secondary ionisation is considered, described by the dotted lines in the right panels.We observe a significant attenuation of the peaks for altitudes lower than 1.29 R p at wavelengths lower than 400 Å.For very short wavelengths, η λ,x e is rather small (∼0.12) and acts as a significant modulation factor, reducing substantially the heating rate arising from those wavelengths.In other words, the contribution of X-rays towards heating is very weak at all altitudes in our investigation when secondary ionisation is taken into account.
In summary, we observe that taking into account photoelectrons leads to a significantly higher ionisation rate deep in the atmosphere at short wavelengths.Concomitantly, this increase implies that less energy is deposited into heat at these low altitudes, which will naturally change the steady state of the planetary atmosphere.We now turn to the investigation of possible feedbacks in the atmosphere, which we approach through numerical simulations.

Characteristics of the escaping atmosphere
We now focus on the self-consistent simulations based on the PLUTO code, which include the hydrodynamics and photochemistry of the gas and, optionally, a parameterised description of secondary ionisation.The implementation of this new model was validated by comparing our results in the limit η λ,x e =1 and Φ λ,x e =0 of no secondary ionisation with a version of the model (solid lines) and with photoelectrons (dotted lines), as a function of wavelength.The difference between the solid and dotted lines at lower altitudes clearly demonstrates the importance of multiple ionisations driven by the photoelectrons.bottom right panels: Contribution to the heating rate at different altitudes, H ′ [erg s −1 cm −3 ] representing the energy available before photoionisation (solid lines), the dashed line representing the heating rate due to photoionisation without the effect of photoelectrons, and the dotted lines the net heating rate taking into account secondary ionisation.Dashed and dotted lines overlap at the higher altitudes.The curves provide valuable insight into the wavelengths that contribute to heating at each altitude and how that energy is affected by photoelectrons.described in García Muñoz (2007) adapted to the physical conditions of the current simulations.The comparison between both sets of simulations proved excellent, with errors in the main properties (mass loss rates, temperature, densities, velocities) of less than a few percent.As a novel feature, we added to our current model the self-consistent effect of photoelectrons, as described in sections 2-3.We focus in this section on a Neptunemass planet (case M0.05 in Table 1).Figure 7 shows the temperature, density, velocity, Mach number, pressure, neutral fraction profiles, along with the ionisation and heating deposition rates as a function of distance to the planet.
We will first describe the case without secondary ionisation (black curves, Φ λ,xe = 0, η λ,xe = 1).The profiles show the usual features that have been reported by other models of close-in exoplanets.In particular, the density (a) and pressure (c) drop rapidly from the bottom to the top of the atmosphere, while the fractional ionisation increases monotonically (b).At 5 R p the gas remains mostly neutral with x e = 0.2.As the flow expands, it accelerates (d), slowly at low altitudes and faster at high altitudes, and reaches the sonic point at about 2.2 R p , marked by a black dot on the Mach profile (f).Above this altitude, the planetary wind becomes supersonic, reaching velocities of about 13 km/s at 5 R p .The temperature profile (g) shows a more complex behaviour.Very close to the bottom of the model, the temperature decreases with altitude because little XUV radiation reaches those layers.The availability of unattenuated XUV radiation causes the rapid increase in temperature (g), ionisation (e) and heating rate (h) that occurs immediately above.Indeed, as we move high in altitude, the temperature rises up to about 7500 K at 8.5 R p .Higher up, the temperature drops due to the expansion of the atmosphere and a lack of energy deposition as the gas becomes fully ionised and optically thin.
When secondary ionisation is taken into account (showed by the red lines on Fig 7), the density (a) and pressure (c) drop faster with altitude, whereas the fractional ionisation (b) increases significantly at all altitudes.The flow velocity (d) also increases faster than in the case without secondary ionisation, yet reaches the sonic point at about the same distance of 2.1 R p (marked by a red dot on the Mach profile).The temperature structure of the atmosphere is also affected when secondary ionisation is considered.Indeed, the temperature (g) peaks closer to the planet at 5 R p and at a lower temperature around 5800K.This cooler temperature is directly associated to the fact that photoelectrons extract some of the available energy to excite and ionise further the gas instead of heating it.Higher ionisation rates (e) cause higher fractional ionisations (b).The drop in the heating rates is particularly significant deep in the atmosphere, where the gas remains mostly neutral, and the heating efficiency is very low (see below).
We calculate the mass loss rate as if the atmosphere was escaping symmetrically over all directions: This is just a convenient assumption adopted by some 1D models that can be easily modified if instead one assumes the gas escapes only through for example the dayside, in which case the quoted mass loss rates should be divided by 4. The calculated mass loss rate for our planet M0.05 is Ṁ = 2.89×10 11 g/s without secondary ionisation and Ṁ = 1.43×10 11 g/s with secondary ionisation, decreasing the mass loss rate by about 50%.This comparison shows that the drop in density when secondary ionisation is included dominates over the increased velocities that are achieved.This 50% drop in the mass loss rates can be connected to the heating efficiency of the gas (η λ,xe in section 3) at the altitudes where most of the stellar energy is deposited, a point that is discussed at length in what follows.A drop in the mass loss rate of that magnitude might have significant implications on atmospheric evolution studies and also on the detectability of diagnostic lines such as Hi Lyα and Hα in transmission spectroscopy, ideas that we will explore in future work.

Heating efficiency
The amount of energy available to excite and ionise further the gas depends on the primary energy E 0 and on the fractional ionisation x e , see section 3. To understand better how the stellar energy is finally deposited into the atmosphere, and the wavelengths and altitudes over which this occurs, the top panel of fig.8 shows the heating efficiency η λ,xe for our M0.05 planet as a function of altitude and wavelength.Following fig.4, η λ,x e = 1 (gray regions) for λ > 520 Å at every altitude due to the low energies of the primary photoelectrons, see section 3.1.The effect of photoelectrons is the most significant close to the surface (blue and green regions in fig 8) where the gas is largely neutral, with x e < 10 −3 .η λ,x e decreases steadily from 0.45 at λ = 520 Å to 0.13 for the most energetic photons (λ < 82 Å).These photoelectrons spend from 55 % to 87% of their energies to excite and ionise further the gas.We note that the heating efficiency for primary photoelectrons produced by radiation λ < 82 Å is particularly low and remains less than 0.7 up to 3.5 R p .Higher altitudes are mildly affected by the presence of photoelectrons, and above 5 R p the fractional ionisa-tion x e is very close to one, rendering the effect of photoelectrons completely negligible (grey area).Above 5 R p the fractional ionisation x e increases steadily and the photons ionising the atmosphere are low-energy, high wavelength photons that do not lead to significant secondary ionisation (see section 3).
The efficiency η λ,x e provides insight into the specifics of wavelength-by-wavelength energy deposition in the atmosphere.To gain insight into the local heating at each altitude, even if losing some information encoded in the dependence of η λ,x e on wavelength, we define in addition at each altitude an effective heating efficiency ηr .It is calculated from the ratio of the actual heating deposition rate of Eq. 3 over the net energy endowed to the photoelectrons after photoionisation: The bottom panel of Fig. 8 shows the variation of ηr with altitude in the atmosphere for the M0.05 case.As altitude increases, the effective heating efficiency ηr increases significantly because of the increase of the fractional ionisation and because most of the nascent photoelectrons have low energies.Above 3 R p , ηr is close to 1 and most of the energy is directly converted into heat without any significant effect of secondary ionisation.Below 3 R p , ηr decreases steadily towards the bottom of the atmosphere, decreases sharply, to reach 0.5 at 1.2 R p and to reach 0.13 at the bottom.This later latter region is particularly important to establish the mass loss rate because it is at these altitudes that the XUV photons become deposited.

Secondary ions per primary photoelectron
Similarly to the heating efficiency, we calculate the number of total (primary + secondary) ions created per primary photoionisation event.The top panel of fig 9 shows 1+Φ λ,xe for the M0.05 planet as a function of altitude and wavelength.
As anticipated, for wavelengths longer than 455 Å, the primary photoelectrons will not produce any secondary ions (grey area in the top panel, see section 3.1).Each primary photoelectron is able to create, on average, one secondary ion for wavelengths between 455 and 306 Å close to the planet when the gas is mostly neutral.For larger energies of the primary photoelectrons, the region where secondary ions are created extends toward higher altitude.Correspondingly, below 2 R p , only shortwavelength photons reach those altitudes where the fractional ionisation is low, and photoelectrons will be able to create numerous secondary ions after primary ionisation.Indeed, for the shortest wavelengths below 82 Å, up to 14 ions can be created per primary ionisation, reducing substantially the heating efficiency and boosting the photoionisation rate.
In order to evaluate the average number of ions created per primary energy E 0 , an effective 1+ Φr is defined at each altitude.It is calculated as: The bottom panel of fig 9 shows 1+ Φr for the M0.05 case.The number of ions created by the primary photoelectron is close to one (i.e.no secondary photoelectrons are created) above 1.2 R p , as expected because at those altitudes only low-energy photons interact with the gas and the fractional ionisation can be relatively high.The most striking feature of this panel is that close to the surface, up to 14 ions are created per primary ionisation.
The significantly enhanced ionisation rate in that region is the ultimate cause of the faster overall ionisation of the atmosphere when secondary ionisation is considered.
In summary, it is in the lower layers of the atmosphere, where only photoelectrons of high E 0 are created, that the heating efficiency becomes the lowest and the number of secondary ions created becomes the highest.The gas ionises faster when the effect of photoelectrons is included, which entails that the amount of heat deposited decreases and the cross-section of the planet to XUV photon decreases as well.The diminution of heat deposition naturally translates into a more compact planetary atmosphere and a reduced net mass loss.Moreover, the decrease in the heating rate we observe is not sufficient to make non-thermal processes dominate for planets on such close-in orbits.

Dependency on planetary mass
The survival of a planetary atmosphere depends on the incoming stellar XUV flux and how this energy is used to overcome the planet's gravitational potential.The low gravitational potential of low-mass planets makes their atmospheres more extended than higher mass planets, with a larger pressure scale-height.Depending on how the neutral density and fractional ionisation profiles change with planetary mass, the relative effect of photoelectrons can be more or less pronounced.We report in table 2 the mass loss rate for all cases presented in this study.A reduction of 43% of the mass loss rate is observed for the case M0.69 and up to 54% in the case M0.02 and a dependency of the planetary is shown: we immediately notice that secondary ionisation is slightly more efficient for low mass planets.
To investigate this aspect, we now turn to analyse the two extreme cases in our sample: a planet with M p = 0.69M J (model M0.69) and a planet with M p = 0.02M J (model M0.02) which is at the mass-limit for Roche lobe overflow.Both models are represented on Fig. 10 with M0.69 in blue and M0.02 in green.In these figures, the dotted lines denote the cases with no secondary ionisation, and the dashed lines the case considering secondary ionisation.In panel (a) the fractional ionisation is shown as a function of neutral density.We observe that, for a given neutral density, the fractional ionisation x e is always smaller for the least massive planet (M0.02) when secondary ionisation is not included (dotted lines).As a result, one naturally expects that secondary ionisation should be more effective for the low mass planet, since its lower ionisation fraction (at a given neutral density) implies that more energy will be used to excite and ionise further the atmosphere instead of heating it.When we include secondary ionisation, this trend is indeed confirmed, and we observe a slightly larger shift of the x e -n Hi profile for the low mass planet (M0.05, blue lines) than for the larger mass planet (M0.69, green lines).
In addition, we observe that in the two cases studied here, the peak of the ionisation rate (panel b) shifts to a higher altitude when secondary ionisation is included and the overall ionisation rate increases slightly.We note that these shifts and increases are more pronounced for the low-mass planet (M0.05, in blue).Parallelly, the energy deposition rate in panel (c) is also significantly impacted in both cases, similarly to the case M0.05 (section 4).Close to the surface, when photoelectrons are taken into account (dashed lines), the heating rate is reduced by one order of magnitude for M0.02 and divided by 2 for M0.69.As a result, the effect of secondary ionisation is found to be important for all the planets studied in this work, and to be the most pronounced for the least massive planets we considered.

Atmospheric escape with different stellar spectra
An accurate description of stellar spectra is fundamental for modelling atmospheric escape.However, the stellar output at XUV wavelengths is quite difficult to constrain and uncertainties remain in the stellar spectra that are adopted by models.One of the difficulties arises from the fact that the high-energy radiative output of stars is variable.The main difficulty, however, is that while the X-ray spectrum is usually measurable, the EUV part is strongly attenuated by the interstellar medium (ISM).Much work has been conducted in recent years to overcome some of these difficulties, either by exploring larger sets of stars at Xray and FUV wavelengths or by complementing the available data with elaborate modelling of stellar coronas.We focus here on a planet similar to the case M0.69 orbiting cool, low-mass stars, and explore the impact of their XUV spectra on the photoionisation and deposition rate into the atmospheres of close-in exoplanets.

Stellar spectral energy of K and M type stars
As seen above, the conversion from stellar photon energy to energy powering the atmospheric escape depends, amongst other factors, on the actual energy distribution of the primary photoelectrons.It is therefore expected that the stellar spectral energy distribution (SED) may have an impact on the macroscopic properties of planetary atmospheres.
For this purpose, we choose a sample of three stars taken from the MUSCLES data survey, GJ1214, GJ832 and HD97658 (France et al. 2016) (https://archive.stsci.edu/prepds/muscles/).The aim of this survey was to obtain the SEDs of M and K dwarfs from 5 Å to 5.5 µm, with special emphasis on the difficult-toconstrain XUV.Typically, the X-ray part of the spectrum are measurements from Chandra and XMM-Newton (Loyd et al. 2016), and the EUV part was reconstructed following the techniques described by (Youngblood et al. 2016), who in turn used an empirical scaling relation based on Lyα flux (Linsky et al. 2013).However, the X-ray spectrum for HD97658 does not  come from XMM/Chandra observations but scaled from X-ray data of HD 85512 according to the ratio of bolometric flux between the two stars because of their similar levels of Fe XII emission relative to the bolometric flux.The synthetic EUV spectrum of the M dwarf GJ 832 was taken from Fontenla et al. (2016).We selected the panchromatic spectrum with adaptive resolution for our sample (Panchromatic SED binned to a constant 1 Å resolution, downsampled in low signal-to-noise regions to avoid negative fluxes).Choosing the spectrum without adaptive resolution does not affect our wavelength range of interest.
Similarly to the Sun, we scaled the stellar fluxes taken at 1 AU by (1AU/R orbit ) 2 , see section 2.2, and degraded their spectra into 20 bins of equal sizes.In order to focus on the shape of the SED, we changed the orbital distance of the planet around each star to ensure that all planets receive the same integrated XUV flux (2174 erg/cm 2 /s) as in the preceding sections.We display in the upper panel of Figure 11 the original MUSCLES spectra, and in the lower panel the spectrum used in our model.The corresponding R orbit are given in table 1.
As seen on the bottom panel of in Fig. 11, the scaled flux of both M dwarfs is higher than the flux of the Sun and of HD97658 between 15 Å and 400 Å, and is lower above 400 Å.Because the effect of photoelectrons is more intense at high energies (i.e.short wavelengths), we expect to have a stronger effect of the photoelectrons for the M dwarfs in our sample.

Heating efficiency and ionisation yield
From our hydrodynamic calculations as described in section 3, the heating efficiency η λ,x e for each star-planet model is shown in the top panels of Fig. 12 organised by increasing effect of secondary ionisation from left to right.Similarly to the solar case (third column), the effect of photoelectrons is the most significant close to the surface (blue and green regions in fig 12) where the gas is largely neutral, with x e < 10 −3 .However, the effect of secondary ionisation is more prominent for the M dwarfs (GJ1214 and GJ832, first and second columns).A larger region of the atmosphere is impacted with η λ,x e < 1.Interestingly, the results for HD97658 (last column) show similar heating efficiencies compared to the solar case, despite significant differences in their binned spectra above 400 Å. Figure 13 represents the heating efficiency ηr as function of altitude for all stars.We observe that GJ1214 and GJ832 have lower ηr for a given altitude (in units of R p ). HD976458 displays a similar ηr compared to our Sun.
Along with the heating efficiencies, the bottom panels of figure 12 show the ionisation yield Φ λ,x e .For the M dwarfs, photoelectrons create, on average, at least one secondary ion onto a larger region of the atmosphere (blue areas) compared to the Sun and to HD97658.The effect of secondary ionisation is the most prominent for GJ1214 at the shortest wavelengths, up to 14 secondary ions are created below 3 R p .
We report in table 2 the mass loss rate for all cases.A reduction of 52% of the mass loss rate is found for the case GJ1214 when secondary ionisation is taken into account.We, therefore, confirm that the effect of secondary ionisation by photoelectrons is stronger for the M dwarfs of our sample.

Conclusion
In this paper, we investigated the role of the XUV flux of a host star on planetary atmospheres, with a focus on the physics of photoelectrons.In order to understand how photoionisation  Notes.Spectral type and M * are taken from (France et al. 2016).M p and R p refer to the M0.069 case.We also display the percentage of the XUV flux that goes into the X-ray flux (λ < 100 Å) and the 1D spherical mass loss rates with and without the secondary ionisation by photoelectrons.
and deposition of energy determine the structure of the planetary atmosphere, we have parametrised the excitation, secondary ionisation, and heating due to photoelectrons in a planetary atmosphere.We implemented this parametrisation in the PLUTO code to perform 1D modelling of escaping atmosphere of exoplanets orbiting a solar like star with a mass range between 0.02 and 0.69 M J assuming a constant planetary density.We studied the impact of secondary ionisation by photoelectrons on those processes, and on the net mass loss rate of the atmosphere.We characterised the importance of the fractional ionisation x e and the energy E 0 of the primary photoelectrons on the gas in the planetary atmosphere.We emphasised the importance of the shape of the XUV flux by considering different stellar spectra of K and M type stars from the MUSCLES data survey (France et al. 2016).
We have seen that the shape of the solar spectra in the XUV range affects significantly the photoionisation and the way the energy is deposited in the atmosphere.Photoelectrons lead to higher ionisation rates deep in the atmosphere, where the fractional ionisation of the gas is low.Consequently, the amount of energy that goes into heating is reduced when secondary ionisation is taken into account, and so does the mass loss rate of the atmosphere.Furthermore, we find that the stellar X-ray flux contributes significantly less to the heating of the atmosphere, and therefore contributes significantly less to setting the strength of the atmospheric escape.
We found that the mass loss rate is decreased prominently from 43% to 54% upon consideration of secondary ionisation for planets with masses of 0.69 M J and 0.02 M J .We find that the diminution of the net mass loss rate is slightly more important for planets with lower gravitational potential than the more massive ones.
The shape of the stellar spectrum is a key parameter to understand the escape of planetary atmospheres.In particular, the spectral energy distribution of the XUV part of a solar-type star can vary greatly due to difference in mass, age, metallicity or magnetic activity.We implemented a sample of stellar spectra of K and M-type stars in our model and studied the effect of XUV shape on secondary ionisation and on the mass loss rate.We found that, for a Jupiter-like planet, the effect of secondary ionisation is stronger for the M dwarfs with lower heating efficiencies and higher ionisation yields compare to the Sun.The contribution of this effect is largely due to wavelengths smaller than 400 Å that are comparatively stronger than the solar spectrum for the two M dwarfs studied in this work.The EUV flux between 100 and 400 Å, which contributes the most to the heating, is higher for GJ1214 and GJ832, meaning that those stars would heat up efficiently the atmosphere of their planets.On the other hand, active stars in the X-ray (λ < 100Å) range will not necessarily accentuate the heating, since most of the energy of the photoelectrons will be used to re-ionise the gas.
To conclude, we stress that secondary ionisation is an important physical process and must be taken into account when modelling an escaping planetary atmosphere in 1D.In a future work, we will assess how photoelectrons impact outflow detection when considering more realistic 2D and 3D geometries.In addition, it remains unclear how secondary ionisation affects the atmospheric outflow when Joule heating is considered.Since more ions are created, it may, in some cases, lead to an increase of the Joule heating (Cohen et al. 2014) and impact the mass loss rate.The presence of a magnetic field or the consideration of a more complex chemical network may therefore also modify the mass loss rate, which we also aim to characterise in a subsequent work.

Fig. 1 .
Fig. 1.Illustration of the model of photoionisation and heating of an atmosphere.The atmosphere is irradiated by XUV stellar photons with an incident flux Fxuv.The gas particles are under the influence of stellar gravity F star , planet gravity F planet and centrifugal forces F cent .

Fig. 2 .
Fig. 2. Solar XUV spectrum.Black: Original solar spectrum downloaded from SOLID at 1 AU.Blue bars: Binned spectrum used in our calculations.

Fig. 3 .
Fig.3.Schematic of the energy channelling from an impacting photon into a photoelectron that can lead to the heating of the gas, its excitation, and secondary ionisation.

Fig. 6 .
Fig.6.Expected effect of photoelectrons on ionisation and heating rates.top panel: Fractional ionisation x e as a function of neutral number density n Hi in the case of a Jupiter-like planet irradiated by a solar spectrum.Five particular altitudes are labelled by coloured squares and correspond to the series of panels below.bottom left panels: Contribution to the photoionisation rate at different altitudes, without secondary ionisation J x Hi [s −1 ] (solid lines) and with photoelectrons (dotted lines), as a function of wavelength.The difference between the solid and dotted lines at lower altitudes clearly demonstrates the importance of multiple ionisations driven by the photoelectrons.bottom right panels: Contribution to the heating rate at different altitudes, H ′ [erg s −1 cm −3 ] representing the energy available before photoionisation (solid lines), the dashed line representing the heating rate due to photoionisation without the effect of photoelectrons, and the dotted lines the net heating rate taking into account secondary ionisation.Dashed and dotted lines overlap at the higher altitudes.The curves provide valuable insight into the wavelengths that contribute to heating at each altitude and how that energy is affected by photoelectrons.

Fig. 7 .
Fig. 7. Evaporating atmosphere of the case M0.05 without (red lines) and with (black lines) the effect of secondary ionisation.Left panels: from top to bottom: density, pressure, ionisation rate and temperature.Right panels: from top to bottom: Fractional ionisation, velocity, Mach number and heating rate.

Fig. 8 .
Fig. 8. Heating efficiencies for the case M0.05.Top panel: heating efficiency η λ,xe as a function of wavelength and altitude.The grey regions indicate where η λ,xe = 1.Bottom panel: heating efficiency averaged over wavelength as a function of altitude.

Fig. 9 .
Fig.9.Ionisation yields for the case M0.05.Top panel: number of total (primary plus secondary) ions created per primary photoelectron as a function of wavelength and altitude for our M0.05 planet.The grey regions indicate where Φ λ,xe = 0. Bottom panel: total (primary plus secondary) ions created per primary photoelectron averaged over wavelength and as a function of altitude.

Fig. 10 .
Fig. 10.Comparison of two different planetary masses, with a Sub-Neptune-like planet of 0.02 M J (M0.02) and a Jupiter-like planet of 0.69 M J (M0.69).Dashed lines and dotted lines represent the cases with and without secondary ionisation, respectively.From top to bottom: fractional ionisation x e as function of neutral density (a), photoionisation rate Jx Hi (b) and heating deposition rate H (c).

Fig. 11 .
Fig. 11.XUV stellar spectra.Top panel: Panchromatic SED downloaded from the MUSCLES Treasury Survey for GJ1214, GJ832 and HD97658 and scaled at a distance of 1 AU.The solar spectrum at 1 AU used in the preceding section is shown by the black line.Bottom panel: Scaled and binned spectra of our sample of stars.After scaling (see text), the integral over wavelength is the same for all four stars.

Fig. 12 .
Fig. 12. Heating efficiencies and ionisation yields for different stars and the Sun.Top panel: Heating efficiencies η λ,xe as function of wavelength and altitude for the four different spectra, GJ1214, GJ832, Sun and HD97658.Bottom panel: Ionisation yields as function of wavelength and altitude.

Fig. 13 .
Fig. 13.Heating efficiency ηr as a function of altitude averaged over wavelength for our sample of stars along with the Sun.

Table 2 .
1D spherical mass loss rates calculated for all planetary systems with and without the secondary ionisation by photoelectrons.
Table A.1.Calculated polynomial coefficient for the heating efficiency η λ,xe .