Evolution of helium triplet transits of close-in gas giants orbiting K-dwarfs

Atmospheric escape in exoplanets has traditionally been observed using hydrogen Lyman-$\alpha$ and H-$\alpha$ transmission spectroscopy, but more recent detections have utilised the metastable helium triplet at 1083$~$nm. Since this feature is accessible from the ground, it offers new possibilities for studying atmospheric escape. Our goal is to understand how the observability of escaping helium evolves during the lifetime of a highly irradiated gas giant. We extend our previous work on 1-D self-consistent hydrodynamic escape from hydrogen-only atmospheres as a function of planetary evolution to the first evolution-focused study of escaping hydrogen-helium atmospheres. Additionally, using these novel models we perform helium triplet transmission spectroscopy. We adapt our previous hydrodynamic escape model to now account for both hydrogen and helium heating and cooling processes and simultaneously solve for the population of helium in the triplet state. To account for the planetary evolution, we utilise evolving predictions of planetary radii for a close-in 0.3$~M_{\rm Jup}$ gas giant and its received stellar flux in X-ray, hard and soft EUV, and mid-UV wavelength bins assuming a K dwarf stellar host. We find that the helium triplet signature diminishes with evolution. Our models suggest that young ($\lesssim 150$~Myr), close-in gas giants ($\sim 1$ to $2~R_{\rm Jup}$) should produce helium 1083$~$nm transit absorptions of $\sim 4\%$ or $\sim 7\%$, for a slow or fast-rotating K dwarf, respectively, assuming a 2$\%$ helium abundance.


INTRODUCTION
Highly irradiated exoplanets undergo extreme atmospheric escape.This occurs mainly through the process of hydrodynamic escape, in which photoionisations of hydrogen (and, to a lesser extent, helium) atoms by X-ray and extreme-ultraviolet (hereafter XUV) photons heat the atmosphere causing its expansion and ultimately a bulk outflow of the material.This escape process can be modelled by treating the atmosphere as a collisional fluid (Yelle 2004;Murray-Clay et al. 2009).In addition to the stellar XUV heating, heating from the planet interior can also contribute to the atmospheric escape mechanism (Ginzburg et al. 2018;Kubyshkina et al. 2020;Kubyshkina & Vidotto 2021;Kubyshkina & Fossati 2022).
Intense atmospheric escape can affect a planet's overall evolution and structure.The level of atmospheric escape in highly irradiated gas giant planets declines substantially with planetary evolution (e.g.Owen 2019;Allan & Vidotto 2019;Pezzotti et al. 2021b).A reduction in the stellar XUV flux received with age and a growing planetary gravitational force are responsible for this.Failure to account for how atmospheric escape varies with evolution can lead to significant inaccuracies in the predicted atmospheric structure and escape.Strong atmospheric escape over the lifetime of a planet is one of the explanations offered for apparent trends in the detected population of exoplanets, namely the 'hot-Neptune desert' (Mazeh et al. 2016;Owen & Lai 2018) and the 'radius valley' (Fulton et al. ★ E-mail: allan@strw.leidenuniv.nl2017).Planetary evolution not only impacts the physical process of escape but additionally its corresponding observable signatures.In Allan & Vidotto (2019), we found that the hydrogen Lyman- and H- signatures of atmospheric escape weaken significantly over the evolution of highly irradiated gas giants with primordial hydrogen atmospheres.
Since our previous evolution study, there has been an explosion of detections of atmospheric escape using a new spectroscopic signature, the metastable helium triplet at 1083 nm (for theoretical studies, see Seager & Sasselov 2000;Oklopčić & Hirata 2018).This signature's ability to be observed using ground-based (in addition to space-based) facilities is greatly responsible for its popularity, and it presents a great advantage over the more traditional hydrogen Lyman- signature, which is observable only from space.Additionally, the Lyman- line suffers from strong interstellar medium absorption and contamination by geocoronal emission, rendering the line core unusable in atmospheric escape analysis.
Of the helium 1083 nm detections, K-type host stars appear to be favoured.This has been explained by their relatively low mid-UV flux, which depopulates the helium triplet state, and a high EUV flux, which populates the state through photoionisations followed by recombinations (Oklopčić 2019).Due to the favourability of K-type hosts, both models and observations of their spectra are important for interpreting escaping helium detections (France et al. 2016;Johnstone et al. 2021;Poppenhaeger 2022;Richey-Yowell et al. 2022).An important consideration is that there is no 'typical' XUV flux that fits all K dwarfs (see, e.g., Loyd et al. 2016;Youngblood et al. 2016;Sanz-Forcada 2022).In cool dwarf stars, the high-energy radiation originates in the chromosphere and corona, and they are ultimately related to the magnetic activity of the star (e.g.Vidotto et al. 2014;Toriumi & Airapetian 2022;Namekata et al. 2023).Because activity decays substantially with age, the XUV flux of K dwarfs is significantly affected with evolution (Johnstone et al. 2021;Pezzotti et al. 2021a).
In the context of predicting helium transit signatures, many of the existing models adopt a one-dimensional (1-D) Parker wind solution (Parker 1958), assuming either an isothermal atmosphere or a constant speed of sound (Oklopčić & Hirata 2018;Lampón et al. 2021;Dos Santos et al. 2022), and thus neglecting the self-consistent calculation of the energetics of the flow.However, such models often result in a degeneracy between the temperature and the mass-loss rate.Linssen et al. (2022) reduce this degeneracy, by linking isothermal Parker wind solutions with the detailed photoionisation code Cloudy (Ferland et al. 1998(Ferland et al. , 2017)), in order to rule out regions of the explored temperature -mass loss parameter space, without requiring the added complexity of self-consistently modelling the escape, by solving the coupled hydrodynamic and radiative processes consistently.The latter is instead the approach of the 1-D models of TPCI (Salz et al. 2015) and ATES (Caldiroli et al. 2021a), which is more similar to our approach.
In this paper, we compute atmospheric hydrodynamic escape solutions using a single, non-isothermal hydrodynamic model, which self-consistently incorporates energetic balance through different heating and cooling mechanisms.As will be detailed, we have modified the hydrogen-only model of Allan & Vidotto (2019) to achieve this.Our model still assumes spherical symmetry (i.e, it is 1-D).Because three-dimensional (3-D) models naturally account for asymmetries arising from a variety of processes, such as interactions with the stellar wind, magnetism and day-to-night side variations (Villarreal D'Angelo et al. 2014, 2018;McCann et al. 2019;Carolan et al. 2020Carolan et al. , 2021)), they are more realistic than 1-D models.A wealth of information on the helium signature of escape has already been gained from 3-D models (Shaikhislamov et al. 2021;Khodachenko et al. 2021a;MacLeod & Oklopčić 2021;Wang & Dai 2021b,a;Rumenskikh et al. 2022;Yan et al. 2022).However, their long computational time leaves 3-D models poorly suited for evolution studies, like ours, that require numerous models at differing ages.
Here, we model the long-term evolution of atmospheric escape of planets orbiting K dwarfs, taking into account both the evolution of the host star and the planet.The main goal is to study the evolution of the helium 1083 nm signature.To the best of our knowledge, this study is the first with this goal.Our model computes self-consistently the ionisation of hydrogen, as well as both the single and double ionisation of helium and its population in the 1 1 , 2 1  and 2 3  states.In section 2, we describe this model.Section 3 outlines our planetary evolution set-up, while section 4 showcases how the physical process of atmospheric escape varies with evolution.In section 5, we present our ray-tracing technique for simulating helium 1083 nm transmission spectroscopic transits and the resulting transit predictions with evolution.Finally, section 6 summarises the main findings of our study.

ATMOSPHERIC ESCAPE MODEL THROUGH PHOTOIONISATION
One of the main drawbacks of describing atmospheric escape using a Parker wind model is that the solution of these types of models are highly sensitive to the assumed isothermal temperature (Parker 1958).Additionally, the escape rates also depend strongly on the assumed base densities.While Parker-type models are extremely useful, allowing for quick predictions for atmospheric escape (e.g.Dos Santos et al. 2022), they omit important physics such as heating and cooling processes, which produce non-isothermal temperature structures, which can then affect the derived signatures of spectroscopic transits.
To model the evolution of atmospheric escape, we update the photoionisation-driven escape model described in Allan & Vidotto (2019, see also Murray-Clay et al. 2009) to now include helium in addition to hydrogen.Original and current model versions are both non-isothermal, 1-D, spherically symmetric and treat the escaping atmosphere as a fluid.They solve the equations of fluid dynamics in a co-rotating frame using a shooting method approach (e.g.Vidotto & Jatenco-Pereira 2006).Additionally, the new helium-incorporated version presented in this work • tracks the state of helium, accounting for transitions between the helium 1 1 , 2 1 , 2 3 , singly and doubly ionised states (displayed in Figure 1).
• Considers heating and cooling due to hydrogen and helium.
• Reads separate X-ray, mid-UV, hard and soft EUV fluxes.
• Considers photons emitted in direct recombinations and radiative decays to the helium (1 1 ) state.
Previously, we considered only photoionisations arising from a given monochromatic flux of EUV photons, implying that the input flux of EUV photons was concentrated at a representative photon energy of 20 eV.Here, we drop the monochromatic approximation and instead utilise four separate energy flux bins, each with their own representative energies (  ).Going beyond a monochromatic flux is necessary for our present work due to the differing minimum photoionisation energies of hydrogen and the considered helium states.Furthermore, the individual flux bins each have their own unique dependencies on age.Table 1 lists the hydrogen and helium photoionisations we consider and their respective cross-sections (  ) and excess kinetic energy factors (  ).The wavelength ranges covered by these chosen flux bins are X-ray (0.517-12.4 nm), hard EUV (hEUV, 10-36 nm), soft-EUV (hereafter sEUV, 36-92 nm) and mid-UV (91.2-320 nm).We set their representative energies   as 248, 40, 20 and 7 eV, respectively.Some of our bin selections are consistent with values used in other models in the literature such as Murray-Clay et al. (2009); Kubyshkina et al. (2018); Wang & Dai (2021a,b).Table 1 also lists parameters for three additional photoionising sources, photons with energies 24.6, 21.2 and 10.3 eV arising from helium state transitions in the planetary atmosphere, rather than originating directly from the assumed stellar source.Photoionisations in the model are discussed in more detail in section 2.3.

Fluid dynamic equations
We model the escaping planetary atmosphere by treating it as a nonisothermal outflow, utilising the equations of fluid dynamics.These equations enforce the conservation of momentum, energy, and mass.In steady state, the momentum equation in spherical symmetry can Table 1.Considered photoionisation processes and their respective excess kinetic energy factors   , cross-section   and their relevant reference.References correspond to: (Spitzer 1978) a , (Osterbrock & Ferland 2006) b , (Brown 1971) c , (Norcross 1971) d , (Verner et al. 1996) e , (Jacobs 1974)  be written as where  is the radial coordinate from the centre of the planet,  and  represent the atmospheric mass density and velocity, respectively, and  is the thermal pressure.The masses of star and planet are  * and  pl , and the orbital distance is .The first term on the right hand side of Equation 1 is the thermal pressure gradient while the second term represents attraction due to gravity.This equation is analogous to the momentum equation for a stellar wind (Parker 1958) with an additional term on the right side due to tidal effects.The tidal term is the sum of the centrifugal force and differential stellar gravity along the line between the planet and star (García Muñoz 2007).Conservation of energy requires that where   is the Boltzmann constant,  is the gas temperature,  = 5/3 and  is the mean particle mass.The term on the left indicates the change in the internal energy of the fluid.The first term on the right represents cooling due to gas expansion (adiabatic cooling), while the second () and third () are heating and cooling terms, which in this work accounts for a number of physical process involving hydrogen and helium species (see sections 2.3 and 2.4).This equation is, again, typically used in stellar wind models, although with different heating and cooling terms (Vidotto & Jatenco-Pereira 2006).
The conservation of mass requires that Our simulated planetary outflow originates from the sub-stellar point of the planet, which is the point closest to the star.We then apply our calculated solution over 4 steradians rendering it an upper limit to atmospheric escape (Murray-Clay et al. 2009;Johnstone et al. 2015b).Therefore, from Equation 3, we have that the mass-loss rate of the escaping atmosphere is  = 4 2 .
We also solve the equation of ionisation balance of hydrogen This equation states that the advection of hydrogen ions (first term) is balanced by the combined rate of photo-ionisation (Φ H 0 , see Equation 10and Table 1) and collisional-ionisation (Ψ H 0 , Table 2) of neutral hydrogen with the radiative recombinations (  H 0 , also see Table 2).Altogether, the fluid dynamic equations (Equations 1-3), the hydrogen ionisation balance equation (Equation 4) and the four additional helium population equations (see next section, Equations 5 -8) form a system of coupled differential equations.There is only one physical solution for this system: a transonic wind beginning at the base of the atmosphere with a subsonic velocity, passing through the sonic point and thereby reaching supersonic velocities (Parker 1958).Contrary to isothermal wind models, we do not know a priori the location of the sonic point nor the initial wind velocity, therefore we follow a numerical approach in solving the coupled differential equations utilising a 'shooting method' combined with a fourth order Runge-Kutta solver and an iteration scheme.Our initial solution assumes a hydrostatic atmosphere for the density, and the shooting method finds the solution that passes through the sonic point.This new solution provides a different density profile than initially assumed, so using this newly found density structure, we again search for the solution that passes through the sonic point through the shooting method.In this way, the density, temperature, velocity and optical depth profiles are updated with each set of iterations of the shooting method, iteratively approaching to a converged solution.We consider this solution to have been reached once the predicted mass-loss rate and terminal velocity of two subsequent model results agree to within 1%.All atmospheric profiles are calculated with a resolution of 10 −6  pl below the sonic point, and 2 × 10 −5  pl above.
The necessary free parameters of our model include the temperature  0 = 1 000 K and density  0 = 4 × 10 −14 g cm −3 at the base of the atmosphere, chosen to be at 1  pl .Allan & Vidotto (2019) and Murray-Clay et al. (2009) found that large variations in these values had a negligible effect on the resulting simulated escape.This is similar in our models as will be further discussed in Appendix A. Our model must also be given initial guesses for the abundances by number at the base of the atmosphere for ionised hydrogen, singlyand doubly-ionised helium as well as its 1 1 , 2 1  and 2 3  states.These initial guesses were 10 −5 , 10 −5 , 10 −7 , ∼ 1, ∼ 0 and ∼ 0, meaning that initially both hydrogen and helium are predominantly neutral, with helium being nearly entirely in its singlet state.To ensure the validity of the hydrodynamic escape we model, once the solution is found, we confirm that the atmosphere remains within the collisional regime with a Knudsen number (ratio of mean free path to atmospheric scale height) < 1.In the non-collisional regime, particle models such as those used in Bourrier & Lecavelier des Etangs (2013); Bourrier et al. (2016); Allart et al. (2018); Spake et al. (2018); Allart et al. (2019) are better suited to model atmospheric escape.Our models are computed all the way to 10  pl , but we note that in reality, the extent of the atmosphere could be larger or smaller, as its boundary should be set by the interactions with the wind of the host star (e.g.Vidotto & Cleary 2020), neglected in our study.

Calculating helium state populations
In addition to Equations 1 to 4, we also solve four continuity equations for the populations of helium, in the helium singlet ground state, He(11 ) the helium triplet state, He(2 3 ) the helium 2 1  state, He(2 1 ) and the singly ionised helium state, He These state that the fractions of helium in its 1 1 , 2 3 , 2 1 , and singly ionised states  1 1  ,  2 3  ,  2 1  and  He + respectively, are determined by the balance of the processes directly populating and depopulating these states.These processes (de-)populating the helium states are defined as follows • Φ -photoionisations (out of the state in parentheses) by capable stellar XUV and mid-UV photons as well as by photons released in the planetary atmosphere through helium transitions, •  -collisional ionisations (out of the state in parentheses), •  -recombinations (into the state in parentheses), •  -collisional excitation with a free electron, •  -collisional excitation with H 0 , •  -radiative decays, •  -charge exchange.
Table 2 lists all of their corresponding rate equations for these helium states in addition to that for hydrogen.Figure 1 visualises the considered paths for the helium states.The helium number densities in each state are obtained by multiplying their respective fraction by the total number density of helium  He .Our model considers free electrons produced from both hydrogen and helium ionisations while maintaining a neutral net charge for the global modleled atmosphere.We also calculate the fraction of doubly-ionised helium, assuming its number density to be  He ++ =  He −  1 1  −  2 3  −  2 1  −  He + .However, including He ++ is found to have a negligible effect on the modelled hydrodynamics and observability of atmospheric escape.
Being a non-metastable and hence short-lived state (Wiese & Fuhr 2009), our model does not solve for the fraction of helium in the 2 1  state.Rather, we follow Oklopčić & Hirata (2018) in allowing electron collisional excitation from 2 3  to 2 1 ,  He(2 3  → 2 1 ) to populate the 1 1  state on account of the rapid decay,  He(2 1  → 1 1 ) .Lampón et al. (2020) and Dos Santos et al. ( 2022) also follow this assumption, which most importantly accounts for the depopulation of the 2 3  state.Additionally, our model considers recombinations into 2 1 ,  − He(2 1 ) , which is also set to populate 1 1  in our model's helium equations.It should be noted that the rapid radiative decay from the 2 1  to 1 1  state  He(2 1  → 1 1 ) releases a 21.2 eV photon (Eikema et al. 1996;Sun & Hu 2020), capable of photoionising hydrogen, helium 2 3  and 2 1  in our model 1 .As will be discussed in section 2.3, our model considers such photoionisations.
We use Case-B coefficients in our calculations of hydrogen and the helium 2 3 , 2 1 , 2 1  states.Direct recombinations (i.e.Case-A − Case-B) are not considered for these mentioned states, given their low energies and the high probability for the photon released in the process to photoionise a species identical to that formed in the recombination, resulting in no net variation in the ionisation fraction.Direct recombinations to He(1 1 ) state however are different in that they release a higher energy 24.6 eV photon (Benjamin et al. 1999;Osterbrock & Ferland 2006), capable of photoionising both the more abundant hydrogen in addition to the various helium states.Accordingly, we calculate both direct and non-direct recombinations for He(1 1 ), in order to obtain the rate of the 24.6 eV photoionisations, as will be discussed in section 2.3.
In short, our model considers hydrogen and helium photoionisation, collisional ionisation, recombination, collisional (de-) excitation, charge exchange and radiative decays.Our model solves the four helium populations equations (Equations 5 to 8) simultaneously with Equations 1 to 4.
To model the escape of a primordial hydrogen and helium atmosphere, the helium to hydrogen number abundance must be known.Models have typically assumed this to be constant throughout their atmosphere, which we also assume.Some typically derived values are: He/H ∼ 0.02/0.98 for HD209458b (Lampón et al. 2020;Khodachenko et al. 2021b); He/H ∼ 0.016/0.984(Shaikhislamov et al. 2021) and He/H ∼ 0.015/0.985(Lampón et al. 2021) for GJ3470b; He/H∼ 0.008/0.992(Lampón et al. 2021) and He/H∼ 0.005/0.995(Rumenskikh et al. 2022) for HD189733b.All of the mentioned he-  2 for the rate equations and their relevant references.lium abundances fall below the solar value of approximately 0.1.Models of WASP-107b were the first to obtain a He/H number abundance close to the solar value, with a predicted helium abundance of 0.075-0.15(Khodachenko et al. 2021a).In our models, we use helium number abundances of 0.02 and 0.1 to investigate how this affects the predicted transits.

Photoionisation and heating processes with hydrogen and helium
Photoionisations are not only important for determining the population of atmospheric hydrogen and helium states.If a given photon has an energy in excess of the photoionisation energy, this excess is transferred to an electron in the form of kinetic energy.This collisional process acts to heat up the atmosphere.The amount of photoionisations and heating that is deposited in the atmosphere depends on the individual optical depths   for each of the considered photoionisations where  is the number density of the absorber (neutral hydrogen, helium singlet or triplet) and   is the absorber-and wavelengthdependent photoionisation cross section (Table 1).The integration is performed for distances  starting at the top of the atmosphere (10 pl , in our models) until the bottom of the atmosphere at 1 pl .
The rate of one of the considered photoionisations is then where   is the energy flux at   and  sp, is a weighting factor we introduce to prevent the same photon from being able to photoionise both hydrogen and helium, where 'sp' refers to the specific species being photoionised.The fraction of photons which photoionise hydrogen rather than helium from a given flux bin with   is (see Equation 2.21 of Osterbrock & Ferland 2006) We use similar weighting factor expressions for each considered photoionisation process in our model.For example, Table 1 shows that the X-ray flux bin can photoionise H 0 , He(1 1 ) and He + .Hence, the weighting factor for an X-ray photon to photoionise He(1 1 ),  He(1 1 S),X−ray in our model is given by This important consideration to prevent from double-counting photons is often not included in hydrogen-helium atmospheric escape models.The rate of the heating associated with a photoionisation is calculated as Combined, the terms    sp,  −   account for the remaining supply of photoionising-capable photons at a given distance into the atmosphere.  = 1 − ( ion /  ) accounts for how much energy the photon with   has in excess of the photoionisation energy cost ( ion =13.6, 24.6, 4.0 4.8, 54.4 eV for H 0 , helium 1 1 , 2 1 , 2 3  and He + , respectively).Values of   are listed in Table 1.Finally, the combined    terms account for the availability of potential absorbers.The total photoionisation heating , included in Equation 2, is then the sum of the heating due to each of the considered species.
In addition to the considered stellar flux, our model also accounts for 24.6 eV and 21.2 eV photons released in recombinations and radiative decay from He(2 1 ) to He(1 1 ).It also accounts for the release of two 10.3 eV (Bergeson et al. 1998) photons released in Table 2.The populating and depopulating processes for the hydrogen and helium states considered in our model.Densities are given in cm −3 and temperature in K. He(2 1 ) * is marked with an asterisk as a reminder that we do not solve for the fraction of helium in the 2 1  state, as it is short-lived (see Section 2.2).   * He(1 1 ) is marked with an asterisk to draw attention to the removal of helium 2 1  and 2 1  recombinations from this term.Letters correspond to the following references: (Brown 1971) a , (Norcross 1971) b , (Verner et al. 1996) c , (Spitzer 1978) d , (Osterbrock & Ferland 2006) e , (Benjamin et al. 1999) f , (Wiese & Fuhr 2009) g , (Eikema et al. 1996) h , (Bergeson et al. 1998) i , (Hui & Gnedin 1997) assumed here the radiative decay from helium's 2 1  to 1 1  state.We assume the photoionisation rates by these photons to be given by the rate of the process from which they are released (× 2 for the two-photon process), as shown in Table 2.In doing so, we assume that all of these photons are re-absorbed by the planetary atmosphere and do not escape.Given that these photons are predominantly released in the optically thick inner-most region of the atmosphere, we consider this assumption to be reasonable.To determine the fractions of which of the viable atmospheric species are absorbed by these photons, we again utilise the weighting factor given by Equation 11.For example, the modelled rate of the released 24.6 eV photons which photoionise hydrogen is  H 0 ,24.6 eV  − He(1 1 ) .As will be shown, this contributes non-negligibly to the overall photoionisation of hydrogen and heating in the inner-most region of the modelled planetary atmosphere.
Cecchi-Pestellini et al. (2006) showed that X-rays can play a large role in heating the upper layers of planetary atmospheres, particularly those which orbit closely to young (high X-ray flux) stars.Owing to their large photon energies, primary X-ray photoionisations are followed by subsequent photoionisations caused by released energetic secondary photo-and Auger-electrons (Güdel 2015).Heating by Xrays is dominated by such secondary photo-electron generation from the K-shells of metals, with oxygen and carbon being the most important (Owen & Jackson 2012).Few works (Cecchi-Pestellini et al. 2006;Shematovich et al. 2014;Locci et al. 2022) have modelled the secondary electron photoionisation cascade.As our model considers only hydrogen and helium and not heavier elements important for such processes, secondary photoionisations are omitted from our model (as is the case for the models of Oklopčić & Hirata 2018;Lampón et al. 2020;Dos Santos et al. 2022).Hence, it should be noted that our models may underestimate the number of photo-inisations by X-rays, particularly when the host star is a young, initially fast rotator.A recent work by Gillet et al. (2023) studies the effects of secondary ionisation by photoelectrons and stresses that such ionisations should be taken into account in 1-D hydrodynamic modelling, reporting a reduction in the predicted mass-loss rate of 43% upon its inclusion for a 0.69 Jup planet orbiting the K1V host HD97658 at 0.074 au.

Cooling processes with hydrogen and helium
In addition to heating, our energy equation (Equation 2) also includes various hydrogen and helium cooling contributions.These are collisional ionisation, collisional excitation, recombination and Bremsstrahlung.The volumetric cooling rates were obtained from Black (1981) and Cen (1992) and are listed in Table 3.To obtain the stated cooling contribution processes (3.4) and (3.10), we rearranged the respective rate equations given in Table 3 of Black (1981) so as to use our own helium triplet number density rather than an approximation given by Equation 11of Black (1981).Rearranging and substituting this equation into their rate equations gives the cooling rate equations listed in Table 3.The sum of all the volumetric cooling rates listed in Table 3 results in , which is included in Equation 2.

Test case: solving helium populations in post-processing
For comparative purposes, we also present a representative model in which the helium equations are solved in a post-processing step (as is the case for the models of Oklopčić & Hirata 2018;Lampón et al. 2020;Dos Santos et al. 2022).If solved in a post-processing step, helium heating and cooling contributions can not be included in the energy equation (Equation 2) as to calculate these contributions requires knowledge of the helium populations.Hence, whether the helium population equations are solved simultaneously with the fluid dynamic equations or in a post-processing step can impact the predicted atmospheric escape, the results of which will be discussed in sections 4.2 and 5.4.
Our post-processing test model • assumes an initially fast rotating stellar host.
• Assumes a constant helium number abundance of 10%.
• Allows heating by photoionisation of hydrogen only.
• Allows photoionisation by stellar XUV and mid-UV photons only, ignoring photons released in the planetary atmosphere as a result of helium transitions.
• Neglects photoionisation weighting factors (see Equation 11), meaning that a single photon could possibly photoionise hydrogen, as well as any of the viable considered helium states in the model.
• Includes only three cooling processes: recombination, collisional ionisation and excitation of hydrogen.
• Considers only electrons arising from hydrogen photoionisations.
While these terms are chosen to be similar to the post-processing 1-D models of Oklopčić & Hirata (2018); Lampón et al. (2020); Dos Santos et al. (2022), the specific modelling differs from each of these works in numerous ways.For example, they assume a Parker wind in their modelling of the atmospheric escape and hence do not calculate the heating due to photoionisations.However, p-winds (Dos Santos et al. 2022) allows the user to input a more complex and selfconsistent atmospheric structure, rather than assume a Parker wind, as done for HD 189733 b in Dos Santos et al. (2023).Additionally, there are some minor variations in a number of the assumed rates between states (see Table 2).

MODEL INPUTS: EVOLUTION OF HIGH-ENERGY STELLAR FLUX AND PLANETARY RADII
Allan & Vidotto (2019) previously showed that the evolution of atmospheric escape for a close-in planet depends on two important factors: (i) as the host star evolves, its activity declines due to spin down (Skumanich 1972;Kawaler 1988;Vidotto et al. 2014), resulting in declining fluxes in the XUV (Jackson et al. 2012;Johnstone et al. 2015a;Johnstone et al. 2021) and (ii) as the planet evolves, cooling causes it to contract with time (Fortney & Nettelmann 2010).
The level of atmospheric escape and consequently the observational signatures of escaping hydrogen in the Lyman- and H- lines were found to vary strongly with the evolution of the modelled planets, with younger planets exhibiting greater escape and deeper absorptions in both lines.This is the result of a favourable combination of higher irradiation fluxes and weaker gravities at young ages.In a continuation of this work, we now study how the helium 1083 nm signature evolves over the lifetime of a planet.
In order to incorporate the decline in XUV flux in this current work, we look to the physical rotational evolution model of Johnstone et al. (2021) that is constrained by observed rotation distributions in young stellar clusters.Their study offers publicly available evolutionary tracks of stellar flux for a wide variety of stellar mass, age and initial rotation.We use their tracks of X-ray, hEUV, sEUV flux and stellar radius evolution for a 0.7- ⊙ star, chosen to be representative of a K-dwarf star as detections of 1083 nm seem to favour such a host.We obtain tracks for both their 'slow' and 'fast' initial stellar rotation definitions corresponding to the 5th and 95th percentiles of their observed 150 Myr rotation distribution.In Allan & Vidotto (2019) and where necessary (such as in the mostly inaccessible EUV wavelength range) stellar spectral models.While other K dwarf spectra were obtained by MUSCLES, we choose HD85512 to perform this normalisation based on its relatively late age.By this age, the slow and fast rotators have converged to the same rotation rate, meaning that normalising the flux tracks at this age retains the original track ratios between slow and fast rotator.We assume an age of ∼ 5600 Myr  [Lower-panel] Photo-ionisation cross sections for hydrogen and considered helium states.Note that while He(2 1 ) photoionisations are not considered in the model, its cross-section profile is shown here for comparison to that of the 2 1  and 2 3  helium states.The four shaded regions indicate the wavelength channels of X-ray, hEUV, sEUV and mid-UV going from left to right.The representative energies of each of these flux bins is marked by the vertical solid white lines in the lower-panel.
for HD 85512 (Pepe et al. 2011) and a distance of 11.28 pc (Gaia Collaboration et al. 2016Collaboration et al. , 2022) ) to calculate its surface flux in each of the considered bins (see star symbols in Figure 3).their near-and far-UV fluxes into a single 'mid-UV flux'.We include the mid-UV flux due to the ability of these photons to photoionise helium out of the triplet state.
The central-panel of Figure 3 shows the evolution of the stellar radius, also obtained from the model of Johnstone et al. (2021).Our evolving planetary radii inputs are shown in the lower-panel of Figure 3.These are the same radii used in Allan & Vidotto (2019), corresponding to a 0.3  jup gas giant planet orbiting a solar-like star at 0.045 au (Fortney & Nettelmann 2010).Note that these models neglect planetary-radius inflation.The close proximity to its host and the planetary mass puts this planet at the upper edge of the hot-Neptune desert (Mazeh et al. 2016).It should be noted that there is an inconsistency between the assumed stellar parameters in our flux (K dwarf) and radius evolution input (which assumes a G dwarf).
In the future, it would be interesting to couple our hydrodynamic escape model to a planetary thermal evolution model as is done in Kubyshkina et al. (2020).This would allow us to also consider the decrease in planetary mass, which could be significant in case of strong atmospheric escape over time.This mass reduction would act to reduce the planetary gravity with evolution.As discussed in Allan & Vidotto (2019), omitting this shrinking mass with evolution in our modelling leads to our models slightly underestimating the true atmospheric escape progressively with planetary evolution, as the planetary gravity is progressively overestimated.This mass loss underestimation is minor due to the stronger dependence of the gravitational force with planetary radius compared to planetary mass, and planetary radius varying more than the mass with evolution.

THE EVOLUTION OF ATMOSPHERIC ESCAPE
Table 4 describes the four sets of evolution-sampled models we consider.This shows the initial stellar rotation, the He/H number abundance and whether or not helium heating and cooling processes are included.Hydrogen heating and cooling processes are included in all models.Helium heating and cooling, as well as the helium populations are self-consistently included in the top 3 models, while the model F10%PostProc does not include any helium heating or cooling terms (in the hydrodynamics model, the presence of the helium particles only affects the mean molecular weight of the gas) -in this case, the helium population is computed in a post-processing step.All models share the following common parameters: stellar mass  * = 0.7 ⊙ , orbital distance  = 0.045 au and planetary mass  pl = 0.3  Jup , while the received stellar flux and the planetary and stellar radii evolve as given by Figure 3.
Figure 4 displays how the predicted mass-loss rate (upper-panel), cumulative mass lost (central-panel) and terminal velocity of the escaping atmosphere (lower-panel) vary as a function of age for each of our model sets.We calculate the cumulative mass lost by integrating the mass-loss rate profile with respect to planetary age.As found in Allan & Vidotto (2019), the diminishing XUV flux required to heat the planetary atmosphere combined with the growing gravitational force due to the shrinking planetary radii lead to the weakening of atmospheric escape as the planet evolves across all models.Variations with evolution are more extreme for the massloss rate compared to the terminal velocity.As discussed below, this is the result of the majority of heating occurring in the sub-sonic region of the atmosphere.
To investigate the evolution of the energetics of our models, the upper-panels of Figure 5 show how the considered photoionisation heating processes differ at young (left) compared to old (right) planetary ages as a function of distance.The lower-panels show the domi-   4. nant cooling processes.The displayed individual heating and cooling contributions correspond to our F2% model (see Table 4). and  correspond to the total heating source and cooling sink contributions (see Equation 2) of this model.In each panel, the outflow velocity of the F2% model remains sub-sonic within the shaded regions, reaching supersonic velocities beyond.For comparison, we show also the total heating and cooling contributions of models F10% and F10%PostProc.It should be noted that the F10%PostProc heating profile could not simply be reproduced by summing the displayed hydrogen-only heating contributions of the F10% model.Firstly, this is because we account for the limited availability of photons as explained previously with Equation 11.Secondly, the heating affects the density distribution of neutrals which in turn affects the heating contributions.
We see in Figure 5 that the same heating process remains dominant above ∼ 1.1 pl throughout the full planetary evolution: heating due to the photoionisation of neutral hydrogen by sEUV photons (ID: 1.1 of Table 1).Below ∼ 1.1 pl , heating from hydrogen photoionisations by hEUV photons (ID: 1.2) dominate over sEUV photons for the young planet, whereas at the oldest age, hydrogen photoionisation by the 24.6 eV photon (ID: 1.4), itself ejected in a direct recombination Table 4. Description of the evolution model sets presented in this paper.Rotation refers to the initial stellar rotation (fast rotators have higher high-energy fluxes during their youth).The third column indicates the helium number abundance, with the remainder of the atmosphere being hydrogen.The listed heating and cooling processes correspond to the row IDs of Tables 1 and 3. Model F10%PostProc does not incorporate helium energetics in the hydrodynamic escape (although the mass of the helium particles affect the mean molecular weight) and instead computes the helium state fractions post-processingly, as described throughout section 2. All models share the following common parameters;  * = 0.7 ⊙ ,  = 0.045 au,  pl = 0.3  Jup .The received stellar flux and radius and the planetary radius vary with evolution following Figure 3 to He(1 1 ), becomes the dominant heater below ∼ 1.1 pl .We found that the contribution of photoionisation of singlet state helium to the heating process, although not dominant, is not negligible for abundances of 2%, and is more important for higher abundances.Sections 4.1 and 4.2 will discuss the resulting effects of altering the assumed He/H number abundance and omitting helium energetics.
Heating arising from photoionisations out of the helium triplet state is negligible.
Clearly, the majority of heating affects the sub-sonic region of the atmosphere.The heating contribution peaks within the sub-sonic region for the displayed young and old models.Accordingly, the temperature profile follows a similar behaviour, reaching peaks of approximately 11 250 K and 7 750 K respectively as will later be shown in Figure 7. Due to the heating being mostly deposited in the sub-sonic region, heating variations with planetary evolution have a greater effect on the mass-loss rate, while having a lesser effect on the velocity of the outflowing atmosphere.This is indeed verified by Figure 4, with a decrease of over one order of magnitude in massloss rate while the terminal velocity decreases by 35% over the full evolution of our F2% model.
The lower-panels of Figure 5 show the dominant cooling processes.Similar to the heating processes, the net cooling peaks in the sub-sonic region of the atmosphere.Adiabatic cooling due to the expansion of the escaping atmosphere is more important than the cooling sink terms (composing  in Equation 2) for the old 5 000-Myr model.Although not shown, this is also true for models F10% and F10%PostProc at the same age.At the young age of 16 Myr however, we see that collisional excitation of neutral hydrogen (dash-dotted green, ID: 3.9 in Table 3) more commonly referred to as Lyman- cooling dominates at distances between 1.1 to 1.3  pl .Outside of this distance cooling due to adiabatic expansion again dominates.In their modelling with ATES, Caldiroli et al. (2021b) show that adiabatic expansion dominates the atmospheric cooling in the case of their low-irradiation model of HD 97658 b, similar to our 5 000 Myr model, while for their high-irradiation planet WASP-43 b, radiative cooling, predominantly Lyman- cooling, dominates below 1.5  pl , similar to our 16 Myr planet.Our model's predicted heating and cooling behaviour also agrees reasonably well with that of Zhang et al. (2022b) obtained with TPCI (Salz et al. 2015) for the young mini-Neptune TOI 560.01.Assuming solar metalicity, they find adiabatic cooling to be more important than radiative cooling at nearly all radii, with radiative cooling overtaking in the ∼ 2−3  pl region of the atmosphere.Their predicted heating is also dominated by neutral hydrogen photoionisations with smaller but non-negligible contributions from helium in agreement with our findings.Below 1.3  pl , line heating and photoionsation of metals dominate their modelled heating.
Figure 6 shows the rates of the processes that directly populate (solid lines) or depopulate (dashed lines) the helium triplet state.Left and right panels correspond to ages 16 and 5 000 Myr, respectively.Again, only the dominant processes are displayed in the figure, the exhaustive list of considered processes is reserved for Table 2.While the relative rates between processes vary slightly with age, it is clear that the same triplet (de-)populating processes dominate throughout the planetary evolution.This also holds true at the intermediate ages not shown.At all distances, the dominant triplet-populating process is the recombination of ionised helium into the triplet state,   [He(2 3 )].
Mostly balancing this process, is the de-excitation of triplet-state helium into the 2 1  state by means of a collision with a free electron ( He(2 3  → 2 1 ) ).The general behaviour of the mentioned processes agree with the 1-D model of Oklopčić & Hirata (2018) and the 3-D models of Khodachenko et al. (2021a) and Rumenskikh et al. (2022) for WASP-107b and HD189733b, respectively.The depopulation of the triplet state due to both the photoionisation from mid-UV and sEUV photons is not a dominant process except in the very outer, more tenuous atmosphere, consistent with the findings of Oklopčić (2019) who shows that for stars cooler than spectral type G, collisional de-excitation from triplet to singlet state dominate.This is due to cooler K-dwarf stars having a smaller ratio of mid-UV : EUV flux, meaning a smaller ratio of triplet depopulating : escape driving flux, compared to warmer G type stars.
The triplet number density profile resulting from the balance of the mentioned populating and depopulating processes is shown in Figure 7, at young (upper-panel) and old (lower-panel) ages.The number densities of the other modelled helium states as well as neutral and ionised hydrogen (dotted profiles) are also included.The atmospheric temperature is shown by the solid profile.The temperature reached by the younger planet's atmosphere is larger due to the higher level of XUV flux received: the atmosphere of the 16-Myr old planet reaches a maximum temperature of 11 250 K, while the atmosphere of the 5-Gyr-old planet reaches a maximum temperature of 7 750 K.The level of hydrogen and helium ionisation is also larger at younger ages on account of the relatively larger incident flux.The number density of helium in the triplet state is greater for the younger planet.This is due to a higher rate of helium recombinations due to a larger availability of ionised helium at younger ages (Figure 6).

Effect of helium abundance on hydrodynamics
To investigate how the abundance of helium affects the hydrodynamic escape models, we run two set of models where we adopt helium number abundances of 0.02 (F2%) and 0.1 (F10%).
Volumetric heating rates (upper-panel) and volumetric cooling rate (lower-panel) profiles at young (left) and old (right) planetary ages.Shaded regions on the left of each panel highlight the sub-sonic region of the atmosphere.In both legend sets, the numbers in parentheses relate each process to their corresponding process number in Table 1 for heating and Table 3 for cooling.In all displayed panels there are processes that fall below the chosen -axis range and hence are not included in the legend.The tables on the other hand list all of the included processes.The displayed individual heating and cooling contributions correspond to our F2% model (see Table 4). and  correspond to this model's net heating source and cooling sink contributions (see Equation 2).For comparison, we show also the  and  of models F10% and F10%PostProc.
In Figure 5, individual heating contributions were shown only for the F2% model, with only the total  profile of F10% on display.A deeper analysis into the heating for the F10% model reveals however that heating due to helium unsurprisingly plays a more significant role for a greater helium abundance (see Appendix B and Figure B1).
The nature of the atmospheric escape with evolution for models F2% and F10% is shown by their mass-loss rate and terminal velocities given in Figure 4 by the solid (F2%) and dash-dotted (F10%) tracks.Clearly, they both lead to very similar mass-loss rates, with both models resulting in similar cumulative loss of 15% (F10%) and 16% (F2%) of the planet's initial mass over its entire evolution.The larger helium abundance produces also a slight decrease of 5-9% in terminal velocity maintained throughout the evolution as seen in the bottom panel of Figure 4. Slower terminal velocities are reached due to the heavier mean molecular weight combined with similar total heating and cooling profiles (see Figure 5).In conclusion, we find   4).The upper-panel corresponds to an age of 16 Myr while the lower-panel shows 5 000 Myr.The planetary parameters are that of our F2% model (see Table 4).Similar figures for our F10% and F10%PostProc are given in the Appendix.
that increasing the assumed helium abundance leads to slightly slower outflows, which undergo slightly less atmospheric mass losses.While the He/H range we consider here appears only to affect the hydrodynamics of escape trivially, there are two important effects resulting from the choice of the abundance.Firstly, the observability of helium transit is strongly dependent on the choice of helium abundance (Section 5.3).Secondly, if we omit the calculation of helium energetics, the choice of helium abundance can significantly affect the atmospheric escape.This will now be discussed.

Effect of helium energetics on hydrodynamics
The inclusion of helium affects the atmospheric dynamics in two ways.Firstly, it raises the mean molecular weight above that of a hydrogen-pure atmosphere.Intuitively, this effect acts to weaken atmospheric escape by raising the gravitational force experienced by the atmospheric material.This escape-reducing effect is straightforward to include in the modelling of atmospheric escape.Secondly, and more complex to model, is the effect of additional heating processes due to the inclusion of helium.In our modelling, we selfconsistently incorporate the helium energetics, i.e. heating due to photoionisation from the 1 1 , 2 1 , 2 3  and singly ionised states, as well as cooling due to collisional excitation, and ionisation, recombination, and Bremsstrahlung in the fluid dynamic equations (Equations 1-3).To investigate the effects that self-consistently including helium energetics has, we compute the model set F10%PostProc which omits such effects for comparison, as described previously in section 2.5.
The top panels of Figure 5 show the total  heating profiles of the model sets F10% (dash-dotted fuchsia) and F10%PostProc (dashdotted purple).Comparing these, we see that close to the planet, the heating rates are comparable, but further away, the heating profile of the F10%PostProc model is considerably less.This leads to similar peak atmospheric temperatures, with models F10% and F10%PostProc reaching 11.2(8.8)kK and 11.3(8.3)kK at ages 16 (5 000) Myr, respectively, but a faster temperature decay with distance when omitting helium energetics (see Appendix C, Figure B2).In order to better understand how the contribution of heating is affected by the specific model set-up, we calculate the net heating contributions for each model by integrating their volumetric heating rates over volume.Doing so reveals a ∼21% heating reduction for the F10%PostProc model relative to the F10% model at the youngest age, growing to 45% for the oldest model age of 5 Gyr.
As reduced atmospheric heating drives less escape, our F10%PostProc model (purple dotted tracks in Figure 4) undergoes less mass loss compared to the F10% model (fuchsia dash-dotted), with respective losses of 13 and 15% of their initial mass over their evolution.The terminal velocity is also reduced by the omission of helium energetics, with reductions between 2.5 to 5.5% throughout the planetary evolution.The main barrier to the inclusion of helium energetics is that the helium populations must be solved simultaneously with the fluid dynamics equations rather than in a post-processing step.This is because knowledge of the helium populations are required for calculating the helium heating and cooling terms, which ought to be included in the energy conservation equation.It is important to note however that even our model without helium heating is more sophisticated than Parker-type wind models, which assume ad-hoc constant temperatures throughout the atmosphere.

The ray-tracing model for transmission spectroscopy
Here, we model transmission spectroscopy of the helium triplet at 1083 nm to investigate how the observability of this signature varies with long-term evolution.Our ray-tracing technique for modelling transmission spectroscopy is a helium-adapted version of a previous model described in Allan & Vidotto (2019, Lyman-𝛼 and H-𝛼) and Vidotto et al. (2018, O I 130.22 nm).Hence, this section gives only a brief summary aimed at documenting the necessary modifications.For a more thorough description of the general model we refer the reader to Allan & Vidotto (2019).
As our ray tracing model is 3-D, we symmetrically fill a 3-D grid centred on the planet with our 1-D hydrodynamic calculations of atmospheric temperature, velocity and density.One of the grid axes is aligned along the observer-star line, so that the grid seen in the plane of the sky is a square of 201 × 201 cells.We also account for the line-of-sight velocity related to the orbital motion of the planet (assuming circular orbit) -this correction is more important for phases farther from mid-transit, as in mid-transit, the line-of-sight velocity of the planet's orbital motion goes to zero.
The helium triplet is comprised of three individual lines, with line-centre wavelengths in air of 1082.909,1083.025,1083.034nm.Hence, we model three individual wavelength-dependent transits for each of the triplet lines following the process described in detail in Allan & Vidotto (2019) for Lyman- and H- lines.We use the NIST database 2 (Kramida et al. 2022) to obtain the following line properties: the Einstein coefficient is   = 1.0216×10 7 s −1 , the mentioned line centre wavelengths and their corresponding oscillator strengths   , 0.059902, 0.17974, 0.29958.With these atomic parameters, we calculate a wavelength-dependent optical depth along a single ray in the direction connecting the observer to the star-planet system.In doing so, we use a Voigt line profile, a convolution of a Gaussian and Lorentzian line profile accounting for both Doppler and natural broadening where  0 is the wavelength at line centre for the specific helium 2 https://physics.nist.gov/asdtriplet line,  =    0 /(4 th ) is the damping parameter,  th = (2 B / He ) 1/2 is the thermal velocity with  He being the mass of atomic helium.The velocity offset from the line centre is Δ = channel −  LOS , where  LOS is the line of sight flow velocity of the escaping wind and channel represents the velocity 'channel' (related to wavelength via Doppler shifts) of the measurement.We slice our velocity calculations in 71 channels from -100 km/s to +100 km/s.In calculating this line profile, we make use of IDL's inbuilt voigt function.From the wavelength-dependent optical depth, we calculate the wavelength-dependent absorption.By integrating this over all considered rays within the stellar disk, we find the total absorption.In this way, we calculate the transit depth contributions due to each of the three lines that make up the helium triplet.Finally, we sum these contributions to obtain the total helium triplet transit depth.
We assume a transit along the centre of the stellar disk (i.e., no impact parameter).We also neglect centre-to-limb variations in the stellar disc.With  * = 0.7 ⊙ ,  = 0.045 au and the age-specific  * shown in Figure 3, we obtain a transit duration between 2 to 2.4 hours for our models (see Equation 17of Allan & Vidotto 2019).

The evolution of helium 1083 nm transmission spectroscopy
Figure 8 shows the plane-of-the-sky extinction computed for each of the individual lines of the helium triplet (shown here at midtransit only).In each panel the stellar disk is outlined by the orange dashed circle.By integrating the extinction, we compute the excess absorption shown in Figure 9, where colour corresponds to various transit phases.At younger ages (upper-panel), the excess transit is larger and the transit depth is greatest at mid-transit as this is when the most dense atmospheric material occults the stellar disk.There is an asymmetry between pre-and post-mid-transit spectra as a result of the Doppler velocities of the atmospheric material due to the orbital motion.This asymmetry is greatest at times furthest from mid-transit as the majority of obscuring material is red-shifted (transit phase=+0.5,fourth contact 'T4') or blue-shifted (transit phase=-0.5,first contact 'T1').The degree of asymmetry would be even greater had we considered the interaction with a stellar wind (Carolan et al. 2021), as the outflowing atmosphere would be asymmetric, likely featuring a comet tail-like structure.
In practice, transit observations are integrated over a certain amount of time to obtain a sufficiently high signal-to-noise ratio.As the helium 1083 nm signature is then averaged out over time, the absorption from such an observation is weaker than the theoretical instantaneous 1083 nm signature at mid-transit.We account for this time-sampling effect by averaging our absorption profiles between T1 and T4, following Dos Santos et al. (2022).The black dashed spectra in Figure 9 show these calculated mean averages.
The upper-panel of Figure 10 shows the (T1-T4 phase-averaged) helium 1083 nm transmission spectra along planetary evolution.The lower-panel offers a closer look at the oldest age spectrum and its constituent contributions from each of the individual helium triplet lines.Figure 11 displays the evolution of the T1-T4 phase-averaged equivalent width (upper-panel) and peak excess absorption (lower-panel) of these helium 1083 nm signatures for all four of our model sets.In our models, thermal broadening is negligible and the broad profiles and large equivalent widths, especially more noticeable at younger ages, are due to Doppler broadening.This is because there is more absorbing material at larger distances (and thus with larger velocities) in the 16 Myr old model than in older systems.Additionally, we see that the helium 1083 nm signature weakens with planetary evolution.The main cause of the larger 1083 nm absorption at younger ages is  4).
the larger amount of material in the triplet state obscuring the stellar disk during planetary transit (see, e.g., the density profiles in Figure 7 and the extinction in Figure 8).Interestingly, we see in Figure 10 very little variation in the 1083 nm absorption profile between ages 43 and 150 Myr, despite the reduction in atmospheric escape.This is because the stellar radius also drops from 0.74 to 0.64  ⊙ , hence the ratio of obscuring He(2 3 ) to the stellar disk area remains similar at both ages.

Effect of helium abundance on 1083 nm observability
When comparing models with different abundances F2% and F10% in Figure 11, we see a large change in the equivalent widths and peak excess absorptions.For example, model-set F2% (F10%) produces an equivalent width and peak excess absorptions of 0.013 nm (0.058 nm) and 6.6% (31.4%) at an age of 150 Myr, while at the later age of 5 000 Myr, the corresponding values are 0.0014 nm (0.007 nm) and 1.5% (7.4%).Although these models exhibited similar atmospheric escape hydrodynamics (section 4.1), our synthetic transit observations are very sensitive to the adopted helium abundance.
The F10% and F10%PostProc model sets show increased 1083 nm absorption going from 43 to 150 Myr while F2% maintains roughly the same absorption between these two ages, as discussed in the previous subsection.The F2% model set exhibits a similar evolution of atmospheric escape to the F10% model, however unsurprisingly the latter has overall more He(2 3 ) extinction and accordingly greater 1083 nm absorption.The larger ratio of obscuring He(2 3 ) material relative to the stellar disk for F10% compared to the F2% model set emphasises the observational effect of the more drastic stellar disk area variation coinciding with a lesser variation in atmospheric escape between 43 and 150 Myr, leading to the bump in Figure 11 at 150 Myr for F10% and F10%PostProc models.

Effect of helium energetics on 1083 nm observability
Figure 11 shows similar transit properties for models F10% and F10%PostProc, with F10%PostProc having smaller and narrower absorptions.Interestingly, the inclusion of helium energetics impacts the resulting 1083 nm signature, however, as for the modelled hydrodynamics (section 4.2), the resulting variation is small over the modelled parameter space.
When excluding helium energetics, the missed heating from helium photoionisation and hydrogen photoionisation by helium produced photons is partially balanced by a larger availability of (stellar) photons capable of photoionising hydrogen, contributing to atmospheric heating (the F10%PostProc sets all photoionisation weighting factors to 1).100% of the received stellar photons (bar those in the mid-UV due to their low energy) are then available to photoionise hydrogen.For a particular photoionisation, there is greater excess kinetic energy to heat the atmosphere than if the same photon was allowed to ionise He(1 1 ).At older ages, the variation between the self-consistent and post-processing model is smaller on account of the reduced level of photoionisations.Transit phases of -0.5 and +0.5 correspond to the time of first and fourth contacts, respectively.The mean average of all spectra between these two phases is shown in dashed black.The upperpanel shows the spectra at a planetary age of 16 Myr while the lower-panel offers a closer look at the spectra at age 5 000 Myr.The assumed planetary parameters are that of our model-set entitled F2% (see Table 4).

Summary of main findings
In this paper, we investigated the evolution of atmospheric escape for highly irradiated 0.3 Jup exoplanets with primordial hydrogen/helium atmospheres orbiting a K-dwarf star at 0.045 au, with the goal of predicting the evolution of their helium 1083 nm transits.Our model self-consistently solves the fluid dynamic equations in addition to the coupled equations for hydrogen ionisation balance and the helium 1 1 , 2 1 , 2 3  and He + states.For the transmission  4).spectroscopy modelling of the helium triplet signature at 1083 nm, we use a ray tracing technique previously applied to Lyman- and H- transmission spectroscopy modelling (Allan & Vidotto 2019).We utilise evolving predictions of stellar flux and radius (Johnstone et al. 2021) and planetary radius (Fortney & Nettelmann 2010) as input for our evolution modeling.We also explored how the He/H abundances and self-consistent inclusion and omission of helium energetics affect both the atmospheric escape and its 1083 nm absorption signature.
Our main conclusions regarding the hydrodynamics of escaping atmosphere are listed below: • The heating processes which dominate are: photoionisation of neutral hydrogen by sEUV (36-92 nm), hEUV (10-36 nm) photons as well as by the 24.6 eV photon produced in the direct recombination to the helium 1 1  state.Although not dominant, heating arrising from the photoionisation of He(1 1 S) by hEUV photons is not negligible for abundances of 2% (Figure 5), and is more important for higher abundances (Figure B1).
• Neglecting helium heating contributions leads to a reduction of 21 and 45% in the net heating contribution at young and old ages, respectively, assuming a helium abundance of 0.1.See Figure 5.  4).
• The majority of the heating affects the sub-sonic region of the atmosphere meaning that heating reductions with planetary evolution affects the mass-loss rate more so than the terminal velocity of the escaping atmosphere.
• Increasing the helium abundance from 2% to 10% results in a slightly slower outflow, with a smaller mass-loss rate (Figure 4).This does not strongly affect the hydrodynamics, but is important in the observability of the triplet (as we will discuss below).
• Omitting helium energetics in the hydrodynamics models leads to a lower mass loss and terminal velocity (Figure 4).
Regarding the population of the helium triplet and its observability, we found the following main conclusions: • The dominant populating process for the helium triplet state is the recombination of ionised helium into the triplet state (Figure 6).
• The dominant depopulating process for the helium triplet state helium is de-excitation into the helium singlet state through a collision with a free electron (Figure 6).
• Although increasing the helium abundance from 2% to 10% does not affect the hydrodynamics, it substantially changes the observability of the transiting atmosphere.At an age of 5 Gyr, the (T1-T4 phase-averaged) peak excess absorption increases from 1.5% to 7.4% (see Figure 11).
When accounting for the evolution of the system, we showed that the diminishing XUV flux required to heat the planetary atmosphere combined with the growing gravitational force due to the shrinking planetary radii leads to the weakening of atmospheric escape as the planet evolves.As a consequence of a slower, less dense and less extended wind, the helium 1083 nm signature weakens overall with evolution (Figures 8,10,11).
We found that the strongest helium 1083 nm absorption occurs at the youngest ages.Our models predict that a very young (≲150 Myr) gas giant closely orbiting a K-dwarf star could produce a helium 1083 nm absorption of ∼4% to ∼7% assuming a helium abundance of 2% (if orbiting a slow or fast rotating star).This weakening with evolution, however, is not necessarily monotonic.For example, the transit depths and equivalent width of the planet at age 150 Myr are larger than those at 40 Myr for models with He/H= 0.1.This is due to a large drop in stellar radius corresponding with small decline of atmospheric escape.

Do we observe larger helium absorption at younger ages?
Our model suggests stronger helium 1083 nm absorption at younger ages.Is this higher absorption seen in the observations?So far, there has been a limited number of helium 1083 nm transit observations of young systems.One difficulty is the greater stellar variability in helium 1083 nm at younger ages.Another difficulty is the small number of exoplanets known to transit a young (< 500 Myr) star, and in particular a K-type star best suited for 1083 nm detections.This currently small sample consists of super-Earths and mini-Neptunes, whose escaping atmospheres are intrinsically harder to detect than for larger planets.
In recent years, helium transit observations of a few young systems resulted in either non-detections or tentative detections possibly of stellar rather than planetary origin, namely of AU Mic b (Hirano et al. 2020), of the V1298 Tau planets (Gaidos et al. 2022;Vissapragada et al. 2021), of TOI 1807band 2076b(Gaidos et al. 2023), and of the 400-Myr old mini-Neptune HD 63433c (Zhang et al. 2022c).Of these, only TOI 1807b and 2076b are thought to have K dwarf hosts, but these planets are rocky planets whose primary H-rich atmosphere might have already been lost.Zhang et al. (2022a,b) recently reported successful helium 1083 nm detections for four young mini-Neptunes orbiting K dwarfs, with mean excess absorptions of ∼1% in each.
It should be noted that the larger absorption predicted at young ages in our evolution models is not at odds with the current observations, as the current sample of observed young exoplanets orbiting K-dwarfs for which helium 1083 nm transmission spectroscopy has been performed does not yet include a gas giant planet.Hence, finding a system in which a larger radius 1 − 2 Jup close-in (< 0.1 au) exoplanet transits a young (< 150 Myr), ideally K dwarf star would be a better test to our evolution models.In a forthcoming work, we apply our modelling techniques for atmospheric escape and helium 1083 nm transmission spectroscopy to smaller, mini-Netpunes, more consistent with the current helium 1083 nm detections of young exoplanets.
//www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Figure 1 .
Figure 1.Schematic displaying the various helium processes considered in our modelling.See Table2for the rate equations and their relevant references.
, we normalised the XUV flux tracks of Johnstone et al. (2015a) so that the solar XUV flux was reproduced at the solar age.Similarly, we now normalise the Johnstone et al. (2021) XUV flux tracks so that they reproduce the XUV flux of spectral type K6 star HD 85512, utilising a spectrum obtained from the Measurements of the Ultraviolet Spectral Characteristics of Low-mass Exoplanetary Systems (MUSCLES) Treasury Survey (France et al. 2016), shown in Figure 2. MUSCLES combines Hubble and ground-based observations

Figure 2 .
Figure 2. [Upper-panel] MUSCLES spectrum of HD 85512 used in to normalise the high-energy evolution tracks of Johnstone et al. (2021), shown in Figure3.The energy of the three considered photons produced due to helium transitions in the planetary atmosphere are marked by the dashed lines.[Lower-panel]Photo-ionisation cross sections for hydrogen and considered helium states.Note that while He(2 1 ) photoionisations are not considered in the model, its cross-section profile is shown here for comparison to that of the 2 1  and 2 3  helium states.The four shaded regions indicate the wavelength channels of X-ray, hEUV, sEUV and mid-UV going from left to right.The representative energies of each of these flux bins is marked by the vertical solid white lines in the lower-panel.

Figure 3 .
Figure 3. [Upper-panel] Evolution of high-energy stellar luminosity, where hEUV, sEUV and X-ray wavelength bins were obtained by normalising the predictions of Johnstone et al. (2021) for a 0.7- ⊙ by the K dwarf HD85512 (star symbols, from Figure 2).The mid-UV luminosity tracks combines the near-UV and far-UV fluxes from Richey-Yowell et al. (2022).[Central-panel] Stellar radius evolution for the same star, also obtained from the model of Johnstone et al. (2021).[Lower-panel] Planetary radius with respect to age for a 0.3- jup planet orbiting a solar-like star at 0.045 au(Fortney & Nettelmann 2010).The vertical grey lines indicate the sampled ages for our various classes of models.

Figure 4 .
Figure 4. Mass-loss rate, cumulative mass lost and the terminal velocity of the escaping atmosphere as a function of planetary age are shown in the upper, central and lower-panels, respectively.The evolution model sets are described in Table4. C

Figure 6 .
Figure 6.Rates of processes directly populating (solid) and depopulating (dashed) the helium triplet state at young (left) and old (right-panel) ages as a function of distance.Table 2 lists each considered transition.Here, only rates within the displayed −axis range are listed in the Figure legend.The assumed planetary parameters are that of our F2% model (Table4).

Figure 7 .
Figure 7. Temperature (left axis, solid line-style) and various number density profiles (right axis, dashed and dotted for helium and hydrogen, respectively).The upper-panel corresponds to an age of 16 Myr while the lower-panel shows 5 000 Myr.The planetary parameters are that of our F2% model (see Table4).Similar figures for our F10% and F10%PostProc are given in the Appendix.

Figure 8 .
Figure 8. Plane-of-the-sky extinction contributions at 16 (upper-row) and 5 000 Myr (lower-row), shown at mid transit.Each column shows the contribution of one of the three lines of the helium triplet.The dashed orange circle marks the stellar disk with respective radii of 0.89 and 0.66  ⊙ at ages 16 and 5 000 Myr.The planetary parameters are those of our F2% model (see Table4).

Figure 9 .
Figure9.Helium 1083 nm transmission spectra as a function of time/transit phase as indicated by the colour bar.Transit phases of -0.5 and +0.5 correspond to the time of first and fourth contacts, respectively.The mean average of all spectra between these two phases is shown in dashed black.The upperpanel shows the spectra at a planetary age of 16 Myr while the lower-panel offers a closer look at the spectra at age 5 000 Myr.The assumed planetary parameters are that of our model-set entitled F2% (see Table4).

Figure 10 .
Figure10.[Upper-panel]The helium 1083 nm transmission spectra averaged over phases between first and fourth contacts.The line-style relates each spectrum to its planetary age.[Lower-panel] A zoom-in of the grey shaded region in the upper-panel.The 5 000 Myr transmission spectrum is shown by the solid line as well as its individual line contributions as shown by the legend.The planetary parameters are those of our F2% model (see Table4).

Figure 11 .
Figure 11.T1-T4 phase-averaged helium 1083 nm equivalent widths integrated over 1082.6-1083.35nm (upper-panel) and peak excess absorptions (lower-panel) as a function of planetary evolution, for our model sets indicated in the legend (see Table4).

Figure B1 .
Figure B1.The same as the upper-panel of Figure 5, now showing individual heating components for our F10% (upper-panel) and F10%PostProc (lower-panel) model.