Detection of Atmospheric Escape from Four Young Mini-Neptunes

Planets between 1 and 4 R_⊕ are frequently found on close-in orbits around nearby Sun-like stars (Zhu & Dong 2021). This class of planets has no comparable analog in the solar system, and there is considerable debate about their origins (Drazkowska et al. 2022). A now-famous "radius gap" or "Fulton gap" (Fulton et al. 2017; Fulton & Petigura 2018) separates the smaller, denser super-Earths (1–1.7 R_⊕) from the larger, less dense mini-Neptunes (2–3 R_⊕). The conventional explanation for this gap is that mini-Neptunes have a primordial hydrogen/helium atmosphere comprising ∼1% of their total mass that significantly inflates their radii, whereas super-Earths have either lost their hydrogen-rich atmospheres or never acquired them in the first place (e.g., Lee & Connors 2021; Lee et al. 2022). If most super-Earths initially formed with hydrogen-rich envelopes, they could have been stripped away by intense X-ray and extreme UV (XUV) irradiation from the young star (photoevaporative mass loss; e.g., Owen & Wu 2017; Mills & Mazeh 2017). Alternatively, the energy powering the outflow may predominately come from the internal heat of formation ("core-powered" mass loss; Ginzburg et al. 2018; Gupta & Schlichting 2019). A third possibility is that some or most mini-Neptunes have no hydrogen-rich envelopes at all but instead formed with substantial water-rich envelopes (e.g., Mousis et al. 2020). In order to differentiate between these competing hypotheses, it is crucial to obtain observations of young mini-Neptunes to ascertain whether and how quickly they are losing hydrogen/helium.

Even though the study of escaping mini-Neptune atmospheres is in its infancy, observations of escaping atmospheres around larger planets have a long history, as do attempts to understand the observations with models. The first detection of Lyα absorption occurred 2 yr after the detection of the first transiting exoplanet (Vidal-Madjar et al. 2003). Because the Lyα core is wiped out by the interstellar medium and only the high-velocity wings can be seen, Lyα observations are only logarithmically sensitive to the mass-loss rate (Owen et al. 2023). Helium absorption observations hold more promise because interstellar absorption is not a problem, and the line probes the inner regions of the outflow. The 1083 nm helium triplet is a sensitive probe of escaping atmospheres (Oklopčić & Hirata 2018; Spake et al. 2018) and can be readily accessed using ground-based telescopes. This technique works best for planets with K dwarf hosts, as these stars are predicted to produce the largest metastable helium population (Oklopčić 2019). Escaping helium has been detected around numerous giant planets, including WASP-107b (Spake et al. 2018; Allart et al. 2019), WASP-69b (Nortmann et al. 2018), and HD 189733b (Salz et al. 2018; Guilluy et al. 2020; Zhang et al. 2022a). Models capable of predicting helium absorption include the Parker wind models of Oklopčić & Hirata (2018) and Lampón et al. (2020), multiple 1D hydrodynamic models (e.g., The PLUTO-CLOUDY Interface, TPCI; Salz et al. 2015a), the 3D Monte Carlo simulation EVaporating Exoplanets (Bourrier & Lecavelier des Etangs 2013; Allart et al. 2018, 2019), and many grid-based 3D hydrodynamic codes (Wang & Dai 2018; Khodachenko et al. 2021; Shaikhislamov et al. 2021; MacLeod & Oklopčić 2022; Rumenskikh et al. 2022). Although these models have been used to reproduce the observed outflows from giant exoplanets, they are largely untested in the mini-Neptune regime.

Earlier this year, we reported the first detections of ongoing mass loss from young mini-Neptunes. In Zhang et al. (2022b), we detected Lyα absorption from the 400 Myr mini-Neptune HD 63433c, which orbits a Sun-like star, but not from the inner mini-Neptune in the same system. This shows that the outer planet has an escaping hydrogen/helium atmosphere, but further interpretation is difficult because the interstellar medium blocks the core of the line, and the nondetection from the inner planet could be due to the absence of a lightweight atmosphere, a magnetic field preventing the outflow from being accelerated to high velocities, or photoionization of escaping hydrogen (Owen et al. 2023). In Zhang et al. (2022c), we reported metastable helium absorption from the outflowing atmosphere of TOI 560b, a 500 Myr mini-Neptune orbiting a K dwarf.

In these studies, our Microthena and TPCI models both struggled to reproduce the observations. The Microthena models could roughly match the depth of the Lyα blue wing absorption for HD 63433c but could not explain either the red wing absorption or the nondetection around HD 63433b. For TOI 560b, both the fiducial Microthena and TPCI models overestimated the observed helium absorption, and the Microthena models predicted a blueshift instead of the observed redshift. We suspect that the poor fit could be due to large uncertainties in planetary and stellar parameters, such as the EUV flux, the atmospheric metallicity, and the stellar wind conditions. It could also be due to limitations of the model, such as the neglect of magnetic fields, which could significantly affect the outflow (e.g., Owen & Adams 2014; Khodachenko et al. 2015; Arakcheev et al. 2017).

To gain a solid understanding of if and how mini-Neptunes turn into super-Earths, a large sample of mini-Neptunes with outflow detections is required, spanning a range of radii, masses, ages, and periods. It is much easier to accumulate a large sample with helium observations than with Lyα because the latter requires Hubble Space Telescope observations of very nearby stars that have favorable radial velocities, and we have not found any additional candidates as favorable as HD 63433. Having a larger sample of helium detections will reveal whether the features we detected in TOI 560b—such as the existence of an escaping hydrogen/helium atmosphere, the redshift of the signal, the asymmetry in the light curve, or the unexpectedly low amplitude of the absorption—are common or rare. Population studies can also reveal patterns in the data that give insight into the physical processes involved. By testing our models on a population level, we can mitigate some of the large uncertainties in individual planet and stellar parameters. By focusing on young planets that nevertheless span a factor of a few in age, we catch young mini-Neptunes during the most critical period of their lives: the period when they could be rapidly losing mass and transitioning into super-Earths.

For these reasons, we undertook a survey of escaping helium observations from eight young (<1 Gyr) mini-Neptunes orbiting K-type stars observable from the Keck telescopes. This paper presents the first four results from the survey. All of our targets are recent discoveries by the Transiting Exoplanet Survey Satellite (TESS). TESS has surveyed the entire sky, finding many more nearby transiting planets than previously known and making a survey of this type possible for the first time. Section 2 presents the observations and data reduction, Section 3 presents the results, Section 4 estimates the mass-loss rate using two different methods and comments on the implications, and Section 5 concludes.

2.1. Targets

TESS has already surveyed nearly the entire sky, discovering many nearby transiting planets amenable to atmospheric characterization. The TESS Object of Interest (TOI) catalog (Guerrero et al. 2021) contains a mix of new planet candidates and confirmed planets, including previously known planets redetected by TESS. To select targets for our Keck survey, we downloaded the catalog from ExoFOP (DOI: 10.26134/ExoFOP5) and imposed several cuts. We required planet candidates to have radii between 2 and 3 R_⊕, periods less than 20 days, declinations greater than −20°, and K dwarf host stars (3900 K < T_eff < 5300 K) with J < 9. For unpublished candidates, we required a rating of "Validated Planet Candidate?" (VPC?) or higher by the TESS Follow-up Observing Program (TFOP) in order to eliminate likely false positives from our sample. This designation means that follow-up observations have at least tentatively confirmed that the transit occurs for the correct star, and that there are no Gaia DR2 stars contaminating the TESS aperture that are bright enough to cause the detection.

Most of the planets and planet candidates we identified did not have well-determined ages, so we ran a Lomb–Scargle periodogram on the TESS light curve and required a robustly detected rotation period in order to limit our sample to young stars. In practice, this means that the star's rotation period must be less than 15 days. We manually examined the TESS light curve to make sure that the rotational modulations looked strong and convincing (coherent and on the order of tens of percent). We also searched for the stars in the ROSAT All-Sky Survey, an X-ray catalog. Detectable X-ray emission is also evidence of youth, and we therefore prioritized targets with detections in this catalog for our survey. In total, we are currently targeting eight planets orbiting seven stars, of which three have ROSAT detections. Although we continue to monitor the list of newly announced TESS candidates, we do not expect the survey sample to expand significantly because TESS has already completed its baseline survey and likely discovered all of the most favorable targets. There are additional favorable candidates remaining in the southern hemisphere, but these are not accessible from Keck.

During the first year of the survey, we observed four targets: TOI 560b, 1430.01, 1683.01, and 2076b. Their properties are summarized in Table 1. Objects TOI 560b and 2076b are both confirmed planets with published discovery papers (Barragán et al. 2022; El Mufti et al. 2021; Hedges et al. 2021, respectively). They are both the inner planets of multiplanetary systems, with one and two outer mini-Neptune companions, respectively. Objects TOI 1430 and 1683 are unpublished candidates, with the former included in the ROSAT All-Sky Survey, indicating its likely youth. Neither have any detected transiting companions. Of the four targets, only TOI 560b currently has a published radial velocity mass measurement. The masses of the other planets are derived from the mass–radius relation of Wolfgang et al. (2016). It is possible that this relation overestimates masses for young planets. Object TOI 560b gives us some reason for optimism, because the radial velocity–derived mass of ${10.2}_{-3.4}^{+3.4}$ M_⊕ is perfectly consistent with the mass inferred from the mass–radius relation, 10 ± 2 M_⊕. Similarly, the relation gives 8.4 ± 2 for TOI 560c, compared to a radial velocity mass of ${9.7}_{-1.7}^{+1.8}$ M_E. We encourage radial velocity follow-up of the other planets to aid the interpretation of atmospheric observations.

Table 1. Planet and Star Properties

Parameter/Planet	TOI 560b a	TOI 1430.01	TOI 1683.01	TOI 2076b b
Rp (R⊕)	2.79 ± 0.1	2.1 ± 0.2 d	2.3 ± 0.3 d	2.52 ± 0.04
Mp (M⊕)	${10.2}_{-3.1}^{+3.4}$	7 ± 2? f	8 ± 2? f	9? f
Teff,*	4511 ± 110	5067 ± 60 c	4539 ± 100 c	5200 ± 70
R* (R⊙)	0.65 ± 0.02	${0.784}_{-0.014}^{+0.018}$	${0.636}_{-0.017}^{+0.031}$	0.762 ± 0.016
M* (M⊙)	0.73 ± 0.02	0.85 ± j 0.10	0.69 ± 0.09 j	${0.824}_{-0.037}^{+0.035}$
Prot,*	12.08 ± 0.11	5.79 ± 0.15	11.3 ± 1.5	6.84 ± 0.58
J*	7.6	7.6	8.8	7.6
Age (Myr)	480–750	165 ± 30 g	500 ± 150 g	204 ± 50
a (au)	0.0604	0.0705	0.036	0.0631
FXUV (erg s−1 cm−2)	5000	6800	12,000	9500
FMUV (erg s−1 cm−2)	3000	27,000	4900	50,000
Teq (K)	714 ± 21	813	927	797 ± 12
P (days)	${6.3980420}_{-0.0000062}^{+0.0000067}$	7.434162 d ± 1.5 ×10-5	3.05752578 d ± 6.3 ×10-6	10.355183 ± 0.000065 i
T0 (BJDTDB − 2,457,000)	${1517.69013}_{-0.00059}^{+0.00056}$	2032.74865 ± 0.0010	2232.2348 ± 0.0011	2654.9837 ± 0.0015 i
a/R*	${19.98}_{-0.63}^{+0.61}$	17.25 ± 0.54	9.99 ± 0.39	17.8 ± 0.4
i (deg)	88.37 ± 0.18	${88.69}_{-0.23}^{+0.30}$	${85.99}_{-0.28}^{+0.29}$	88.9 ± 0.11 h
Transit duration τ14 (hr)	${2.143}_{-0.027}^{+0.029}$	${2.71}_{-0.04}^{+0.07}$	${1.28}_{-0.04}^{+0.06}$	3.251 ± 0.03

Notes.

^a Barragán et al. (2022). ^b Osborn et al. (2022). ^c Gaia Collaboration et al. (2018). ^d Own fits to TESS data. ^f Mass–radius relation of Wolfgang et al. (2016). ^g Gyrochronology using TESS-derived rotation period. ^h Hedges et al. (2021). ⁱ TOI 2076 has TTVs. See Appendix C. ^j TESS Input Catalog v8.2.

Download table as:  ASCIITypeset image

For TOI 1430.01 and 1683.01, which do not have published discovery papers, we took additional steps to validate their planetary nature. We checked that the even transits had a consistent transit depth with the odd transits. Object TOI 1683.01 was verified (in part) by the Dragonfly Telephoto Array, a robotic telescope stationed in New Mexico that is comprised of 48 individual telephoto lenses acting in concert similar to a 1.0 m refractor. With Dragonfly, we detected an on-time and on-target 1.6 part per thousand transit event, excluding the nearby Gaia stars as potential causes of a false-positive detection by TESS. Photometric data collected by the Las Cumbres Observatory have only tentatively detected the transit of TOI 1430.01 due to its shallow transit depth. We searched the Gaia DR2 catalog to ensure that no bright stars are within 21'' (1 TESS pixel) of the target, as this could contaminate the light curve. For TOI 1430.01, Gaia detected six stars within this radius, but their total flux in the Gaia G band is only 0.03% that of TOI 1430. Near TOI 1683, Gaia detected one star with a flux 0.025% of that of TOI 1683. Gaia has an angular resolution of 04, corresponding to 17 au at the distance of TOI 1430 and 21 au at the distance of TOI 1683. We therefore conclude that, when combined with our detection of helium absorption, these tests suggest that these two candidates are likely to be real planets and treat them as such in the subsequent analysis.

For these two planets, we obtained transit parameters by fitting the TESS light curves using the same methodology as Dai et al. (2021). To summarize, we downloaded the TESS photometry from the Mikulski Archive for Space Telescopes (MAST). We fit each light curve with batman (Kreidberg 2015), with free parameters including the out-of-transit flux level, stellar density, quadratic limb-darkening coefficients, period, transit midpoint, R_p /R_s , and orbital inclination. For both planets, we have multiple sectors of TESS observations spanning 2 yr, giving us a very accurate ephemeris that is sufficient to predict transit times corresponding to our helium observations to within 2–3 minutes. The phased light curves and transit fits are plotted in Appendix A.

We obtained age estimates for TOI 1430 and 1683 using gyrochronology. Specifically, we used Schlaufman (2010), who related age to stellar mass and rotation period, and Mamajek & Hillenbrand (2008), who related age to B – V color and rotation period. We obtained the rotation period from a Lomb–Scargle periodogram of the TESS simple aperture photometry (SAP) fluxes and took the masses and B – V colors from the TESS Input Catalog (Stassun et al. 2019). These stellar masses are computed from the effective temperature, which in turn is obtained from spectroscopy (if available) or dereddened colors. For TOI 1430, we obtain 156 ± 30 Myr with the first method and 174 ± 8 Myr using the second; for TOI 1683, we obtain 550 ± 150 and 460 ± 110 Myr. We adopt the average of the two methods and the more conservative error bars.

2.2. Keck/NIRSPEC

We used Keck/NIRSPEC to observe one transit for each of the four planets. All spectra were observed in the Y band with the 12'' × 0432 slit. The exposure time was 60 s frame^–1, and we used an ABBA nod to subtract background, giving us a typical observing efficiency of 77%.

Table 2 shows the target-specific details for each observation. We include TOI 560b for comparison, with the full details for this target provided in Zhang et al. (2022c). All 4 nights of observations were affected to varying degrees by weather. During the observations of TOI 560b, 1430.01, and 1683.01, transparency was poor throughout the night. The 1683.01 half-night also suffered from poor seeing. When we began observing at the beginning of the night, the seeing was 1'', but it started rising around 12:44 UTC, eventually reaching values as high as 3''. The combination of a faint host star, clouds, and slit losses from poor seeing resulted in a low signal-to-noise ratio (S/N) for this target. For TOI 2076b, a combination of telescope problems, fog, and unexpected transit timing variations (TTVs; see Appendix C) prevented us from observing the first two-thirds of the transit. However, the remainder of the night had stable seeing and no detectable clouds.

Table 2. Keck/NIRSPEC Observations

Parameter/Planet	TOI 560b	TOI 1430.01	TOI 1683.01	TOI 2076b
Date (UTC)	2021-03-18	2021-09-21	2021-10-01	2022-03-16
Time (UTC)	12:06–16:07	05:04–09:58	10:31–15:51	12:06–16:07
S/N	150	160	67	190
Preingress baseline (hr)	0.7	0.7	2.4	0
Transit fraction observed (%)	100	100	100	31
In-transit duration (hr)	2.1	2.7	1.3	1.0
Postegress baseline (hr)	0.6	1.4	1.6	3.0

Download table as:  ASCIITypeset image

To analyze the data, we used the same pipeline described in Zhang et al. (2022c). Among other steps, this pipeline produces A-B subtraction images to remove the background, uses optimal extraction on the differenced images to obtain 1D spectra, fits a template to the spectra to obtain wavelength solutions, and runs molecfit to remove telluric absorption features. The pipeline calculates excess absorption during the transit as a function of time and wavelength, which can then be collapsed along either axis to give the average in-transit excess absorption spectrum or the band-integrated light curve.

In addition to fitting telluric features, molecfit also fits for the FWHM of the line-spread profile, but it is not stable enough to reliably fit the FWHM for every spectrum. Instead, we run it on every spectrum and take the median of the FWHMs across the night for every target except TOI 1683, which has a low S/N. We obtain 3.4 pixels, corresponding to R = 32,000 at the position of the helium line. We then rerun all the molecfit fits with the FWHM fixed to 3.4 pixels. This value is close to what we obtain for each individual night: 3.65 ± 0.38 (TOI 560), 3.36 ± 0.17 (TOI 1430), and 3.40 ± 0.12 (TOI 2076), where the reported error is the standard deviation of the FWHM across the night.

Figure 1 shows the stellar spectrum in the vicinity of the helium line both before and after correction for tellurics by molecfit. The strongest stellar line is Si i at a rest-frame vacuum wavelength of 10830 Å; the second strongest line is the metastable helium line at 10830 Å. For TOI 1430 and 2076, the helium line was far from any tellurics. For TOI 1683, telluric absorption begins to pick up only at the reddest edge of the helium line, 20 km s⁻¹ from the center. As we will see, the observed planetary absorption signal is blueshifted, making the telluric absorption even less important. It is only for TOI 560 that a telluric line falls right on top of the helium line. Due to extraordinary luck, the telluric absorption was minimal because this night happened to be exceptionally dry; molecfit reports a precipitable water vapor column of 0.15 mm, compared to 0.9 mm for TOI 1430, 1.0–1.4 mm for TOI 2076, and ∼1 mm for TOI 1683. We discuss the telluric correction for TOI 560 in greater detail in Zhang et al. (2022c) and are confident that it does not greatly affect our conclusions.

Zoom In Zoom Out Reset image size

2.3. X-Ray Observations

Figure 2 shows the excess absorption from each of the four planets as a function of time and wavelength relative to the average flux outside of the white-light transit. Figure 3 shows the average in-transit excess absorption spectrum in the planetary rest frame, while Figure 4 shows the band-integrated light curve in a 3 Å bandpass centered on the main helium peak.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The four absorption signals share both striking similarities and striking differences. Despite having different radii, orbital periods, ages, and masses, all four planets show an average in-transit absorption depth of 0.7%–1.0%. None of the planets show an evolving velocity that tracks the radial acceleration of the planet. This might cause one to suspect that stellar activity is the origin of the putative planetary signals, but the band-integrated light curves strongly argue for a planetary origin; they decrease at the white-light ingress and recover by the white-light egress for all planets except TOI 2076b. One final common feature of all of our observations is unexplained variations in the light curve, such as the rise prior to ingress for TOI 560b, the dip just after ingress for TOI 1430.01, the dip well before ingress for TOI 1683.01, and the decline for TOI 2076b 3 hr after egress. We attribute these features to stellar activity because our targets are young stars, and the stellar helium line is a well-known tracer of chromospheric activity. It is also conceivable that some of this variability is due to absorption from the escaping planetary atmosphere, which may form extended and complex structures (see, e.g., the 3D simulations of McCann et al. 2019).

The differences between the planets are as interesting as the similarities. One of the most puzzling features of the TOI 560b absorption signal, and the hardest to explain with our 3D models, is the redshift (Zhang et al. 2022c). Young stars are expected to have a strong stellar wind pushing the outflow away from the star, which is toward the observer during transit, and previous helium observations of giant exoplanets—including HD 189733b (Salz et al. 2018; Guilluy et al. 2020; Zhang et al. 2022a), WASP-107b (Allart et al. 2019; Kirk et al. 2020), WASP-69b (Nortmann et al. 2018), GJ 3470b (Palle et al. 2020), HD 209458b (Alonso-Floriano et al. 2019), and HAT-P-11b (Allart et al. 2019)—saw a blueshift. The sole exception is HAT-P-32b (Czesla et al. 2022). Our new observations show that the redshift we observed was unusual, as our other three planets all show net blueshifts. Also potentially exceptional with TOI 560b was the time evolution of the absorption; in Figures 2 and 4, it is evident that absorption peaks after mid-transit. Object TOI 1430.01 does not show this characteristic, but more data are required to evaluate whether TOI 1683.01 and 2076b do.

If TOI 560b has its quirks, the other planets are no less peculiar. Object TOI 1430.01 shows what appears to be preingress absorption. Without any detected velocity evolution matching the planet's orbit, it is difficult to determine whether this preingress absorption might be due to stellar variability or a combination of planetary absorption and stellar variability. Real preingress absorption is not implausible, as it has been seen for HAT-P-32b (Czesla et al. 2022) and can be produced by an up-orbit stream (e.g., Lai et al. 2010; Matsakos et al. 2015; McCann et al. 2019). We are not aware of any published 3D hydrodynamic simulations that predict strong preingress helium absorption, but MacLeod & Oklopčić (2022) predicted weak preingress absorption in their simulation with an intermediate-strength (∼10× solar) stellar wind. Incidentally, in their simulation with a strong stellar wind (∼100× solar), the helium absorption has considerably less velocity evolution than the planet itself, reminiscent of our observations. Finally, TOI 2076b is unique among our four planets in showing extensive absorption for at least half an hour after white-light egress, with the band-integrated light curve only declining 50 minutes after egress.

Object TOI 2076b is also unique in showing both a clear secondary peak and a primary peak, both blueshifted. The peak ratio is 0.40 ± 0.05, and there is a distinct valley in between the two peaks. The peak ratio is 0.125 in the case of a perfectly optically thin outflow and 1 for a perfectly optically thick outflow; thus, the observed peak ratio indicates an effective optical depth in between these two extremes. We do not detect a secondary peak for the other three planets, indicating that their outflows are likely closer to the optically thin limit. There are relatively few detections of secondary peaks in published helium observations of giant planets (e.g., HD 189733b, HAT-P-32b), and to our knowledge, no other planet has such a prominent secondary peak separated from the main peak by such a prominent valley. This is partly due to the high S/N of our detection; with higher-S/N data, HAT-P-32b (Czesla et al. 2022) may turn out to have a similar peak ratio and valley depth. Salz et al. (2018) reported a similar peak ratio of 0.36 ± 0.03 for HD 189733b, corresponding to an optical depth in the main peak of 3.2, although the peaks are not as clearly separated.

Coincidentally, Gaidos et al. (2023) observed the same transit of TOI 2076b from the same mountain using the high-resolution InfraRed Doppler spectrograph on the 8.2 m Subaru telescope. They obtained a lower S/N than we do, with much higher correlated noise. Their results are consistent with ours: 1% absorption in the helium line during transit, with some postegress absorption. However, they tentatively interpreted their results as stellar activity rather than planetary absorption. Given that we see planetary absorption of similar amplitude and width in every planet of our survey, we disagree with this interpretation, but more transits will be needed to conclusively disambiguate between these two possibilities.

Table 3 shows the quantitative properties of the outflow: peak absorption, redshift, peak ratio, and equivalent width. Defining these properties is not trivial because the line profile is unknown. We attempted to fit two Gaussians to the excess absorption spectrum, one for each absorption peak, but the fit was not good and would have resulted in misleading values if used for quantitative purposes. Instead, we take a data-driven approach. The peak absorption is taken from the highest data point. The redshift is similarly taken from the position of the highest data point and assigned an error of 1 pixel (TOI 1683.01) or half a pixel (all other planets). The peak ratio is calculated by dividing the mean of the three points closest to the secondary peak by the mean of the three points closest to the primary peak. The primary peak is identified by the point of highest absorption, while the secondary peak is assumed to be blueshifted or redshifted by the same amount as the primary peak. Because we do not know how much of the absorption at the secondary peak is actually due to highly blueshifted gas absorbing at its primary peak, these ratios are at all upper limits. For TOI 2076b, we also obtain a lower limit by subtracting off the depth of the valley. The equivalent width is calculated from the light curves in Figure 4. For each planet except TOI 2076b, we reject the two lowest in-transit points to account for outliers and average the next four lowest in-transit points. For TOI 2076b, for which we have limited in-transit coverage, we average the two lowest in-transit points. These are next to each other and are the closest points to mid-transit, giving us increased confidence that the points are not outliers. Finally, we calculate the width of the absorption by dividing the equivalent absorption by the height. The width is the most consistent metric and is within 10% of 1.1 Å for all planets.

Finally, Figure 5 shows the maximum excess absorption in units of scale height and compared to the EUV flux. This way of comparing planets has been used in many previous papers (e.g., Nortmann et al. 2018; Orell-Miquel et al. 2022). The uncertainties are large because the planetary masses are not well known, but the mini-Neptunes tend to be higher on this metric than the gas giants.

Zoom In Zoom Out Reset image size

Our next step is to translate the magnitude of the measured helium absorption signal into a present-day mass-loss rate for each planet. Unfortunately, this rate is not directly derivable from the data. The observed absorption comes from a trace species, helium in the triplet ground state, comprising less than a millionth of the atoms in the outflow and whose abundance is determined by photoionization and collisional de-excitation, among other processes, and therefore changes with distance from the planet. The properties of the outflow depend on the metallicity of the atmosphere (e.g., Zhang et al. 2022c), which could plausibly range from 1× to 100× solar for mini-Neptunes. The stellar wind and planetary and interplanetary magnetic fields influence the outflow geometry, which in turn alters the magnitude and shape of the absorption signature. The wind and magnetic field have highly uncertain strengths, especially because the wind is highly variable. For these reasons and others, we have previously struggled to match our Lyα absorption data for HD 63433b/c (Zhang et al. 2022b) and the helium absorption data for TOI 560b (Zhang et al. 2022c) with either the 1D hydrodynamics code TPCI (Salz et al. 2015a) or the 3D hydrodynamics code Microthena (Wang & Dai 2018). Without an accurate physical model of the outflow, it is impossible to derive a precise mass-loss rate.

Instead, we calculate order-of-magnitude estimates of the mass-loss rate using two methods. The first method, introduced in Zhang et al. (2022c), estimates the mass of metastable helium crossing in front of the star from the maximum equivalent width of the absorption, assuming the outflow is optically thin. It then converts the mass of metastable helium into total mass, assuming that 10⁻⁶ of the helium atoms are in the metastable state, and that 10% of all atoms are helium while the other 90% are hydrogen. The 10⁻⁶ ratio is typical of a K-type host star, as shown by our 1D and 3D simulations of HD 63433b/c and TOI 560b and the independent models of Oklopčić (2019; their Figure 2). To obtain a mass-loss rate, this total mass is divided by a timescale on which the mass disappears from view, which we take to be the radius of the star divided by a typical sound speed (10 km s⁻¹).

The second method we use to estimate the mass-loss rate is to put the Parker wind model of Oklopčić & Hirata (2018) in a Bayesian framework and run nested sampling retrievals on our data sets. The model requires the stellar spectrum, which we reconstruct using the methodology of Zhang et al. (2022c). Importantly, the EUV reconstruction is highly uncertain because the stellar EUV is only measurable from space, and no space telescopes currently have EUV capabilities (France et al. 2022). For example, France et al. (2022) compared four different methods of reconstructing the EUV spectrum for Proxima Centauri and found that they were discrepant by 3–100×, depending on the wavelength; an order-of-magnitude difference in EUV flux can result in a 3× difference in helium absorption depth (Oklopčić 2019).

Leaving aside the EUV uncertainty, the Parker wind forward model has three free parameters: the log of the mass-loss rate, the temperature, and the velocity offset, all of which are given uniform priors. The first two parameters are required to compute the Parker wind, while the third is meant to account for the moderate redshifts and blueshifts seen in the data, which could be due to a combination of the nonspherical nature of photoevaporation (only the dayside is irradiated), the stellar wind, and a small eccentricity of the planetary orbit. The Parker wind model outputs optical depth as a function of impact parameter and wavelength, which we then combine with the stellar radius, stellar limb-darkening coefficients, and transit parameters to derive a prediction of the excess absorption as a function of wavelength and time. The log-likelihood is computed by comparing this prediction to the observed absorption and error, and the nested sampling code dynesty (Speagle 2020) repeatedly calls the log-likelihood function to compute the posterior distribution.

This model assumes that the outflow is isothermal and spherically symmetric and that it consists only of hydrogen and helium in a 10:1 ratio. These assumptions are of dubious validity. For TOI 560b, our fiducial 3D simulations showed an outflow that rose in temperature from 1000 K at launch to almost 8000 K at the edge of the interface with the stellar wind, while the 1D solar metallicity simulation showed an outflow that declined from 10,000 K at maximum to 3000 K at 15 R_p . The TOI 560b outflow is also not symmetric, as can be seen from Figures 2 and 4, and neither is the TOI 1430.01 outflow, if the preingress absorption is not in fact stellar variability. Lastly, mini-Neptunes can plausibly have very high atmospheric metallicities. For all of these reasons, we do not expect the Parker wind model to be much better at estimating the mass-loss rate than the order-of-magnitude method.

All timescale estimates are shown in Table 4. These were obtained by assuming that 1% of the planetary mass is in the envelope, with the planetary mass being taken from Table 1. After estimating the timescale using both methods, we computed the mass-loss efficiency, as defined by η:

$$\begin{eqnarray}&&\displaystyle \frac{{dm}}{{dt}}=\displaystyle \frac{\eta }{4}\displaystyle \frac{{R}_{\mathrm{XUV}}^{3}{L}_{\mathrm{XUV}}}{{{GM}}_{p}{a}^{2}.}\end{eqnarray} \tag{ 1 }$$

Table 4. Estimated Timescale to Lose 1% of Total Planet Mass

Planet	OOM (Gyr)	Parker (Gyr)	${T}_{\mathrm{Parker}}({kK})$	ηOOM	${\eta }_{\mathrm{Parker}}$
TOI 560b	0.9	${0.12}_{-0.04}^{+0.05}$	9.8 ± 0.7	0.5	3.9
TOI 1430.01	0.5	0.10 ± 0.03	6.7 ± 0.3	0.5	2.8
TOI 1683.01	0.7	${0.61}_{-0.27}^{+0.54}$	6.7 ± 1.4	0.2	0.3
TOI 2076b	0.6	${0.70}_{-0.11}^{+0.16}$	5.0 ± 0.4	0.4	0.3

Note. Despite the naming, both methods should be considered order-of-magnitude estimates only (see Section 4).

Download table as:  ASCIITypeset image

To roughly calculate R_XUV, we follow Wang & Dai (2018) in assuming that the EUV photosphere is at ρ = 10⁻¹³ g cm⁻³. Assuming that the atmosphere is isothermal at the equilibrium temperature, and that the white-light radius corresponds to P = 100 mbar, we can derive the radius of the EUV photosphere (which we assume to be similar to the XUV photosphere):

$$\begin{eqnarray}&&{\rho }_{\mathrm{phot}}=\displaystyle \frac{P\mu }{{k}_{B}{T}_{\mathrm{eq}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{{{GM}}_{c}\mu }{{R}_{p}{k}_{B}{T}_{\mathrm{eq}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{EUV}}=\displaystyle \frac{{R}_{p}}{1+{\beta }^{-1}\mathrm{ln}({\rho }_{\mathrm{EUV}}/{\rho }_{\mathrm{phot}})}.\end{eqnarray} \tag{ 4 }$$

Adopting this procedure, we find XUV radii 28%–43% higher than the white-light radii. Changing the pressure of the white-light photosphere by a factor of 10 only changes this result by a few percent, but violating the isothermal assumption would result in larger changes.

Using Equation (1), we calculate the mass-loss efficiency for each planet and report it in Table 4. Just like the mass-loss timescales, these efficiencies should not be trusted to better than a factor of a few, and efficiencies greater than 1 are particularly unrealistic. However, taking into account the large uncertainties, our observational values are sensible. Caldiroli et al. (2022) used the 1D hydrodynamic code ATES to calculate the efficiency of photoevaporation from a variety of planets with different gravitational potentials and XUV irradiations. For planets similar to ours in gravitational potential and XUV irradiation, they found η ≈ 0.7, comparable to the values we computed from our observations. Our observations allow us to rule out scenarios in which the planets have a very high gravitational potential, corresponding to very low efficiencies (∼3 × 10⁻⁴).

Our detection of escaping helium from all of the first four planets in our survey strongly suggests that most mini-Neptunes have hydrogen/helium atmospheres, disfavoring the water-world hypothesis for their low densities. The high mass-loss rates we infer, which are sufficient to strip the envelope in hundreds of megayears, show that the radius gap separating super-Earths and mini-Neptunes must evolve over a planet's lifetime. However, they do not necessarily disprove the primordial origin of the radius valley, as Lee & Connors (2021) found that photoevaporation would fill in a primordial radius valley instead of deepening it. Further work is necessary to investigate the consistency of our observations with the primordial radius valley hypothesis.

4.1. Photoevaporation or Core-powered Mass Loss?

4.2. On Masses and Eccentricities

We surveyed four nearby mini-Neptunes orbiting young K-type stars and detected helium absorption from all of them. The helium absorption signals show some consistent properties across the four planets, but each planet has unique features as well. Our simple models suggest that the observed outflow velocities and temperatures appear to be consistent with expectations for photoevaporation, with an inferred mass-loss rate within a factor of a few of hydrodynamical predictions. The mass-loss timescales for the observed outflows are on the order of several hundred megayears. This suggests that these mini-Neptunes will all likely lose their hydrogen/helium envelopes, evolving into super-Earths.

If we wish to develop a more detailed understanding of mass-loss processes on these planets, there are several promising paths forward. On the observational side, additional helium observations of these planets are needed in order to map out the kinematic structures of their outflows in more detail, mitigate the effect of stellar variability on the measured helium signals, and characterize the magnitude of any time variability in the outflow properties. Observations with Lyα could be used to probe neutral hydrogen absorption far from the planet and at high velocities, providing complementary constraints on the outflow properties, although many factors could result in the nondetection of even a strong outflow. If we could measure the planetary masses for TOI 1430.01, 1683.01, and 2076b via radial velocities or TTVs, it would also greatly facilitate our theoretical understanding of their outflows because it would help reduce the dimensionality of the model parameter space and minimize degeneracies with other model parameters, such as the assumed EUV flux and atmospheric composition. By expanding the sample of young mini-Neptunes with outflows detected in the helium line, we can also obtain a better statistical picture of their outflow properties. This will allow us to search for correlations between specific features (like redshifted/blueshifted absorption) and other planet properties. We are currently pursuing all of these approaches. For example, we have HST programs to look for Lyα absorption from TOI 560b and 560c, with plans to write more proposals if these programs return promising data.

On the theoretical front, there is also much work to be done. In our previous studies of the mini-Neptunes HD 63433b/c (Zhang et al. 2022b) and TOI 560b (Zhang et al. 2022c), we modeled our observations with the 1D hydrodynamic code TPCI (Salz et al. 2015a) and the 3D hydrodynamic code Microthena (Wang & Dai 2018) and found that neither code was able to to fully reproduce the observed signals. We encourage other teams with 3D hydrodynamic codes to try modeling the four outflows reported here. It is also important to consider whether it is realistic to expect any model to match the observations, given the uncertainties in the composition of the outer atmosphere, the amount of turbulence in the outflow, the XUV flux, the planetary and interplanetary magnetic fields, the velocity of the stellar wind, and the potentially large time variability in all of these quantities due to stellar activity. If this variability turns out to be a limiting source of uncertainty, expanding the sample of mini-Neptunes with helium measurements would help to mitigate against these effects by enabling broader population-level comparisons between data and model, as well as an empirical search for correlations between outflow properties and planetary or stellar parameters.

The helium data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

Based on observations obtained with XMM-Newton (observation ID 0882870701 for TOI 1430, 0882870501 for TOI 1683), an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. The grant number is 80NSSC22K0742.

Funding for the TESS mission is provided by NASA's Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and the TESS Science Processing Operations Center. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Software: numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), dynesty (Speagle 2020), SAS, exoplanet (Foreman-Mackey et al. 2021a, 2021b), PyMC (Salvatier et al. 2016), theano (Theano Development Team 2016), celerite2 (Foreman-Mackey 2018).

All TESS data used in this paper can be found in MAST: 10.17909/b5pm-bb37.

In Figure 6, we show the phase‐folded TESS light curves for TOI 1430.01 and 1683.01, which we analyzed ourselves because they are not yet published elsewhere in the literature.

As part of this paper, we analyzed the X-ray measurements of TOI 1430 and 1683 taken by XMM-Newton's three EPIC cameras. The light curves are plotted in Figure 7. The measured spectrum of TOI 1430 and the theoretical model that best explains it are plotted in Figure 8. We fit a one-component model for TOI 1683 and a two-component model for TOI 1430 because the S/N is far higher for the latter than for the former.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The physical parameters we obtained from our Markov Chain Monte Carlo fits to the X-ray data are shown in Table 5.

Table 5. Model Parameters for XMM-Newton Data

Parameter	Value (TOI 1430)	Value (TOI 1683)
Metallicity	${0.29}_{-0.16}^{+0.03}$	${0.13}_{-0.07}^{+0.35}$
kT1 (keV)	${0.16}_{-0.02}^{+0.31}$	${0.27}_{-0.03}^{+0.05}$
EM1 (cm−3)	${8.1}_{-2.1}^{+6.6}\times {10}^{50}$	${17}_{-11}^{+7}\times {10}^{50}$
kT2 (keV)	${0.68}_{-0.04}^{+0.06}$
EM2 (cm−3)	${14}_{-4}^{+8}\times {10}^{50}$
Fluxa (erg s−1 cm−2)	${1.37}_{-0.06}^{+0.15}\times {10}^{-13}$	${0.32}_{-0.04}^{+0.05}\times {10}^{-13}$

Note.

^a Derived, not a fit parameter. For the range 5–100 Å (0.124–2.48 keV).

Download table as:  ASCIITypeset image

When we proposed for our TOI 2076b helium observations, we used the Hedges et al. (2021) ephemeris to predict transit times, but Osborn et al. (2022) subsequently discovered TTVs in the three-planet system. As a result, the Hedges et al. (2021) ephemeris is far less accurate than its error bars imply. Since Osborn et al. (2022), TESS has reobserved the system and captured three more transits of b, one more transit of c, and one more transit of d. The three transits of b are the three immediately following our Keck/NIRSPEC transit.

Accurately predicting the transit time for our Keck observation is critical to the interpretation of the helium outflow. Ideally, we would use a Bayesian method to model the TTVs of all three planets simultaneously and the ensemble of TTV models consistent with the data to generate predictions of the transit time. However, Osborn et al. (2022) reported that none of their three models are robust; one relies on the dubious assumption that the planets are not in resonance, and the other two arrive at either implausibly low or implausibly high values for the planetary masses, depending on our choice of priors. To estimate the transit time, we therefore consider three possible approaches. First, we fit a linear ephemeris to the three transits in the recent TESS sector and extrapolate back one epoch; second, we fit a linear ephemeris plus a sinusoidal TTV model to all transit times for b; and third, we fit a full TTV model using the ensemble of measured transit times for all three planets.

All three approaches require us to derive transit times from TESS data. We do this using the exoplanet package (Foreman-Mackey et al. 2021a, 2021b) by adapting their tutorial entitled "Simultaneous Fitting of a Transit with Stellar Variability." Specifically, we download the SAP fluxes and model the stellar variability using the quasiperiodic RotationTerm kernel from the celerite2 package (Foreman-Mackey 2018), which is a sum of two simple harmonic oscillator terms: a primary term at the stellar rotation period and a less coherent (lower quality factor, Q) secondary term at half the stellar rotation period. Unlike in the tutorial, we set the priors on log(Q) and log(ΔQ) to be Gaussians with a mean of 12 and standard deviation of 5, instead of a mean of zero and standard deviation of 2, because visual examination of the light curve shows highly coherent rotation across the entire month for all sectors of observation.

In addition, unlike in the tutorial, we fix all stellar and transit parameters except the transit times to the median values derived by Osborn et al. (2022). This is for two reasons: first, Osborn et al. (2022) included non-TESS photometry for all planets, so their values may be more robust, and second, we wanted to minimize the number of free parameters to avoid getting caught in a local minimum, which sometimes happened when we experimented with additional parameters. Our final model includes seven free parameters describing the Gaussian processes model (mean, white-noise level, rotation period, rotational variability σ_rot, log(Q), log(ΔQ), and relative amplitude of secondary term), in addition to individual transit times for each planet (14 transits in total, five in the new sector). We use PyMC (Salvatier et al. 2016) to perform No-U-Turn Hamiltonian Monte Carlo sampling with two chains with 1500 tune-up steps and 1000 sampling steps each and use the built-in tools to check for convergence. Table 6 shows the transit times we measure from our Gaussian processes fit; non-TESS times from Osborn et al. (2022) are added for completeness.

When we use only the new sector of TESS data, compute a linear ephemeris from its three transits of b, and extrapolate one epoch back to the Keck observation epoch, we obtain a transit time of 2,654.9837 ± 0.0015 (BJD_TDB – 2,457,000). When we compute all transit times of all planets from all TESS data, add in the non-TESS timings reported by Osborn et al. (2022), and perform a linear least-squares fit of the b transit times to a linear ephemeris plus a sinusoid, we obtain a transit time of 2,654.9841. The sinusoid is assumed to have a superperiod of 1/∣2/P_c − 1/P_b ∣ = 713.74 ± 0.14 days, consistent with the 713.1 ± 2.7 day superperiod found by Osborn et al. (2022). For our third method, we combine the N-body code TTVFast (Deck et al. 2014) with the nested sampling code dynesty (Speagle 2020) to run two full TTV models on all transit times of all planets, following the methodology of Hadden & Lithwick (2016). Briefly, we assume zero mutual inclination and have five parameters per planet (mass, ecos(ω), esin(ω), period, and mean anomaly). The mean anomaly is reparameterized as the epoch zero transit time to minimize the degeneracy with period. In the "low-mass" TTV model, we adopt a log-uniform prior on the masses from 1 to 20 M_⊕ and a uniform prior on the eccentricities from 0.001 to 0.2; in the "high-mass" TTV model, we adopt a uniform prior on the masses and a log-uniform prior on the eccentricities, with the same bounds. After nested sampling finishes, we draw randomly from the samples, taking into account the weights, and rerun TTVFast to predict the transit time for b at the desired epoch. We obtain 2,654.9830 ± 0.0024 for the low-mass model and 2,654.9837 ± 0.0019 for the high-mass model. Since all of our estimates for the transit time are within 1.6 minutes of each other, we are confident that we correctly calculated the transit time despite the TTVs. We adopt as our final answer the result of the first method, 2,654.9837 ± 0.0015. Unfortunately, this timing is 51 minutes earlier than what we expected when we proposed the Keck observations. In combination with the telescope and weather problems, this meant that we were unable to observe the first two-thirds of transit.

In running the TTV models, we had hoped to compute not just the transit time corresponding to the Keck observations but also the masses of the planets. Unfortunately, our fits were underconstrained, and we found that our masses were highly sensitive to our choice of priors. The best-fit TTV models using both priors are shown in Figure 9. We obtained ${1.8}_{-0.6}^{+1.9}$, ${3.1}_{-1.0}^{+2.5}$, and ${2.3}_{-1.0}^{+2.8}\,{M}_{\oplus }$ using the low-mass prior and ${11.3}_{-4.2}^{+5.1}$, ${11.4}_{-2.7}^{+2.5}$, and ${6.4}_{-4.3}^{+6.6}\,{M}_{\oplus }$ using the high-mass prior. The low-mass fit favors eccentricities of a few percent, while the high-mass fit favors eccentricities of a few tenths of a percent. The large uncertainties on the masses, together with the strong dependence on the choice of prior, mean that additional transit and radial velocity observations are needed to fully characterize this system.

Zoom In Zoom Out Reset image size

As a final step, we improved upon the ephemerides of Osborn et al. (2022) by deriving new periods and transit midpoints. For systems with TTVs, there is no single correct definition of period. Possible definitions include (1) the coefficient found from a linear least-squares fit of transit time versus epoch, (2) the coefficient found from a fit of transit time versus epoch of the form $T=PE+{T}_{0}+A\sin ({\omega }_{{\rm{s}}{\rm{u}}{\rm{p}}{\rm{e}}{\rm{r}}}PE)+B\cos ({\omega }_{{\rm{s}}{\rm{u}}{\rm{p}}{\rm{e}}{\rm{r}}}PE)$, and (3) the period implied by the position and velocity of the planet at a given time if the other planets did not exist. Definitions 1 and 2 are the same if the TTVs are perfectly sinusoidal and the observations evenly sample the superperiod. Definition 3 depends on the time chosen. For our TTV fits, we choose a time of 1743, just before the first transit of b. In Table 7, we present the periods and mean transit times that follow from definition 1.