TOI-2374 b and TOI-3071 b: two metal-rich sub-Saturns well within the Neptunian desert^††thanks: This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.

TESS observed TOI-2374, with TESS Input Catalog (TIC) ID 439366538, in sector 1 (UT 2018 July 25–August 22) with a 30 minute cadence and in sector 28 (UT 2020 July 30–August 26) with a 10 minute cadence, both on camera 1 CCD 4. TOI-2374 b was detected in the Full-frame images (FFIs) by the MIT’s Quick-Look Pipeline (QLP, Huang et al., 2020) and alerted as a TESS object of interest (TOI) on 28 October 2020²²2https://exofop.ipac.caltech.edu/tess/target.php?id=439366538 (Guerrero et al., 2021). The detection gave a period of $4.31362\pm 0.00001$ d, an epoch of $2458326.555387$  days in BJD, and a depth of $9.490\pm 0.008$ parts-per-thousand (ppt).

TOI-3071 (TIC ID 452006073) was monitored by TESS in sectors 10 (UT 2019 March 26–April 22) and 37 (UT 2021 April 2–28) with camera 2 CCD 2 with a 30 and 10 minute cadence FFI, respectively. The QLP detected the candidate and it was alerted as a TOI on 4 June 2021³³3https://exofop.ipac.caltech.edu/tess/target.php?id=452006073. The detection gave a period of $1.2671\pm 0.0003$ d, an epoch of $2459331.9673\pm 0.0017$  days in BJD, and a depth of $2.580\pm 0.002$ ppt. It was subsequently observed in sector 64 (UT 2023 April 6-May 4) on camera 2 CCD 1 as a 2-minute cadence target. The Science Processing Operations Center (SPOC) at NASA Ames Research Center conducted a transit with a noise-compensating matched filter (Jenkins, 2002; Jenkins et al., 2010, 2020) and detected the transit signature of TOI-3071 b. An initial limb-darkened model fit was performed (Li et al., 2019) and a suite of diagnostic tests were conducted help make or break the planetary nature of the signal (Twicken et al., 2018). The candidate transit signature passed all these diagnostic tests including the difference image centroiding test, which constrained the location of the host star to within 2.1 +/- 2.6 arc sec of the transit source.

For observations of TOI-2374 from TESS sectors 1 and 28 and observations of TOI-3071 from sectors 10 and 37, we downloaded the light curves computed by the QLP and used the Kepler Spline Simple Aperture Photometry (KSPSAP) flux, which is detrended by applying a high-pass filter to remove low-frequency variability from stellar activity or instrumental noise (Huang et al., 2020). We then used the lightkurve python package (Lightkurve Collaboration et al., 2018) to remove all data points flagged as being affected by excess noise (corresponding to stray light from Earth or the Moon in the camera field of view or scattered light). For TOI-3071 observations from sector 64, we downloaded the photometry provided by the SPOC pipeline, and used the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP), from which common trends and artefacts, as well as an estimated contamination from blended sources, have been removed by the SPOC Presearch Data Conditioning (PDC) algorithm (Twicken et al., 2010; Smith et al., 2012; Stumpe et al., 2012; Stumpe et al., 2014). No further detrending of the light curves was deemed necessary as they are relatively flat across the whole time series. The median-normalised light curves for both targets were used in the joint model outlined in section 3.3.

TIC contamination ratios of 6.65 and 0.62 are reported for TOI-2374 and TOI-3071 respectively. These values are computed as the nominal flux from the contaminants divided by the flux from the source. The large values reported for both targets are most likely due to close sources as can be seen in the TESS Target Pixel Files images shown in Figure 1 (created with tpfplotter⁴⁴4Publicly available at https://github.com/jlillo/tpfplotter, Aller et al. 2020). The TPFs show several sources of potential contamination around the target stars (particularly TOI-2374, which has a very close source two magnitudes brighter). Lightcurves from QLP and SPOC data processing pipelines are approximately corrected for blending from nearby stars. The SPOC pipeline computes a crowding metric using a series of assumptions described in Section 2.3.11 of Stumpe et al. (2012). The QLP uses a different set of assumptions described in step 7 in Huang et al. (2020). To further measure the extent of this blending and contribute to the correction, follow up ground-based photometry was necessary; the details of these observations are described below. Additionally, the TPF for TOI-2374 shows that both the target and the bright companion have the same proper motion direction so they might be bounded companions.

The TESS pixel scale is $\sim 21\arcsec$ pixel^-1 and photometric apertures typically extend out to roughly 1 arcminute, generally causing multiple stars to blend in the TESS photometric aperture. This is precisely the case for TOI-2374 and TOI-3071 (Fig.1). To determine the true source of the TESS detection, we acquired ground-based time-series follow-up photometry of the field around both TOI-2374 and TOI-3071 as part of the TESS Follow-up Observing Program (TFOP; Collins et al., 2018)⁵⁵5https://tess.mit.edu/followup. The follow-up light curves are also used to confirm the transit depth and thus the TESS photometric deblending factor and refine the TESS ephemeris. We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen, 2013), to schedule our transit observations. The observations are summarised below and confirm the transit-like events detected by TESS to be occurring within the small TOI-2374 and TOI-3071 follow-up photometric apertures.

We obtained radial velocity (RV) measurements of TOI-2374 and TOI-3071 with the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph mounted on the ESO 3.6 m telescope at the La Silla Observatory in Chile (Mayor et al., 2003a). Under the HARPS-NOMADS large programme (ID 1108.C-0697, PI: Armstrong), a total of 21 spectra were obtained between 9 April 2022 and 7 July 2022 for TOI-2374 and 14 spectra between 19 March and 2 April 2022 for TOI-3071. The instrument (with resolving power $R=115\,000$) was used in high-accuracy mode (HAM). The exposure time for the TOI-2374 observations was 1800 s. This resulted in a signal to noise ratio (SNR) of between 22 and 35 per pixel. The exposure time for the TOI-3071 observations was 2400 s, resulting in a SNR of between 10 and 25 per pixel. The data were reduced using the standard offline HARPS data reduction pipeline, and a K5 template (TOI-2374) or G2 template (TOI-3071) was used to form a weighted cross-correlation function (CCF) to determine the RV values (Baranne et al., 1996; Pepe et al., 2002). The line bisector (BIS) and full-width at half-maximum (FWHM) were measured using methods by Boisse et al. (2011).

Three bad data points were removed from the TOI-2374 dataset. One of them (BJD = 2459711.8) was reduced incorrectly and the other two (BJD = 2459712.8 and BJD = 2459716.8) contained large errors due to bad weather during the observations. The achieved instrumental precision for the rest of the RV measurements is $\approx 4$ m s^-1. The observations are presented in tables 4 and 5.

We observed TOI-2374 with the Planet Finder Spectrograph (PFS; Crane et al. 2006; Crane et al. 2008, 2010) on the Magellan/Clay telescope at Las Campanas Observatory in Chile. We obtained seven observations of this target between 2021 Sep 14 and 2021 Nov 17 through an iodine cell, with relatively short exposure times of 5 minutes. These observations were made in 3x3 binning mode to reduce readout noise. We also obtained an iodine-free template spectrum at higher signal to noise, which we used to extract relative RVs. The RVs were extracted using the pipeline described in Butler et al. (1996), achieving a typical instrumental RV precision of $\approx 3$ m s^-1. We also quantified potential variations in the spectral line profiles using bisector inverse slopes (BIS) measurements. These were derived using the procedures described in Hartman et al. (2019), applied to the iodine-free orders of each spectrum. The BIS measurements are shifted such that the median bisector span for each order is zero. All measurements are presented in Table 6.

We obtained two reconnaissance spectra on UT 2021 July 21 and UT 2021 August 5 of TOI-2374 using the Tillinghast Reflector Echelle Spectrograph (TRES, Fűrész, 2008) located at the Fred Lawrence Whipple Observatory (FLWO) in Arizona, USA. TRES has a resolving power of $\sim$44,000 and operates in the wavelength range 390-910 nm. The spectra were extracted as described in Buchhave et al. (2010).

The stellar spectroscopic parameters ($T_{\mathrm{eff}}$, $\log g$, microturbulence, [Fe/H]) were estimated using the ARES+MOOG methodology from the respective combined HARPS spectrum of each star. The methodology is described in detail in Sousa et al. (2021); Sousa (2014); Santos et al. (2013). We used the latest version of ARES ⁷⁷7The last version, ARES v2, can be downloaded at https://github.com/sousasag/ARES (Sousa et al., 2007, 2015) to consistently measure the equivalent widths (EW) on the list of iron lines. For TOI-3071 we used the list of lines presented in Sousa et al. (2008), while for TOI-2374 we used the list of lines presented in Tsantaki et al. (2013) which is more appropriate for stars cooler than 5200 K. In the analysis we use a minimization process to find the ionization and excitation equilibrium to converge for the best set of spectroscopic parameters. This process makes use of a grid of Kurucz model atmospheres (Kurucz, 1993) and the radiative transfer code MOOG (Sneden, 1973). We also derived a more accurate trigonometric surface gravity using recent GAIA data following the same procedure as described in Sousa et al. (2021). We estimated rotational velocity $v\sin i$ values by performing spectral synthesis (with the same code and model atmospheres as for stellar parameters) of two iron lines in the 6705Å region. We fix the macroturbulence velocity with the empirical formula by Doyle et al. (2014), which provides v_mac=4.6 km s^-1 for TOI-3071. Cool stars are outside of such calibration, hence we adopted a value v_mac=2 km s^-1 for TOI-2374. Both stars seem to be slow rotators with $v\sin i$ = 2.0 km s^-1 for TOI-3071 and $v\sin i$ $<$ 0.5 km s^-1 for TOI-2374. The stellar parameters are presented in tables 1 and 2.

Independently, we derived the stellar atmospheric parameters for TOI-2374 from the TRES spectra using the Stellar Parameter Classification (SPC, Buchhave et al., 2012) tool. SPC cross correlates an observed spectrum against a grid of synthetic spectra based on Kurucz atmospheric models (Kurucz, 1992). With this method we obtained $T_{\mathrm{eff}}=4937\pm 50$ K, $\log g=4.61\pm 0.10$ and $[\mathrm{Fe/H}]=0.24\pm 0.08$, which are in agreement (within errors) with the results derived from the HARPS spectra.

Stellar abundances of the elements were derived using the same tools and models as for stellar parameter determination under the assumption of local thermodynamic equilibrium. For the derivation of chemical abundances of refractory elements we closely followed the methods described in e.g. Adibekyan et al. (2012); Adibekyan et al. (2015); Delgado Mena et al. (2014); Delgado Mena et al. (2017). Abundances of the volatile elements, C and O, were derived following the method of Delgado Mena et al. (2021) and Bertran de Lis et al. (2015). Nine out of the 14 individual spectra of TOI-3071 were contaminated around the 6300.3 Å oxygen line and had to be discarded before being coadded, leading to a lower S/N spectrum and thus a higher error in oxygen abundance. These two elements could not be derived for TOI-2374 because their lines are very weak in cool stars spectra. All the [X/H] ratios are obtained by doing a differential analysis with respect to a high S/N solar (Vesta) spectrum from HARPS. The abundances of the elements are shown in table 7 in the appendix.

From the parameters obtained in the spectroscopic analyisis we estimated the stellar radius and mass using the calibrations presented in Torres et al. (2010) (see tables 1 and 2). In addition, we performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia DR3 parallax, in order to determine an empirical measurement of the stellar radius (Stassun & Torres, 2016; Stassun et al., 2017; Stassun et al., 2018). We pulled the $JHK_{S}$ magnitudes from 2MASS, the W1–W3 magnitudes from WISE, and the $G_{\rm BP}G_{\rm RP}$ magnitudes from Gaia. We also utilized the absolute flux-calibrated Gaia spectrum where available. Together, the available photometry spans the full stellar SED over the wavelength range 0.4–10 $\mu$m (see figures 4 and 5). We performed a fit using PHOENIX stellar atmosphere models (Husser, T.-O. et al., 2013), adopting from the spectroscopic analysis the effective temperature ($T_{\rm eff}$), metallicity ([Fe/H]), and surface gravity ($\log g$). We fitted for the extinction $A_{V}$, limited to the maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998). The resulting fits (figures 4 and 5) have $A_{V}=0.02\pm 0.02$ and $0.20\pm 0.03$ for TOI-2374 and TOI-3071, respectively, with a reduced $\chi^{2}$ of 1.3 and 1.4, respectively. Integrating the (unreddened) model SED gives the bolometric flux at Earth, $F_{\rm bol}=4.981\pm 0.058\times 10^{-10}$ erg s^-1 cm^-2 and $3.291\pm 0.038\times 10^{-10}$ erg s^-1 cm^-2, respectively. Taking the $F_{\rm bol}$ together with the Gaia parallax directly gives the bolometric luminosity, $L_{\rm bol}=0.2867\pm 0.0034$ L_⊙ and $2.539\pm 0.032$ L_⊙, respectively. The Stefan-Boltzmann relation then gives the stellar radius, $R_{\star}=0.774\pm 0.025$ R_⊙ and $1.393\pm 0.030$ R_⊙, respectively. These estimates differ by $2.6\sigma$ and $1.7\sigma$ with the ones obtained from the spectroscopic analysis, respectively. Given that the two sets disagree by less than $3\sigma$ and that the former analysis works with the full spectrum, we decided to use the ones estimated from the calibrations in Torres et al. (2010) as input for the joint model described in section 3.3.

We used the relations stated in Noyes et al. (1984) to derive chromospheric activity index $\log R^{\prime}_{\rm HK}$ values from the HARPS measurements of the Mount-Wilson S-index. With these values we obtained estimates of the rotational period of the stars $P_{\rm rot}$. For both stars, $P_{\rm rot}$ is consistent with the maximum rotation period calculated with the $v\sin i$ estimates, $P_{\rm rot}/\sin i>70$ d for TOI-2374 and $33\pm 8$ d for TOI-3071. We used Eq. 3 and Eq. 16 in Barnes (2007) to calculate the gyrochronological stellar ages from the rotation periods, giving 4.34 $\pm$ 0.67 Gyr and 4.9 $\pm$ 0.8 Gyr for TOI-2374 and TOI-3071, respectively.

Moreover, we obtained two other independent estimates of the stellar ages $\tau$: first, we used the web interface PARAM 1.3 ⁸⁸8http://stev.oapd.inaf.it/cgi-bin/param_1.3, which estimates the basic intrinsic parameters of stars given their photometric and spectroscopic data following the method described in da Silva et al. (2006). This resulted in stellar age estimates of 4.8 $\pm$ 4.2 Gyr and 0.8 $\pm$ 0.6 Gyr for TOI-2374 and TOI-3071, respectively. As expected, the age determination based on evolutionary models is extremely uncertain for main-sequence stars. The age of TOI-2374 formally agrees with the one from gyrochronology. On the other hand, PARAM provides a much younger age for TOI-3071. We then used the chemical abundances of some elements to derive ages through the so-called chemical clocks (i.e. certain chemical abundance ratios which have a strong correlation for age). We applied the 3D formulas described in Table 10 of Delgado Mena et al. (2019), which also consider the variation in age produced by the effective temperature and iron abundance. The chemical clocks [Y/Mg], [Y/Zn], [Y/Ti], [Y/Si], [Y/Al], [Sr/Ti], [Sr/Mg] and [Sr/Si] were used from which we obtain a weighted average age of 1.1 $\pm$ 0.2 Gyr for TOI-3071. This age is again much younger than the age from gyrochornology. The ages from chemical clocks for cool stars such as TOI-2374 must be taken with caution since the calibrations are obtained using hotter stars. Therefore, we choose to apply the 2D formulas considering the abundance ratios and the metallicity (see Table 8 in Delgado Mena et al., 2019) which provide a weighted average age estimate of 3.3 $\pm$ 0.7 Gyr. Ages from each individual clock are shown in Table 8. We include in Table 1 the geometric average of the three stellar age estimates for TOI-2374. The error is given by the deviation of the three values. In Table 2 we report two values for the age of TOI-3071, corresponding to the two disagreen results: the gyrochronological age and the weighted average of the ages derived with PARAM and with the chemical clocks.

To model the evaporation histories of the two planets, we assumed a simpler two-layer internal structure with a rocky core surrounded by a pure H/He envelope. This model can be described with four quantities: the rocky core mass M${}_{\text{core}}$ and core radius R${}_{\text{core}}$, and the gaseous envelope thickness R${}_{\text{env}}$ and envelope mass fraction $\text{f}_{\text{env}}$, defined as the ratio of envelope mass to planet mass $\text{M}_{\text{env}}/\text{M}_{\text{p}}$. In order to determine these parameters, we adopted the empirical mass-radius relation by Otegi et al. (2020) for the rocky core, and the envelope structure model by Chen & Rogers (2016) for the gaseous atmosphere, which is based on MESA simulations. Uisng this model, we thus determined envelope mass fractions of $38\pm 5$ % for TOI-2374 b and $26\pm 6$ % for TOI-3071 b. These correspond to core (heavy-element) masses of M${}_{\text{core}}=34.7\pm 3.7$ for TOI-2374 b, and M${}_{\text{core}}=50.5\pm 4.9$ for TOI-3071 b, both consistent with the heavy-element masses determined in Sect. 4.2 within $1\sigma$.

Simulating the evolution of a planet’s atmosphere under photoevaporation requires knowledge of its X-ray irradiation history. We can estimate the X-ray and EUV emission (together, XUV) of a star from its spin period, as the two are related via rotation–activity relations (e.g. Wright et al., 2011; Pizzolato et al., 2003). We adopted the stellar rotation evolution models by Johnstone et al. (2021), which simulate the evolution of the angular momentum of stars constrained by observations of young star clusters. Moreover, to estimate the EUV component, we adopted the scaling relation by King & Wheatley (2021), which are based on Solar TIMED/SEE data. We plot the evolution of the spin period and corresponding XUV emission of TOI-2374 and TOI-3071 on Figures 16 and 17, respectively.

Regarding TOI-2374, we found the rotational models by Johnstone et al. (2021) are consistent with the measured age and spin period of the star.

For TOI-3071, we considered the two age estimates separately: an age of 1 Gyr from the chemical clocks and PARAM isochrones, and an age of 5 Gyr from gyrochronology (see Sect. 3.1.3).

We found the choice of age had a strong impact on the X-ray emission history of the star. The younger age is more consistent with a short spin period of 2-4 d, whereas the older age from gyrochronology is more consistent with a slower spin period of 20-30 d, as the star has had more time to spin down. The spin period of the star, however, is constrained to 26 d, and so adopting the younger age would imply that the star must be unusually slowly rotating and X-ray faint (Fig. 17, red line).

For that reason, we adopted two scenarios for the X-ray emission history of TOI-3071: (1) the high activity scenario where we adopted the predicted evolution from the models by Johnstone et al. (2021), consistent with the older gyrochronological age (see Fig. 17, black line); and (2) the low activity scenario, where we adopted a fainter emission history model consistent with the younger age from the 3D chemical clocks (Fig. 17, red line).

We then simulated the past and future evaporation histories of the planets with the photoevolver code (Fernández Fernández et al., 2023). The simulation is built upon three ingredients: (1) the XUV emission history of the star, which determines the X-ray irradiation on the planet, (2) a mass loss formulation, which determines the amount of gas lost due to XUV irradiation, and (3) an envelope structure model, which relates the mass of the gaseous atmosphere to its thickness. Therefore, on each time step, the X-ray luminosity of the star at that age is used to calculate the mass loss rate, which is then used to recalculate the mass and size of the planet.

We adopted the models by Johnstone et al. (2021) for the stellar XUV emission history (described above), and the model by Chen & Rogers (2016) for the structure of the gaseous atmosphere, described in Section 4.3.1. Regarding the mass loss formulation, we adopted the model by Kubyshkina et al. (2018). We evolved the planets’ atmospheres using the 4th order Runge-Kutta algorithm (Dormand & Prince, 1980) as the integration method with a variable time step no larger than 1 Myr. We ran the simulation backwards to the age of 10 Myr and forwards to 10 Gyr. In the case of the F-type TOI-3071, we stop the simulation at 6 Gyr, the approximate age at which the star leaves the main sequence (Choi et al., 2016).

Moreover, for each of the two planets, we also explored a range of possible evaporation histories taking into account the uncertainties on the planet’s mass, radius, and XUV irradiation. To maximise the evaporation rate, we picked a radius on the higher end of its $1\sigma$ uncertainty and a mass on the lower end of its $1\sigma$ uncertainty, which minimises the density and maximises the escape rates. We also picked the highest XUV emission history allowed by the model’s $2\sigma$ spread, as shown on Figures 16 and 17. Likewise, to minimise the evaporation rate, we picked the opposite parameters: values for the mass and radius within their $1\sigma$ errors that maximise density under the faintest XUV irradiation history allowed by the model’s $2\sigma$ spread.

The evaporation histories of TOI-2374 b and TOI-3071 b are plotted on Figure 18 (left and right panels, respectively).

We found that TOI-2374 b is stable against evaporation, as its envelope mass fraction changes only by a few percent across its lifetime. The evolution of its radius is thus largely driven by thermal contraction. We also constrained its current mass loss rates to $0.4-0.8\times 10^{8}$ g s^-1.

On the other hand, TOI-3071 b is more susceptible to escape despite its higher mass. This is due to being twice as close to its star compared to TOI-2374 b, as well as the higher X-ray luminosity of its F-type host star. Under the higher X-ray irradiation history consistent with the gyrochronological age (5 Gyr), our simulations show that the planet is reaching the point of losing its atmosphere completely when its star starts evolving off the main sequence (see Fig. 18), with an envelope mass fraction of 10-20% at 6 Gyr down from 40-60% initially. Under this X-ray irradiation history, the planet currently experiences mass loss rates of $7.5-40.7\times 10^{8}$ g s^-1.

However, under the fainter irradiation history consistent with the much younger age derived from chemical clocks (1 Gyr), TOI-3071 b experiences lower escape rates throughout its life and is stable against evaporation. We constrained its current mass loss rates to $0.9-2.6\times 10^{8}$ g s^-1, one order of magnitude lower compared to the high activity scenario.

In Sect. 4.2, we presented more detailed internal structure models to estimate the bulk metallicity of the two planets. These models cannot predict the atmospheric composition, which can only be determined by observations. However, since both are metal-rich, they are also likely to have high-metallicity atmospheres. Such high-metallicity atmospheres experience much lower escape rates as heavier species require more energy to remove (Wilson et al., 2022; Owen & Jackson, 2012). Since we only considered a pure H/He composition in our atmospheric evolution models, our results in Fig. 18 constitute an upper bound to the evaporation histories of the planets.

For TOI-2374,b, the atmosphere is stable against evaporation throughout its life, even when assuming a pure hydrogen envelope. TOI-3071,b experiences moderate escape rates and is significantly affected by evaporation only as its host star transitions out of the main sequence, risking total atmospheric loss. However, heavy-mass planets like TOI-3071 b could present a metallicity gradient in their gaseous envelopes, with a relatively pure H/He upper atmosphere that progresses towards metal-rich layers deeper in the atmosphere. If the planet experienced stronger escape rates earlier in its life, either from photoevaporation or Roche-lobe overflow, the pure H/He upper layers of the envelope would be stripped off first, exposing the deeper metal-rich layers which would then curtail atmospheric escape (Owen & Murray-Clay, 2018; Hoyer et al., 2023).

Therefore, the presence of heavy species in the hydrogen-rich envelope of TOI-3071 b at the present time would lead to diminished escape rates, making its atmosphere more resistant to evaporation and thus likely stable. As a result, we find the two planets are not significantly affected by photoevaporation and are thus able to hold onto their gaseous envelopes throughout their lifetimes.