The XUV-driven escape of the planets around TOI-431&$\nu^2$ Lupi

One of the leading mechanisms invoked to explain the existence of the radius valley is atmospheric mass loss driven by X-ray and extreme-ultraviolet irradiation, with this process stripping the primordial envelopes of young, small planets to produce the observed bimodal distribution. We present an investigation into the TOI-431 and $\nu^2$ Lupi planetary systems, both of which host planets either side of the radius valley, to determine if their architectures are consistent with evolution by the XUV mechanism. With $\textit{XMM-Newton}$, we measure the current X-ray flux of each star, and see evidence for a stellar flare in the TOI-431 observations. We then simulate the evolution of all of the transiting planets across the two systems in response to the high-energy irradiation over their lifetimes. We use the measured X-ray fluxes as an anchor point for the XUV time evolution in our simulations, and employ several different models of estimating mass loss rates. While the simulations for TOI-431b encountered a problem with the initial calculated radii, we estimate a likely short ($\sim$ Myr) timespan for primordial envelope removal using reasonable assumptions for the initial planet. $\nu^2$ Lupi b is likely harder to strip, but is achieved in a moderate fraction of our simulations. None of our simulations stripped any of the lower density planets of their envelope, in line with prediction. We conclude that both systems are consistent with expectations for generation of the radius valley through XUV photoevaporation.


INTRODUCTION
The radius(-period) valley is a dearth of planets with radii of 1.5 to 2 R ⊕ at short orbital periods up to a few tens of days.First predicted almost a decade ago (Owen & Wu 2013;Lopez & Fortney 2013), it was discovered observationally by the California-Kepler Survey following their reassessment of Kepler system properties with improved stellar parameters (Fulton et al. 2017).Further analysis revealed a negative slope as a function of period (e.g.Martinez et al. 2019), including that performed with a smaller sample of systems with more precise characterisation (Van Eylen et al. 2018).The two populations of planets straddling the valley are thought to be the smallest planets able to retain a gaseous envelope (those just above the valley), and bare rock with at most a small atmosphere of heavier elements (those below the valley).
The two leading proposed mechanisms for producing the valley are core-powered mass loss and photoevaporative escape driven by X-ray/extreme-ultraviolet (EUV; together, XUV) light.Both of these produce the bare rock population via the stripping of primordial gaseous envelopes.In the core-powered mass loss mechanism, energy stored in the planet from its formation gradually radiates out into the envelope as the rocky core cools.The thermal energy supplied to the atmosphere in this way has been shown to be sufficient to ★ E-mail: kinggw@umich.edudrive substantial mass loss over the first Gyr or so of a planet's lifetime (Ginzburg et al. 2016;Gupta & Schlichting 2019, 2020;Gupta et al. 2022).These studies can reproduce both the bimodality and negative slope in radius-period space that are observed in the exoplanet population.Alternative explanations for the valley such as different core compositions for the two populations (Venturini et al. 2020), and a primordial generation at the epoch of gas accretion (Lee & Connors 2021;Lee et al. 2022) have also been proposed, and warrant further study.However, in this work, we focus on the XUV-driven escape mechanism.
In the XUV-driven escape mechanism, XUV photons are absorbed high in the atmosphere of the planet, efficiently heating this region, which sets up and drives a hydrodynamic outflow.The XUV irradiation is highest in the early life of the system and reduces in intensity as the star undergoes magnetic braking (e.g.Micela et al. 1985;Güdel et al. 1997;Micela 2002;Feigelson et al. 2004;Jackson et al. 2012;Johnstone et al. 2021), so one expects the resulting mass loss to also reduce as the system ages.The first 100 Myr, when the XUV irradiation is at its highest level, has therefore been widely thought to be the dominant epoch for evolution due to XUV photons, and the timescale on which the radius valley is imprinted into the population (e.g.Lopez & Fortney 2013;Lammer et al. 2014;Owen & Wu 2017).This would provide a possible method of distinguishing which mechanism, if either, dominates the generation of the radius valley, if one had a large enough population of exoplanets at all young ages Table 1.Adopted stellar parameters for the two systems we investigated.Parameters are taken from Osborn et al. (2021) and Delrez et al. (2021) for TOI-431 and  2 Lupi, respectively, except for distances which are from Gaia DR3 parallaxes (Gaia Collaboration et al. 2022) up to 1 Gyr.However, King & Wheatley (2021) showed that the reduction in EUV irradiation may be significantly slower than X-rays.This in turn could prolong the timescale on which sub Neptune-sized planets are significantly evolved via XUV-driven escape up towards that expected for core-powered mass loss.Other methods of trying to distinguish between the mechanisms have also been proposed, such as examining the radius valley in a 3D parameter space (Rogers et al. 2021).
The NASA TESS mission has discovered many multi-planet systems over the past few years which contain at least one planet just above the radius valley, and at least one planet just below it (e.g.Günther et al. 2019;Cloutier et al. 2020;Demory et al. 2020;Hawthorn et al. 2022).Such systems are ideal laboratories to test the mechanisms theorised to explain the valley by examining whether they have properties consistent with a particular mechanism(s).In the case of photoevaporation, one advantage of having the planets in the same system is that they will have experienced an identically evolving XUV irradiation history (e.g.Owen & Campos Estrada 2020).
Among the TESS discoveries of systems with planets either side of the radius valley are TOI-431 and  2 Lupi (HD 136352), both of which host three known planets.Two of the TOI-431 planets have been observed to transit: b and d (Osborn et al. 2021).The system is unusual in that planet d transits despite being exterior to the nontransiting planet c.Planet b orbits on a 0.5 d period, and its density is consistent with it being a rocky core with no envelope.Planet c has a minimum mass similar to b, while planet d has a density consistent with having a gaseous envelope.The three planets known to orbit the naked-eye brightness  2 Lupi star were originally discovered via radial velocities (Mayor et al. 2011;Udry et al. 2019).Transits of planets b and c were later detected by TESS (Kane et al. 2020), and d was also found to transit by CHEOPS (Delrez et al. 2021).Planet b is consistent with being a rocky core, while planets c and d have densities that suggest they have a gaseous envelope.The visual brightness of  2 Lupi means that the system is an excellent target for atmospheric characterisation of its planets with transmission spectroscopy.Furthermore, with both systems residing relatively nearby to the Solar System, they are possible future targets for Ly- or helium triplet observations targeting escaping atmospheres.In our work, we adopt parameters for the TOI-431 system from Osborn et al. (2021), and for the  2 Lupi system from Delrez et al. (2021).We display some key characteristics for the stars and planets in Tables 1 and 2, respectively.
The aim of this work is to test whether the TOI-431 and  2 Lupi systems are consistent with the XUV-driven photoevaporation mechanism for generating the radius gap.We measured the X-ray emission from the two stars with XMM-Newton, which is important for two reasons.First, knowledge of the XUV irradiation is an important input to interpreting observations attempting to detect extended or escaping atmospheres of the planets in these systems.Second, direct measurement is important for systems of interest as empirical relations that can be used to estimate X-ray emission are subject to variation up to an order of magnitude (e.g.Wright et al. 2011Wright et al. , 2018;;Jackson et al. 2012).Following analysis of our observations, we ran simulations for the lifetime evolution of the transiting planets in the systems under the action of XUV-driven photoevaporation, using the measured fluxes as an anchor point for the escape-driving XUV emission evolution.Using the simulations, we assessed whether XUV irradiation could have sculpted the systems into what we observe today.

OBSERVATIONS
We observed TOI-431 and  2 Lupi with XMM-Newton on 2021 Aug 28 and 2021 Sep 12, respectively (OBSIDs 0884680101 and 0884680201; PI: King).In both cases, the prime instrument was the EPIC-pn camera.
We operated the simultaneous Optical Monitor (OM) observations in fast mode with a single near-ultraviolet (NUV) filter, in order to monitor the stars for short term variation, such as flares.The TOI-431 observations were made with the UVW1 filter (effective wavelength 291 nm, width 83 nm), while the  2 Lupi observations were made with UVW2 (effective wavelength 212 nm, width 50 nm).
Both targets were successfully detected in the EPIC-pn observations.We analysed the observations using the Scientific Analysis System (sas 19.1.0),using the standard process, as detailed on the "SAS Threads" webpage1 .Elevated background due to Solar protons, a problem well-known to affect XMM-Newton observations (Walsh et al. 2014), affected small portions of both observations in the EPICpn camera.As we have done previously (e.g.King et al. 2018), we double the threshold suggested in the SAS Threads for cutting events due to high background, as this results in less cuts, and therefore more usable counts, without adversely affecting the results.For TOI-431, there were two high background peaks in the final quarter of the observation, while the  2 Lupi observations contain a single peak right at the beginning of the observation.Both EPIC-MOS cameras were unaffected by this issue throughout both observations.
In our data reduction of the OM data, we used the same processes as described in King et al. (2018): following running the standard chains for image and fast mode, we corrected the fast mode time series using the image mode count rates.Both observations exhibit relatively high count rates of over 100 s −1 , though both are well within the valid range for coincidence loss correction, given the frame times of 7.5 ms and 5.5 ms for TOI-431 and  2 Lupi, respectively.A description of the brightness limits of the OM can be found in section 3.5.5 of the XMM-Newton Users Handbook 2 ; the raw fast mode count rates (i.e.before corrections for e.g.coincidence loss) for comparison with table 21 and the surrounding descriptions are 60 s −1 and 70 s −1 for TOI-431 and  2 Lupi, respectively.
Both sets of OM images exhibit commonly seen artifacts.For TOI-431, there was a faint, extended loop of emission extending from the bottom right of the location of the target.This is caused by a chamfer in the housing of the detector reflecting light from a bright source just outside the detector field of view (Mason et al. 2001).The overall contamination of the target from this effect is low.The number of counts per pixel in the loop is much less than one per cent that in the central region of the target's point spread function (PSF), and the majority of the loop does not overlap the source.For  2 Lupi, there is a faint signature of the detector's modulo-8 pattern in the source PSF (Mason et al. 2001).The algorithm that removes this signature from the background of OM images is imperfect for bright sources due to coincidence losses, though photometery is unaffected due to the number of photons being conserved (Breeveld et al. 2010).

X-ray and NUV light curves
In Figs. 1 we show the observed XMM-Newton light curves for TOI-431, and for  2 Lupi in Fig. 2. The X-ray light curves depict the coadded count rates across the three EPIC cameras -pn, MOS1, and MOS2 -though for both targets the observed counts are dominated by EPIC-pn.Also included are the NUV light curves plotted from the data taken with the OM.
The TOI-431 X-ray light curve in Fig. 1 is split across the top two panels, separating the soft (0.166-0.65 keV) and hard (0.65-2.4 keV) energies.The orbital phases plotted are for planet b, where midtransit is at a phase of 1.Though the observation covers a primary transit (grey shaded region), there is no transit visible in the X-ray light curve.This is unsurprising given the small size of planet b and the count rate errors at this cadence of 10-20 per cent.The phases covered by the observation for the other transiting planet d are 0.437-0.466.
The light curves in both energy ranges show evidence of stellar X-ray variation.During the last quarter of the observation, the count rate increases by about 50 per cent on the previous quarter in both the soft and hard bands.While there are two peaks in the high-energy background due to Solar protons (discussed in Section 2) around 24 and 30 ks into the observation, these are short, lasting less than a where  and  are the soft and hard band count rates, respectively.The resulting plot shows no evidence of variation in the hardness ratio throughout the observation, including during the possible flaring epoch.The fact that the increase is similar in both energy bands point to the variation possibly not being the result of a flare on the host star, since stellar flares typically shows a hardening of the Xray emission (e.g.NUV (UVW2) Figure 2. X-ray and NUV light curves for  2 Lupi.The X-ray data are plotted in the top panel, binned to a cadence of 2 ks, and for the full 0.166-2.4keV energy band.They are the total from coadding all three EPIC cameras together.
In the bottom panel is the NUV data measured with the OM, using the UVW2 filter.
2015).However, the NUV light curve in the bottom panel of Fig. 1 also shows an increase in count rate at the end of the observation.In this case, the rise is 2 s −1 (1.7 per cent), and appears to peak some 50 minutes after the start of the bump in X-rays, though there is some overlap between the increases.The two highest points in the X-ray light curve cover the phase range 0.947 -1.042, while including the point either of this extends the range to 0.900 to 1.089.
In comparison, the phase range covered by the last three points in the OM light curve constituting the rise is 1.008 -1.080.Additionally, delays in the peak of flares between different wavelengths on the order of tens of minutes would not be unprecedented (e.g.MacGregor et al. 2021).It is unfortunate that the OM observation ended when it did, because observations of the following hour could have allowed us to unambiguously identify whether or not it was a true stellar flare by looking for the characteristic flare decay shape.In its absence, and in the absence of a hardening of the X-ray emission, we can only conclude the ramp up to possibly be due to a stellar flare, with other scenarios possible too.One possibility is the rotation into view of a particularly active region of star.In that instance, differing optical thickness between X-ray and NUV could explain the delay in the increase between the two wavelengths.We do not split the  2 Lupi X-ray light curve in the top panel of Fig. 2 into separate soft and hand bands, as there were too few hard energy photons through the observation for it to be useful on its own.The observation did not cover a transit of any of the three planets.Planet b's orbital phases during the observation are plotted on Fig. 2. For planets c and d, the orbital phases were 0.6841-0.6925and 0.2908-0.2930,respectively.
As for TOI-431, the  2 Lupi light curve also shows evidence of stellar variation.The most notable of this is in the last 3 ks of the observation, where the count rate jumps up from being consistent with zero to over 12.5 ks −1 .This could be the beginning of a stellar flare, but again, without capturing the full flare in the observation, it is difficult to speculate further.With an age of 12.3 Gyr, one might not expect to observe many flares from this star, as flare occurrence decreases with age (e.g Skumanich 1986;Davenport et al. 2019).
Additionally, the NUV light curve with the OM shows no rise at this point in the observation, and instead appears to decrease slightly.Overall, the OM light curve is largely unremarkable and exhibits little variability.

X-ray Spectra
We fit the EPIC-pn observations in Xspec 12.11.1.We did not include the EPIC-MOS data in the fits for either star because the signal was small compared to the EPIC-pn data, and so its inclusion did not improve our model fits.
The TOI-431 data was binned to a minimum of 25 counts per bin, while the  2 Lupi data was binned to 10 counts per bin due to the lower number of overall counts.We fit the spectra with APEC models, which describe an optically thin, collisionally ionised plasma (Smith et al. 2001).For TOI-431, we had to use three temperature components to obtain a good fit, whereas for  2 Lupi a single temperature was sufficient due to the much fewer spectral counts.We also included a TBABS term to account for interstellar absorption (Wilms Note that we use the following energy band definitions: X-ray: 0.2-2.4keV EUV: 0.0136-0.2keV et al. 2000), fixing the hydrogen column density to an assumed local density of 0.1 cm −3 (Redfield & Linsky 2000).Abundances were fixed at Solar values, according to Asplund et al. (2009).We used Cash statistics in our fitting (Cash 1979), and used Xspec's built-in MCMC sampler in order to assess the 1- uncertainties on the model parameters and fluxes.
Figures 3 and 4 show our extracted spectra and best fitting models for TOI-431 and  2 Lupi, respectively, while in Table 3 we give our best fit model temperatures, , emission measures, EM, and unabsorbed fluxes at Earth,  X .Using our fluxes, we compute the X-ray luminosity,  X , of the stars, and combine this with the X-ray-EUV scaling relations from King et al. (2018) to estimate the EUV luminosities,  EUV .Finally, we estimate the current XUV irradiation rate at each of the planets in the two systems by scaling to their semimajor axes.
Both spectra are dominated by softer energies, especially in the case of the old star  2 Lupi for which the detected emission above 0.6 keV is consistent with zero.This is not so surprising since the coronae of older stars have been seen to exhibit much lower emission above this energy compared to younger stars (e.g.Güdel et al. 1997), and given the observed relationship between the emission measure weighted average coronal temperature and X-ray surface flux (which in turn decreases with age) (Johnstone & Güdel 2015).For TOI-431 however, there is a small bump in the spectrum between 0.5 and 0.9 keV, characteristic of Fe L-shell emission.Unlike for much younger stars where this bump can dominate, the overall spectrum is still dominated by the soft end, suggestive of a star that is ≳ 1 Gyr in age.This is consistent with the previously inferred age of 1.9±0.3Gyr (Osborn et al. 2021), estimated from its chromospheric activity data and the Mamajek & Hillenbrand (2008) empirical relations.

ATMOSPHERIC EVOLUTION SIMULATIONS
Using our measured fluxes as an anchor point, we ran a series of simulations to investigate the evolutionary histories of the five transiting planets in the two systems in our sample.We do not consider TOI-431 c as it is not transiting, and as such we have no measurement of the current radius of the planet.For the planetary evolution, we employed the 1-D stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019;;Jermyn et al. 2023) version 12778.The code we implemented is exactly that used by Malsky & Rogers (2020) and Malsky et al. (2023), based on Chen & Rogers (2016), except for the final step in which the planet is evolved in time with XUV irradiation.We edited this final step to change the mass loss rate prescription, for which we use three different methods described below, and the time evolution of the XUV emission from the star, which we calculated using the Mors code (Johnstone et al. 2021).
The evolutionary tracks calculated by the Mors code provide a more complex XUV luminosity progression than assuming a constant saturated period followed by a power law decline (e.g.Jackson et al. 2012).For stars a few hundred Myr old, there is evidence of increasing dispersion in the observed  X / bol values, which later reconverge at an age of about a Gyr (see also e.g.Tu et al. 2015), something which Mors aims to model.We generated an evolution track for each star based on the stellar mass, then scaled the median track to the measured  XUV from Section 3.2.We note that the measured  XUV for  2 Lupi was a factor of three smaller than the median track calculated by Mors.This could be related to the relatively uncertain age of the star, it being in a quiet activity state during our observations, or just a demonstration of the scatter in the XUV time evolution relations.We also scaled the Mors  bol track for both stars to their measured values, although the effect on the simulations from this correction is very small.
We adopt the same atmospheric boundary conditions as Malsky et al. (2023): a grey Eddington T() relation, and model the atmosphere up to an optical depth of 2/3.Following Malsky & Rogers (2020), we calculate the incident bolometric flux and the column depth of irradiation for each planet.Furthermore, we use the gaseous mean opacity tables from Freedman et al. (2014).
In King et al. (2022), we predicted the mass loss rates and future evolutionary tracks for four transiting planets in the 670 Myr old Praesepe open cluster.This time, we focus on reconstructing the past mass loss history for older field age systems.This meant defining a grid of starting masses and envelope mass fractions, displayed for each planet in Table 4.The grid choices were based on the current known mass characteristics of the planets.The range of simulated masses are wider for the two planets consistent with being bare rock at the present age, TOI-431 b and  2 Lupi b, since the aim of those simulations are somewhat different, as explained in Section 4.2.
We employed a few different methods of estimating the mass loss rates as part of our simulations, in order to compare the methods.The first of these is energy-limited mass loss, first formulated in the context of the Solar System planets by Watson et al. (1981), and which has since seen widespread use for exoplanets over the past two decades (e.g.Lecavelier Des Etangs 2007; Sanz-Forcada et al. 2011;Salz et al. 2015;Louden et al. 2017;Poppenhaeger et al. 2021).The energy-limited mass loss rate,  En , is given by (Baraffe et al. 2004;Erkaev et al. 2007) where  is the heating efficiency,  is a correction factor to the distance above the optically-measured radius of the planet at which XUV photons are on average absorbed, and  is a Roche-lobe correction factor (Erkaev et al. 2007).We assume a canonical  value of 0.15 (e.g. Watson et al. 1981;Salz et al. 2015;Kubyshkina et al. 2018a), and a  value of 1.The latter is in reality a lower limit, as XUV photons will be typically absorbed higher up in the planet's atmosphere.
The second method of estimating mass loss we explored was the Ates photoionisation hydrodynamics code (Caldiroli et al. 2021).Performing the full hydrodynamic calculation for each step on each track would not have been computationally possible.We do, however, run Ates in order to obtain current mass loss rates for the planets in Section 4.3.To simulate the grid of mass loss histories within a reasonable timeframe, we instead use an approximation to their code provided in a follow up paper (Caldiroli et al. 2022).The approximation gives both the mass loss rate and effective efficiency function,  eff , based on fitting an analytic expression to the results of the Ates code.Note that  eff is not the same at the heating efficiency  defined above, as  eff also encapsulates the effects of the unknown  parameter.
The final method was based upon the hydrodynamic modelling of Kubyshkina et al. (2018a).They produced a grid of models to estimate the mass loss rate, given some star and planet characteristics.As with ATES, the interpolator they provide in Kubyshkina & Fossati (2021) was too computationally expensive to couple with MESA, and so we instead use their "hydro-based approximation."This was detailed was described in Kubyshkina et al. (2018b), with the relation we used given in their equation 9.This is an analytic expression fitted to values they obtained from their grid of models.
We began our evolution models at 5 Myr old, an age by which the majority of protoplanetary discs have dissipated (e.g.Haisch et al. 2001;Mamajek 2009;Ribas et al. 2015).We do not directly consider the "boil-off" phase in our models, which is a more rapid process thought to occur on kyr to Myr time-scales immediately after the disc dissipates, and before the onset of meaningful XUV photoevaporation (e.g.Stökl et al. 2015;Owen & Wu 2016;Ginzburg et al. 2016;Fossati et al. 2017;Rogers et al. 2023).During this phase, the accumulated planetary atmosphere can contract significantly.The release of binding energy during contraction, together with the atmosphere no longer being protected from the bolometric stellar radiation by the disc, can drive substantial mass loss in this short epoch.We note that the Kubyshkina et al. (2018a) model does partially cover this effect, and so our results based on that work could be affected by this phase of evolution.

Results for Gaseous Planets
Three of the planets in our sample have densities consistent with retaining a substantial gaseous envelope at their current age: TOI-431 d,  2 Lupi c, and  2 Lupi d.We assessed which of the mass loss tracks, described above, resulted in radii and masses that are consistent with the currently measured mass and radius values to within 1-.Starting Envelope Mass Fraction No. of Methods Consistent  2018a) models (Kubyshkina et al. 2018b).Some simulation tracks lie outside the limits of the plot, but none of those produced any consistent tracks.

TOI-431 d
In Fig. 5 we show a heatmap of the starting grid of masses and envelope mass fractions for TOI-431 d, highlighting which of those initial conditions led to consistent results for that planet.We label each mass loss rate calculation method, and assign a shade displaying how many of the methods worked.In Fig. 6, we show the evolution of TOI-431 d's key characteristics over time for these "successful" evolutions, as well as the evolution of the high-energy emission from the host star.
The results demonstrate we can successfully reproduce TOI-431 d with XUV-driven escape, and with all three mass loss methods used.The results do however differ somewhat between the methods.From Fig. 6, we can see that the energy-limited mass loss rates for the consistent tracks are far lower than the other two methods.This is not to say that there were not other simulations for which the energylimited method gave higher mass loss rates, or simulations with the other methods which did not give lower rates, but none of them were able to produce TOI-431 d as we see it today.This is also manifested in the heatmap in Fig. 5, where the successful starting envelope mass fractions for the energy-limited method are below all of the successful values for the other two methods.
The hydro-based approximation is perhaps at the other extreme.A wide range of starting envelope mass fractions can end up with a TOI-431 d-like planet at the current age.Looking at Fig. 6, it seems that in those cases with much higher starting values, the envelope mass fraction quickly decreases and all of the curves converge to a value less than 10 per cent.The calculated mass loss rates that seem to be driving this, as with the energy-limited conclusions above.The early mass loss rates for some tracks are an order of magnitude greater than for ATES, and two orders of magnitude greater than energylimited, allowing for much greater atmospheric stripping than the other methods.
To explore this further, in Fig. 7 we plot the first 500 Myr of evolution of the mass loss rate and the restricted Jeans escape parameter, Λ, from equation 2 of Kubyshkina et al. (2018b) for both ATES and hydro-based approximation.For the most direct comparison possi- ble, only the four tracks where both methods reproduced the observed planet mass and radius.This plot demonstrates that the hydro-based approximation mass loss rates are higher than ATES for these planets for about the first 100 Myr.Some of this could be due to the consideration of the effect of boil-off in the Kubyshkina et al. (2018a) models.
The initial Λ values of 20 lie at the upper limit of where Kubyshkina et al. (2018a) expect boil off to be relevant.However, this elevated mass loss rate lasts longer than the expected ∼Myr timescale on which boil-off is thought to occur (e.g.Owen & Wu 2016), and the reason why the rates remains elevated above ATES for much longer is less clear.We also note Fig. 6 shows that by a few Gyr in age the mass loss rate for the successful hydro-based approximation tracks are actually all below the successful ATES tracks.
Again, all this is not to say that the other methods did not have simulations with such high mass loss rates early on too, but it seems none of them were able to produce planets with radii and mass consistent with measurements at the current age.The fact that the hydro-based approximation is such an outlier suggests that dramatic XUV-driven atmospheric loss, removing up to 25 per cent of the planet's starting mass by the current age of 1.9 Gyr, is perhaps in reality unlikely.

𝜈 2 Lupi c & d
For the other two low density planets in our sample,  2 Lupi c and d, none of our simulations reproduced the radius at the current age.However, inspection of the results revealed numerous tracks did not reach the system age.For these tracks, MESA was unable to find an acceptable model and so ended the evolution early, often a few Gyr before the 12.3 Gyr age estimated for the system.Upon investigation, we found that the likely explanation was the evolutionary models reaching the boundaries of what MESA is capable of simulating.Mass Loss Rate (g/s) Figure 7.The first 500 Myr of evolution of the mass loss rate and the restricted Jeans escape parameter, Λ (equation 2 of Kubyshkina et al. 2018b), for the four tracks where both ATES and the hydro-based approximation both reproduce the observed planetary mass and radius within 1- at the age of the system.
The SCVH tables used by MESA in calculating the equation of state at the low densities and temperatures required for these planets hold for log  = log  − 2 log  + 12 ≤ 5.0 (Saumon et al. 1995).The tracks that ended early were found to have log  max > 5 in their final few steps before terminating.The early end of some simulations can also be seen in some of the consistent results for TOI-431 d plotted in 6, where many of the tracks also end before the 10 Gyr time limit for that planet's simulations.The difference is that TOI-431 is estimated to be much younger, and so we can perform the consistency test at 1.9 Gyr before the simulations fail later at ≳ 8 Gyr.
The key takeaway from the simulations of these two planets is that none of the evolutions stripped the planet of its envelope, or was on track to do so.The lack of consistency among those tracks that reached the system age arose because the planet was too large, indicative of an envelope that was too thick, not too thin as to be stripped or nearly so.Moreover, the lowest envelope mass fraction reached in any simulation was 4.69 per cent for planet c, and 4.93 per cent for planet d.Both were for simulations with a starting envelope mass fraction of 0.05, indicating that little mass loss had taken place to reach these values.Given all simulations reached an age well past the early epoch of the highest mass loss rates, it is clear that none of those that ended early would have stripped the planet by the current age, or indeed by the time its star evolves off the main sequence.Even in the absence of our ability to determine consistency at the system age, this is a strong hint that these planets are consistent with the XUV photoevaporation for sculpting the valley.
We can draw a few conclusions about what starting conditions, if any, could have given consistent results for these two planets.In Figs. 8 and 9, we show heatmaps of which simulations finished before 11 Gyr.Since none of those which did successfully run to the system age gave mass and radii that were both consistent, we can rule out any of the grid positions except those highlighted in these heatmaps as successful starting conditions.Almost all of the tracks that ended early have envelope mass fractions below 0.10, in line with that noted by Malsky & Rogers (2020) when encountering the same issue with the simulation of cooler planets.Additionally, while none of the radii were consistent within 1- among the simulations that reached the system age, many of the masses were.For both planets, consistent masses were found for starting masses towards the lower end of the grid, with none higher than 12.25 M ⊕ for planet c, and none higher than 9.75 M ⊕ for planet d.The dependence of these starting mass consistency limits with mass loss method and starting envelope mass fraction were weak in both cases.This suggests that any consistent results for lower starting envelope mass fraction simulations that did not reach the system age would likely have fallen within limits.
As a further check, we also ran "hybrid" simulations for these two planets, wherein we evolved the planet for 100 Myr with MESA, and then the rest of the simulation using the empirical relation from Chen & Rogers (2016) to recalculate the planet's radius at each timestep.This method yielded consistent masses and radii at the system age for a few of the simulations for each planet: for c, four simulations all starting at an envelope mass fraction of 0.07, with masses starting between 11.0 and 11.75 M ⊕ ; for d, eight simulations all starting at an envelope mass fracton of 0.05, and starting masses between 8.0 and 9.75 M ⊕ .These starting conditions are in line with what we inferred above from MESA alone.

Results for Stripped Planets
TOI-431 b and  2 Lupi b both have small radii and high densities, and are therefore consistent with already having been stripped of any primordial envelope they may have had.In this instance, we ran the simulations over a wider grid.We explored the limits of what initial conditions would lead the planet to lose its atmosphere.In these simulations, we assumed that any atmosphere of heavier elements that remains now increases the mass and radius by a negligible amount.Therefore, we only allowed starting conditions that implied an initial core mass that is within 3- of the current measured planetary mass.

TOI-431 b
None of the simulations for TOI-431 b advanced more than a few Myr in the maximum 5000 steps we initially allowed for the evolution.The timesteps in MESA are dynamic, and so we interpreted the code forcing the them to be small being due to rapid changes to the planet over short periods of time.We further investigated by running a subsample of the original grid for 200,000 steps, which also failed to advance beyond the first few Myr after our starting age.Inspection of the tracks revealed very large initial planetary radii (> 50 R⊕) and brisk further expansion across the small elapsed time the simulations were able to cover.
The failure of our simulations for this planet are likely reflective of the extreme environment in which this planet resides, with an orbital separation of just 0.011 au, and its relatively low total mass.In the final setup stage for our models in MESA before the planet is evolved with XUV-driven escape, the planet is gradually irradiated to the level of the bolometric flux of the star, calculated from its  eff , over 6 Myr (Malsky & Rogers 2020) The small orbital separation therefore leads to this final setup stage puffing up the planet significantly under the intense bolometric irradiation, and when we attempt to then evolve the planet under XUV irradiation, it undergoes further rapid expansion.
In lieu of usable results from our simulations, we can still make some inferences if we assume some reasonable characteristics for a TOI-431 b planet.We took a starting mass of 4 M ⊕ , close to the centre of our test grid.We then also assumed a bulk density of The squares are coloured according to the number of mass loss calculation methods that suffered an early end to the simulation at that starting point in the grid, while the text labels are as in Fig. 5.We used lower case letters here, as well as a different colour bar, to distinguish from that plot however, as they display very different things.No. of Methods Ended <11 Gyr 0.5 g cm −3 (implying an initial radius of 3.5 R⊕); a conservatively high initial value when compared to the simulations of  2 Lupi b in Section 4.2.2,where the initial densities simulated by MESA are much lower.We calculated the XUV flux at the planet from the first Mors step for TOI-431 (4.4 × 10 6 erg s −1 cm −2 ), leading to an energy-limited escape mass loss rate of 2×10 14 g.This is sufficient to remove a 20 per cent initial envelope in 0.6 Myr, assuming a constant mass loss rate.Therefore, while we cannot simulate the system in MESA, this quick calculation demonstrates the ease with which a primordial envelope may have been removed, due the short orbital separation.

6 8
Starting Mass (M ⊕ ) No. of Methods Consistent

𝜈 2 Lupi b
The results for  2 Lupi b are very different to TOI-431 b.In the majority of the simulations, the MESA evolution made it to the estimated system age.We tested the simulated tracks for consistency with the measured mass and radius, but loosened the restrictions as compared to the tests for gaseous planets.Firstly, we tested the closest point in the track to the system age, regardless of if the track ended early.We did not want to skip planets stripped well before 12.3 Gyr, not least because their radius and mass have likely changed little between the point the envelope was removed entirely and the present day (unlike gaseous planets for which the track ended early with an appreciable atmosphere remaining).Second, we only required that the mass and radius be within 3- of the measured values, rather than 1.Our primary aim with these simulations is investigating whether the planet can be stripped by XUV-driven escape, as opposed to reproducing the remnant planet core exactly.
In Figs. 10, we display a heatmap of the starting conditions which resulted in planet with mass and radius consistent with  2 Lupi b within 3-, similar to 5 for TOI-431 d.There are two extra grey regions on this plot.The dark grey shows combinations of the starting conditions that implied a core mass more than 3- different from the measured mass today, and were thus excluded from our analysis.The light grey region shows where MESA failed to reach a stable solution in the earlier steps of setting the planet model up (which uses existing code that we did not edit in our study), and so the evolution under XUV irradiation was not able to be tested.For the tracks that were found to be consistent, we also plotted their evolution of the key planetary characteristics in Fig. 11.
The vast majority of the consistent tracks were found using the hydro-based approximation, while a few were found with ATES.Starting masses between 5 and 6 M⊕ were favoured, although almost all between 4 and 5 M⊕ failed to find a stable solution in setting up the model in MESA and so we not able to be tested.The successful starting envelope mass fractions were over a wider range, with viable simulations for all values tested below 0.15, although both successful ATES tracks started at 0.01, the lowest initial value.There is no obvious relationship between starting envelope mass fraction and starting mass in terms of which combinations give consistent results.Energy-limited escape was not able to reproduce the planet for any starting conditions we evolved under XUV-driven escape.However, with the starting masses between 4 and 5 M⊕ largely unable to be tested, we cannot completely rule out the possibility of the method being able to reproduce  2 Lupi b.This is again though suggestive of the method underestimating mass loss rates as compared to those based on the hydrodynamic models, such as the other two methods we employ here.
As an additional test of which starting conditions could lead to a stripped planet, we also investigated simply setting a threshold envelope mass fraction at the end of the simulations.To determine this threshold, we examined the last step of each of the consistent tracks within 3- from above, and found the highest envelope mass fraction among them: 0.0002.If we assume anything below this is stripped (or imminently will be), a few extra starting conditions result in a stripped planet, as shown in Fig. 12.The typical ranges of the two starting conditions is largely unchanged, with the extra "successful" tracks in this case most filling in gaps within the consistent region from Fig. 10.The energy-limited method still fails to remove the envelope of the planet models under this threshold definition of stripping.
Finally, we tested what happened if we adopted an XUV history of  2 Lupi that was at the median level for a star of its mass.The median Mors track suggests the star was under-bright in X-rays for its age in our observations.With four test cases that did not strip for any method in our main results (pairings of starting masses of 6, 6.5, 7.0, and 7.5 with starting envelope mass fractions of 0.05, 0.10, 0.15, and 0.21, respectively).The hydro-based approximation strips all four of these cases, while ATES strips all but the highest mass case.Energy-limited still fails to reproduce the stripped core composition.In general, these results suggest that if  2 Lupi happened to be at a lower than average brightness during our observations, then many more starting conditions may have been possible than suggested by Fig 10.

Current Mass Loss rates
In Table 5 we give current mass loss rates of the three transiting planets in our sample that have densities consistent with them being gaseous.We compared the different mass loss rate estimation methods used in Section 4, along with two extra methods: the interpolator for the Kubyshkina et al. (2018a) models (Kubyshkina & Fossati 2021), and running the full hydrodynamic Ates code (Caldiroli et al. 2021).For  2 Lupi d, the code failed to converge in a reasonable amount of time, and so we do not give a value for this method.The failure to converge may be linked to its current low irradiation level of 14.1 +3.0 −3.4 erg s −1 cm −2 , meaning the planet may be in the Jeans Starting Mass (M ⊕ )  escape regime, and thus modelling it with hydrodynamic escape at current age may be inappropriate.We did not calculate current mass loss rates for TOI-431 b and  2 Lupi b because our mass loss rate methods all assume a hydrogendominated atmosphere.These two planets are thought to not retain any envelope and thus the applicability of these methods at the current epoch is questionable.
Looking at Table 5, we can compare the results across the different methods, particularly those based on the same models.The mass loss rates given by the hydro-based approximation,  HBA , are all systematically smaller than those given by the interpolator tool,  Int , by factor of between two and five.From fig. 3 of Kubyshkina et al. (2018b), this factor appears to be in line with expectation.Overall, one would expect the interpolator to be a better approximation to the Kubyshkina et al. (2018a) models, as it provides a local interpolation to that region of the grid, as opposed to a single analytic expression across all models.Meanwhile, those calculated with the Ates code,  ATES , are encouragingly matched very well by the analytic fit to those models: the approximate effective efficiency function,  AEEF .Finally, all five methods agree that TOI-431 d is undergoing the most mass loss currently, and  2 Lupi d the least.

TESTING THE SYSTEMS WITH EVAPMASS
As both TOI-431 and  2 Lupi both have planets above and below the radius valley, we used the EvapMass code (Owen & Campos Estrada 2020) as an additional test of the systems' compatibility with the photoevaporation mechanism.This code takes the parameters of a rocky planet in the system and calculates the minimum mass of a second, gaseous planet to be consistent with the photoevaporation mechanism.The code exploits the fact that the two planets have experienced the same XUV emission from the star, scaled to their respective orbital distances.In their study, Owen & Campos Estrada (2020) explored 73 systems to determine if they are consistent with the photoevaporation model, finding only two systems that appear inconsistent.
For TOI-431, we performed the test with planets b and d as the stripped and gaseous planets, respectively, and found them to be consistent with photoevaporation.However, the constraint provided is not very informative.Planet b is so small and close to the star that it provides very little constraint on planet d.The minimum mass required for consistency is just 0.099 M ⊕ at 95 per cent confidence, compared to the true mass of planet d which is 9.9 ± 1.5 M ⊕ .
With  2 Lupi, we applied the code to two planet pairs: b and c, and b and d.Planet d's required minimum mass is a similarly low, 0.59 M ⊕ at 95 per cent confidence.For planet c, the constraint is somewhat more informative, with a minimum mass of 2.7 M ⊕ .However, the measured mass of this planet is again far above this minimum confidence limit, at 11.24 +0.65  −0.63 M ⊕ .This method demonstrates that the two systems in this work are consistent with the photoevaporation model, in agreement with the results of our simulations in Section 4, though its minimum mass constraints are not strong.The EvapMass code is more informative for systems with a gaseous planet interior to the stripped one.One such example is K2-3, for which the measured parameters are not consistent with photoevaporation, according to the minimum mass limits calculated by the code (Diamond-Lowe et al. 2022).

CONCLUSIONS
We have examined the TOI-431 and  2 Lupi planetary systems, both of which host planets either side of the radius valley, in the context of XUV photoevaporation.Using XMM-Newton EPIC observations of the systems, we measured the stars' X-ray fluxes.These measurements will prove useful in interpreting any future observations targeting the detection of escaping or extended atmospheres of the planets in these systems, such as at Ly-.The X-ray and NUV light curves for TOI-431 show evidence of a possible stellar flare with a rise in both at the end of the observation.There is, however, no evidence of the emission hardening at that epoch, meaning that we cannot conclusively determine if the increase in X-ray and UV flux can definitively be attributed to a stellar flare.
Using the current X-ray flux as an anchor for the XUV evolution of each host star, we simulated the lifetime evolution of a H/He envelope for each of the five transiting planets across the two systems.Building on the work of King et al. (2022), we simulated the past evolution as well as the future, starting at an age of 5 Myr.For the XUV time evolution, we generated evolution curves with the Mors code (Johnstone et al. 2021), and anchored them to the XUV fluxes we measured with the XMM-Newton observations.To get the planet structure, and therefore radius, at each step, we used MESA to evolve our planets, building on previous work by (e.g.Chen & Rogers 2016;Malsky & Rogers 2020).We used multiple methods of estimating the mass loss rate throughout the observations, comparing the predicted mass and radius to those observed today.
Our results demonstrate TOI-431 d can be easily reproduced by all mass loss methods we employed, though the viable starting conditions varied between these methods.Energy-limited favoured a smaller starting envelope mass fraction of 0.07, while the hydrobased approximation to the Kubyshkina et al. (2018a) models suggested it could be as high as 0.25, though most successful models were found just either side of 0.1.Our simulations for  2 Lupi c and d encountered issues with the equation of state tables used by MESA, causing numerous tracks to terminate prior to the older system age of 12.3 Gyr.However, none of the simulations were able to strip the envelope off either planet, or were on track to do so.We conclude that the survival of their envelopes to the present day, as shown by their measured densities, is consistent with the XUV-driven escape model.
Our simulations for TOI-431 b also encountered issues due to high bolometric flux at the planet in the final setup stage, leading to very large radii.However, a quick calculation suggests that removing a 20 per cent initial primordial envelope could be achieved on Myr timescales due to the very small orbital separation and thus intense XUV irradiation the planet would have experienced immediately following disk dissipation.
Our simulations were also able to completely remove the primordial envelopes of  2 Lupi b, for at least some combinations of starting mass and envelope mass fraction.Complete stripping is in line with expectations from its measured density.While none of the energylimited method simulations resulted in a complete removal, the other two methods were able to strip the planet in some cases, Furthermore, there were a number of lower starting mass simulations that failed in MESA before we reached the XUV evolution stage, and so were not able to be tested.
Our simulations, along with use of the EvapMass code (Owen & Campos Estrada 2020), show that XUV-driven escape can successfully reproduce the system architectures seen, with the innermost planet stripped, and the outer transiting planets able to retain at least some of their primordial envelope.These systems therefore appear to be consistent with the XUV mechanism for forming the radius valley.Whether it is the dominant mechanism, for these systems or indeed in general, remains an open question.

Figure 3 .Figure 4 .
Figure 3. X-ray spectrum for TOI-431, observed by the EPIC-pn camera.For details of the overlaid model fit, see the main text.

Figure 5 .
Figure 5. Heat map for the TOI-431 d simulations, displaying the starting conditions which result in evolutionary tracks where the radius and mass are consistent with the measured values to within 1- at the current estimated age.The squares are coloured according to the number of mass loss calculation methods that produce a consistent result.The labels display which method(s) produced a consistent result: E -energy limited; A -Ates approximation; Hanalytic hydro-based approximation to theKubyshkina et al. (2018a) models(Kubyshkina et al. 2018b).Some simulation tracks lie outside the limits of the plot, but none of those produced any consistent tracks.

Figure 6 .
Figure6.The subset of evolutionary tracks for TOI-431 d where the planet's radius (marked as a black dotted line in the top right panel) and mass matched the measured values at the current estimated age (marked as a solid black line in all panels) to within 1- of both.The top left panel shows the evolution of the X-ray and EUV luminosities in time by plotting their respective ratios with the bolometric luminosity .The other five panels show how key characteristics of the planet, and the escape of material from it evolve in time across the three mass loss rate methods we used.

Figure 8 .
Figure8.Heat map for the  2 Lupi c simulations, displaying those evolutionary tracks ended before an age of 11 Gyr.The squares are coloured according to the number of mass loss calculation methods that suffered an early end to the simulation at that starting point in the grid, while the text labels are as in Fig.5.We used lower case letters here, as well as a different colour bar, to distinguish from that plot however, as they display very different things.

Figure 10 .
Figure10.Heat map for the  2 Lupi b, displaying the starting conditions which result in evolutionary tracks where the radius and mass are consistent within the measured values to within 3- at the closest timestep to the current estimate age.Note that in some cases, the tracks end many Gyr before the estimated age.Light grey regions are those where MESA failed to reach the evolution step.Dark grey regions in the top left and bottom right are excluded as they imply a core mass that is not within 3- of the current measured mass, as we assume the core mass is unchanged during the evolution.The text labels indicating the consistent models are as in Fig.5.

Figure 11 .
Figure11.The subset of evolutionary tracks for  2 Lupi b where the planet's radius (marked as a black dotted line in the top right panel) and mass matched the measured values at the current estimated age (marked as a solid black line in all panels) to within 3- of both.The layout and labels are as in Fig.6.

Figure 12 .
Figure 12.As Fig. 10 but showing all tracks whose final envelope mass fraction was below 0.0002 (see main text), regardless of the age of the stopping point of the simulation.

Table 2 .
Delrez et al. (2021)parameters for the transiting planets in the systems we investigated.Parameters are taken fromOsborn et al. (2021)andDelrez et al. (2021)for the TOI-431 and  2 Lupi systems, respectively.The numbers in brackets for the  orb values are the uncertainties on the final two decimal places.

Table 3 .
Best fit temperatures, emission measures, and unabsorbed fluxes at Earth for the fits to our X-ray spectra of TOI-431 and  2 Lupi.We also give luminosities and fluxes at the planets calculated using from  x .

Table 4 .
Choices for the starting mass and envelope mass fraction grid for the planets in the simulations.Note that for TOI-431 b and  2 Lupi b, we only use combinations of the starting mass and envelope mass fraction grids listed that imply a core mass within 3- of the current measured mass of the planet. *
The Ates code did not successfully converge for  2 Lupi d (see main text).