key: cord-0554288-x6vbsivm
authors: Kasper, David; Bean, Jacob L.; Oklopvci'c, Antonija; Malsky, Isaac; Kempton, Eliza M.-R.; D'esert, Jean-Michel; Rogers, Leslie A.; Mansfield, Megan
title: Non-detection of Helium in the Upper Atmospheres of Three Sub-Neptune Exoplanets
date: 2020-07-25
journal: nan
DOI: nan
sha: cd0ec199e24efb6bec090f6c27ad3c097e3d1d5d
doc_id: 554288
cord_uid: x6vbsivm

We present a search for helium in the upper atmospheres of three sub-Neptune size planets to investigate the origins of these ubiquitous objects. The detection of helium for a low density planet would be strong evidence for the presence of a primary atmosphere accreted from the protoplanetary nebula because large amounts of helium are not expected in the secondary atmospheres of rocky planets. We used Keck+NIRSPEC to obtain high-resolution transit spectroscopy of the planets GJ1214b, GJ9827d, and HD97658b around the 10,833 Ang He triplet feature. We did not detect helium absorption for any of the planets despite achieving a high level of sensitivity. We used the non-detections to set limits on the planets' thermosphere temperatures and atmospheric loss rates by comparing grids of 1D models to the data. We also performed coupled interior structure and atmospheric loss calculations, which suggest that the bulk atmospheres (winds) of the planets would be at most modestly enhanced (depleted) in helium relative to their primordial composition. Our lack of detections of the helium triplet for GJ1214b and GJ9827d are inconsistent with the predictions of models for the present day mass loss on these planets at more than 5 sigma confidence. Higher signal-to-noise data would be needed to detect the helium feature predicted for HD97658b. We identify uncertainties in the EUV fluxes of the host stars and the lack of detailed mass loss models specifically for cool and metal-enhanced atmospheres as the main limitations to the interpretation of our results. Ultimately, our results suggest that the upper atmospheres of sub-Neptune planets are fundamentally different than those of gas giant planets.

The basic nature of the atmospheres of sub-Neptune size planets, whether they are primary or secondary, is currently a major question in the field of exoplanets. The bimodality of the small planet radius distri- Figure 1 . The targets of our project in the context of the broader population of closein exoplanets. Left: Mass-radius diagram showing solar system planets (black stars), exoplanets with precisely measured parameters (grey circles, data from Otegi et al. 2020) , and our targets (red circles). The lines are theoretical models for pure compositions of iron, silicate, water, and hydrogen (Seager et al. 2007) . Right: Frequency of closein planets as a function of planet size from Fulton & Petigura (2018) . The grey part of the line indicates the regions of parameter space with poor constraints. Kreidberg et al. 2014; Knutson et al. 2014; Benneke et al. 2019; Tsiaras et al. 2019) . Furthermore, even if we could determine that these planets' atmospheres were hydrogen-dominated this would still not be definitive proof that they were accreted from the primordial nebula because hydrogen can be outgassed from the interiors of rocky planets Rogers & Seager 2010a; Rogers et al. 2011; Chachan & Stevenson 2018; Kite et al. 2020) .

In this paper we pursue a new way to investigate the atmospheric compositions of the mysterious sub-Neptune size exoplanets. We aimed to detect helium in these planets' upper atmospheres (i.e., the thermospheres and exospheres) using transmission spectroscopy. The thermosphere of a planet is the outermost bound atmospheric layer, above which is the exosphere and the transition to interplanetary space. The upper atmospheres of close-in exoplanets are expected to be extended due to substantial ongoing atmospheric escape driven by absorption of UV flux from their host stars. The advantage of targeting the thermospheres and exospheres for sub-Neptune size planets is that they should extend to altitudes well above the aerosol layers that blocks transmission spectroscopy observations of their bulk atmospheres. That is, we shouldn't be blocked from viewing the upper atmospheres by the aerosols that have caused previous transmission spectroscopy observations to yield flat spectra. If helium is present in the atmospheres of these planets then it should be abundant in the thermospheres and exospheres because it is a relatively light species.

Our observations target the 10,833Å He triplet feature, which has recently emerged as a way to probe the escaping atmospheres of larger planets (Spake et al. 2018; Nortmann et al. 2018; Allart et al. 2018 Allart et al. , 2019 Mansfield et al. 2018; Salz et al. 2018; Kreidberg & Oklopčić 2018; Alonso-Floriano et al. 2019; Ninan et al. 2020; Kirk et al. 2020; Gaidos et al. 2020; Palle et al. 2020; Vissapragada et al. 2020; dos Santos et al. 2020 ). The He triplet feature arises from a metastable transition. It has a strong absorption cross section that compensates for the low density of the upper atmosphere and allows detection even at very low partial pressures.

Our idea is similar to searching for hydrogen gas in the upper atmosphere using Lyman α (Bourrier et al. 2017) , but with the advantage that a detection of helium would have a less ambiguous interpretation. The direct detection of helium in the atmosphere of a planet would be the smoking gun evidence for the presence of a primary atmosphere accreted from the primordial nebula. This is because no other formation mechanism (e.g., sublimation of ices or outgassing of rocky material) can give rise to the presence of large quantities of helium in a planetary atmosphere . Helium is thus actually a cleaner diagnostic of formation history than the now standard transmission spectroscopy scale height test (Miller-Ricci et al. 2009 ). Furthermore, the helium infrared triplet that we target doesn't suffer from absorption by the interstellar medium like Lyman α.

We present Keck+NIRSPEC transit spectroscopy observations of the sub-Neptune size exoplanets GJ 1214b (Charbonneau et al. 2009 ), HD 97658b Dragomir et al. 2013) , and GJ 9827d (Niraula et al. 2017; Rodriguez et al. 2018) to search for the 10,833Å He triplet feature. Of these, only GJ 1214b has a pre-viously published result on the He triplet. Using low resolution (R ≈ 500) archival IRTF/SpeX observations, Crossfield et al. (2019) were able to place a modest upper limit of 2.1 R p at 95% confidence for absorption in the line core. Our observations reach much higher sensitivity for this planet due to the higher spectral resolution and signal-to-noise, and the lower levels of systematic noise in our data. The adopted properties of our targets are given in Table 1 and placed in the context of the population of sub-Neptune size planets in Figure 1 . These planets are archetypal intermediate-size planets as they sit squarely in the middle of the degenerate mass-radius space where there are multiple plausible models for their internal structure (e.g., Adams et al. 2008; Rogers & Seager 2010b ). Our targets were chosen for having the highest expected signal-to-noise for planets in this part of parameter space (formally, planets having M p < 10 M ⊕ and R p > 2 R ⊕ ) and being observable from Mauna Kea.

GJ 1214b in particular has been the subject of intense scrutiny that has revealed an astonishingly high and opaque aerosol layer that is challenging to explain Morley et al. 2015; Gao & Benneke 2018; Adams et al. 2019) . HD 97658b has also been observed to have a featureless transmission spectrum (Knutson et al. 2014; Guo et al. 2020 ). These two planets remain the best targets for transmission spectroscopy observations in this part of the parameter space, and GJ 9827d is not far behind. TESS is not expected to find better planets for this study Kempton et al. 2018) .

The planets we observed differ in terms of their host stars and system multiplicity, thus providing an interesting comparison set. GJ 1214b orbits an M dwarf, while HD 97658b and GJ 9827d orbit early and late K dwarfs, respectively. Looking at planets orbiting different stellar types is important for this study because both photoevaporative mass loss and the population of the metastable level that gives rise to the 10,833Å He triplet feature depend on the high-energy irradiation received from the host star. Oklopčić (2019) have suggested that K stars in particular have the ideal balance between extreme-ultraviolet (EUV) and mid-ultraviolet (mid-UV) flux levels. Beyond the stellar hosts, no other planets are known in the GJ 1214 ) and HD 97658 systems. However, GJ 9827d is part of a system with two other transiting planets on interior orbits. The density of the closest-in planet (planet b) is consistent with being rocky (Prieto-Arranz et al. 2018; Teske et al. 2018; Rice et al. 2019) , which is highly suggestive of photoevaporative mass loss from initially volatile-rich Southworth (2011) worlds. GJ 9827c is between planets b and d, but its mass is poorly constrained.

The paper is laid out as follows. We describe our observations and data analysis in §2 and §3. Since we do not detect the He triplet feature for any of the planets we perform a suite of atmosphere ( §4) and interior structure ( §5) model calculations to interpret the upper limits from our data. We discuss the implications of our results in §6, and we conclude in §7 with a summary.

We observed transits of our targets using the NIR-SPEC instrument (McLean et al. 1998 ) on the Keck II telescope on Mauna Kea, Hawai'i. An observing log is given in Table 2 . Our program was conducted following the NIRSPEC hardware upgrades described by Martin et al. (2018) . We used the 0.432 × 12 slit, which delivers a nominal resolving power R ≈ 25,000. We used the NIRSPEC-1 (Y-band) filter setting without the optional THIN blocking filter to avoid fringing effects at the recommendation of the Keck instrument scientist. We took observations in a standard ABBA nodding sequence with a throw of 6 . Dark, flat, and arc-lamp calibration frames were taken at the beginning and end of the observation periods.

We observed one transit each of GJ 1214b and HD 97658b, and two transits of GJ 9827d. One additional planned transit observation of HD 97658b in April 2019 was lost due to a telescope shutdown caused by the COVID-19 pandemic. We observed our targets continuously starting well before transit ingress and ending well after transit egress except for minor interruptions due to telescope and software glitches and major interruptions due to bad weather during the HD 97658b transit. For the HD 97658b transit, the dome had to be closed for 97 minutes during ingress and then closed for the night just after mid-transit. Nevertheless, we obtained a series of spectra both before and during transit for HD 97658b, and the data quality for this target is useful due to the very bright host star (m J = 6.2) and the large collecting area of the 10 m Keck telescope. The weather conditions were generally clear for the other transits. All observations were taken at airmass values less than 2.0, with most at less than 1.5.

We analyzed the data for this project using custom routines that have been previously used to reduce ground-based low-and high-resolution spectroscopy data for atmospheric studies and radial velocity measurements (Bean et al. 2010a (Bean et al. ,b, 2011 . The routines were originally written in IDL and were ported to Python, version 3.7, for this project. In addition to our own data, we also reduced the Keck+NIRSPEC WASP-107b data described in Kirk et al. (2020) as a check of our pipeline. This planet was the first to show the helium feature, and the feature has been detected for this planet using HST +WFC3, CARMENES, and NIR-SPEC, with consistent results from all the instruments (Spake et al. 2018; Allart et al. 2019 ). The helium feature appears in two of the NIRSPEC echelle orders (70 and 71), and we limited our analysis to just these orders. The scripts used for this paper are available from the authors upon request.

Processing the calibration frames -Our NIRSPEC data reduction starts with creating a master dark image for the flat field exposures by averaging exposures taken with a cold blank in the filter wheel blocking the light path. The dark frames were obtained using the same exposure time as the flat fields. We then subtracted this master dark image from the flat field exposures to remove the signals from the bias (pedestal) and dark currents. The flat field exposures were then averaged to create a single master flat field image. We removed the spectroscopic signature from the master flat for each order considered by averaging the rows in the cross dispersion direction, applying a median filter to smooth the data, fitting a high order polynomial, and then dividing out the fit. We then normalized the flat by dividing by the median of all the illuminated pixels in an order. In most of the averaging steps we used iterative outlier rejection with cutoffs ranging from three to five standard deviations. We created an initial bad pixel mask from the flat normalization process that was later used in the spectral extraction.

Processing the science frames -The science frames were first processed by subtracting the closest-in-time opposite nod pair frame taken with the same exposure length. This simultaneously removed the bias and dark currents, and the sky background (both continuum and emission line). Telluric emission lines can be a pernicious problem for data taken with fiber-fed spectrographs to search for the helium feature because the source and background are scrambled and the spatial information is lost (e.g., Ninan et al. 2020; Palle et al. 2020 ). Our slit-fed observations sacrifice stability (see later in this section) but make accurate sky subtraction straightforward.

After subtracting the nod pairs we divided the science frames by a normalized flat field. Then we traced the orders and used an optimal extraction algorithm (Horne 1986 ) to extract 1D spectra from the images. Our optimal extraction algorithm uses a spatial profile weighting that is determined from the data and includes iterative identification and masking of bad pixels and cosmic rays. We used an aperture radius of 10 pixels around the cen- Figure 2 . Example telluric removal around the 10,833Å feature using Molecfit. The continuum normalized GJ 1214 spectrum is shown prior to (black) and post correction (blue). The Molecfit best-fit transmission curvethe correction function for telluric absorption -is shown in dashed green. The gray region is used to correct instrumental drift between frames and is notably free of significant telluric or planet features that could compromise the crosscorrelation technique.

ter of the spectral trace for the extraction, and we verified that the results did not depend on our choice for this parameter.

Wavelength calibration -To provide an initial wavelength calibration of the data, we used a second order polynomial to fit the positions of emission lines bracketing the 10,833Å feature along each order in nightly NeArXeKr arc-lamps exposures, with six lines in order 70 and five lines in order 71. We did this separately for the A and B nod positions but we found that this wasn't necessary because the slit is well aligned with the columns on the detector. As a reminder, even though these data are from ground-based observations, the wavelengths are measured in vacuo because the detector is in a cryo-vacuum dewar. Contrary to Kirk et al. (2020) , we did not find that our initial wavelength solution based on the arc lamp data had any obvious distortions. During this step we also confirmed that the data have a resolution R ≈ 25,000 by measuring the width of the arc lamp lines.

Removing telluric absorption lines -Following Kirk et al. (2020) , we ran Molecfit Kausch et al. 2015) separately for the two orders in each frame to remove telluric absorption lines. We did this before applying barycentric corrections to the wavelength scale because the telluric features appear in the rest frame of the observatory. We found that by-hand Figure 3 . Drift correction for GJ 1214b data. The cross-correlated shift in km s −1 for each frame from the first frame is plotted as a function of phase. Care was taken in picking the wavelength region to perform the crosscorrelation on to be near the helium feature but avoid potential contamination from residual telluric features and planetary features.

identification of clean telluric feature regions was necessary to optimize the Molecfit model parameters. The exact number of clean telluric features varied by target. For example, we fit 13 features in order 70 and 8 features in order 71 in the GJ 1214b dataset to optimize the model parameters.

We allowed the H 2 O abundance to vary and assumed the nominal values of the CH 4 and CO 2 abundances. Molecfit divided out the best-fit telluric model determined for each order and frame. An example spectrum before and after telluric subtraction for GJ 1214 is shown in Figure 2 . We found, as expected, that the coincidental overlap of telluric absorption features (in the Earth frame) and the helium triplet feature (in the system frame) was different for each planet and observation due to differences in the systemic velocities and the motion of the observatory along the line of sight.

Correcting for radial velocity shifts -After the telluric correction, we applied shifts to the wavelength scales of each frame to put them in the rest frames of the corresponding systems (i.e., accounting for the motion of the observatory relative to the barycenter of the systems). The adopted systemic velocities used for this correction are given in Table 3 . Note that there is a sign error for the radial velocity of the GJ 1214 system given by Charbonneau et al. (2009) . The correct systemic velocity should be +21.3 km s −1 (i.e., the star is moving away from the solar system, private communication from E. Newton). Following this step we normalized the con- tinuum in the data by fitting a line to relatively clean spectral regions within 15Å to either side of the helium feature and then dividing it out. We noticed that the time series of spectra for each observation do not match up perfectly after ostensibly putting them in the system rest frame, which we attribute to a combination of variations in the illumination of the entrance slit and instability of the spectrograph optics. We corrected for this effect by cross correlating the spectra against the first spectrum taken during each observation and then shifting the spectra according to the measured values. An example of the drift in the spectra during the GJ 1214b observation is given in Figure 3 . The typical rms of the offset in the timeseries is 1.3 km s −1 , which corresponds to 0.46 pixels.

Our final analysis step was to create transmission spectra of the planets by looking for excess absorption around the helium feature during the transits. We calculated a master out-of-transit spectrum by taking the weighted mean of the continuum normalized spectra with the data labeled according to the ephemerides given in Table 3 . For the targets of this study (i.e., the ones we observed) we assumed data during ingress and egress were in transit, but for WASP-107b we used the exact frames that Kirk et al. (2020) identified in their analysis for direct comparison. Note that the handful of frames that were obtained during telescope or software glitches were excluded from our analysis.

Individual transmission spectra for each planet were calculated as T = 1 − F ini /F out for each in transit spectrum. In this convention for the transmission spectrum the continuum is around zero and excess planetary absorption gives positive numbers. These individual transmission spectra were then shifted into the planet's rest frame according to the planetary orbital motion as calculable from Tables 1 and 3 before co-adding to make a single planetary transmission spectrum.

This data processing was performed independently for NIRSPEC orders 70 and 71. The two resulting transmission spectra for each data set were then resampled to a common wavelength grid and averaged to give a combined spectrum. Figure 4 shows the final averaged transmission spectrum for WASP-107b compared to the results of Kirk et al. (2020) , which were derived for order 70 only. We achieve excellent agreement with their detection of a large signal, which gives confidence in our pipeline. Figure 5 shows the measured transmission spectra for our targets including example forward models for helium absorption in the planets' upper atmospheres (the models are described in detail in §4). As mentioned in the Introduction, we do not detect the signature of helium absorption in any of our targets. The scatter in the data is generally consistent with the expectations from the photon-limited error bars except near the positions of telluric features where correlated residuals can be seen from imperfect subtraction. Similarly, prior to averaging, the data in order 71 also had lower signal-tonoise and larger residuals because Molecfit had more trouble fitting the tellurics in that order as compared to order 70. Some nights are also better than others in terms of telluric absorption. For example, the second night for GJ 9827d yielded higher signal-to-noise data and cleaner residuals from the telluric subtraction. We use these data to set limits on the helium absorption and the properties of the planets' atmospheres in the following sections.

As another test for helium absorption, we plotted the time series of data in the system barycenter frame relative to the master out-of-transit spectrum and searched for signals that were not perfectly centered on the transits. Spectral absorption from the upper atmospheres of transiting planets do not have to correspond to the white light transit of the bulk planet. For example, GJ 436b shows a long tail of neutral hydrogen absorption at Lyman α (Ehrenreich et al. 2015) , and Kirk et al. (2020) found that the helium absorption from WASP-107b began after ingress and continued after egress. Figure 6 shows the time-series data for GJ 1214b with the expected trace due to the planet's orbital motion. Neither GJ 1214b nor any of the other planets in our sample showed evidence for a signal not centered on the transit.

We used the obtained upper limits on the He 10,833Å absorption signal to constrain the atmospheric properties of GJ 1214b, HD 97658b, and GJ 9827d. Using methods similar to those described in Oklopčić & Hirata (2018) and Mansfield et al. (2018) , we modeled the upper atmosphere of each planet as a spherically symmetric hydrogen-helium envelope extending to altitudes of several planetary radii (up to the Roche radius of the system). The assumed atmospheric density profile was that of an isothermal Parker wind, in which the atmosphere is close to hydrostatic at low altitudes, but has a radial velocity component, increasing with altitude. An isothermal Parker wind model has two free parameters, the temperature of the upper atmosphere (thermosphere) and the density normalization, which can be linked to the atmospheric mass-loss rate and thus measures the total mass escaping from the planet per unit time.

We performed radiative transfer calculations for each model atmosphere to obtain non-local thermodynamic equilibrium populations (as functions of altitude) of neutral and ionized hydrogen, and helium atoms in the ground, the metastable (2 3 S), and the (singly) ionized state. For GJ 1214 and HD 97658, we use the spectral energy distributions of these stars obtained from the MUSCLES survey Youngblood et al. 2016; Loyd et al. 2016) . For GJ 9827, we reconstruct the EUV part of the spectrum using the observed Ly α flux (private communication from I. Carleo) and the scaling relations from Linsky et al. (2014) ; the spectrum at longer wavelengths is modeled after the MUSCLES spectrum of HD 85512, a star of the same spectral type as GJ 9827. More details on the radiative transfer methods, including various reaction rate coefficients used in our calculations, can be found in Oklopčić & Hirata (2018) .

Using the obtained 1D density profile of metastable helium and assuming a spherically symmetric planetary atmosphere, we computed a mid-transit transmission spectrum at wavelengths around 10,833Å for a grid of model atmospheres for each planet. Each grid spans a broad range of values in thermospheric temperature (2,000 -8,000 K) and mass-loss rates (10 6 -10 11 g s −1 ). Our nominal model grid assumed roughly solar composition (10% helium number fraction) because this is the expected composition of the wind (see §5 and 6). We also calculated model grids for GJ 1214b with the helium fraction by number ranging from 5% to 20%, i.e. from subsolar to supersolar values, to explore how changing the helium abundance in the wind would impact our results.

We compared each simulated spectrum to the data to determine which parts of model parameter space could be excluded by our non-detection of the helium feature. To do this we first transformed the transmission models into their equivalent excess absorption. Then we broadened the models to the R = 25,000 instrumental response and resampled them to the wavelengths of the combined (order 70 and 71) observed data set. We then identified the optimal window around the 10,833Å feature for each model by maximizing the cumulative distribution function of the model absorption (defined at each pixel point) over the standard deviation of the observed data (defined in a constant ±1.5Å region around the feature). The window was centered on the middle of the two stronger lines of the triplet feature to maximize the contrast. The observed data region to use for testing the models was thus chosen as a compromise between sampling the true spread in the data assuming the null hypothesis and retaining information in the relevant region assuming the contrary. Following the optimization of the test window, the corresponding (maximal) rejection of the model given the data was found by its relation to the cumulative distribution function. Figure 7 shows the map of the statistical deviation from our data for the GJ 1214b model grids. As expected, the strength of the triplet feature is strongly correlated with the helium abundance in the wind. Therefore our non-detection rules out a larger part of the parameter space for high helium abundances, but is correspondingly less constraining for lower helium abundances. Figures 8 and 9 show the same maps for solar composition model grids for GJ 9827d and HD 97658b, respectively. We compared our model grid to the second observation of GJ 9827d because it was the higher quality dataset of the two we obtained.

Notably, the contours for GJ 1214b in Figure 7 reverse direction at high mass loss rates. This is because at a fixed temperature and at relatively low mass loss rate, increasing the rate increases the overall density, and consequently, the density of the metastable helium, which causes stronger absorption. However, after reaching a certain mass-loss rate value, the atmospheric density becomes so high that the EUV photons (which ionize helium and thus populate the metastable state) cannot penetrate all the way to the bottom of the atmosphere, so the lowest atmospheric layers become more and more depleted of metastable helium as the product of the mass-loss rate, and thus density, increases. 

Fractionation between different chemical species could lead to atmospheric loss with a wind composition that is different from the bulk of the envelope. The planets in our sample are expected to be in the hydrodynamic mass loss regime where the large-scale loss of hydrogen drags heavier species along (Salz et al. 2016, and see §6) . Nevertheless, diffusive separation that takes place between the homopause (where the atmosphere is no longer expected to be well-mixed due to eddy diffusion) and the sonic point (where the flow composition becomes fixed) can reduce the amount of heavier species that are present in the wind. Helium is at risk of fractionation because it is four times heavier than the atomic hydrogen that makes up the bulk of the wind. For example, the energy-limited escape of hydrogen from modern Earth is incapable of pulling along helium and other heavy species (Catling & Kasting 2017) . Mass fractionation in a hydrodynamic wind is enhanced for lower temperatures and lower mass-loss rates. Hu et al. (2015) suggested that mass-dependent fractionation could lead to preferential loss of hydrogen and a corresponding enhancement in helium abundance in the bulk atmospheres of Neptune-size planets. If our targets are in a regime where the relative abundance of helium in their bulk was increasing, it would necessarily mean that the winds are actually depleted in helium, thus making our observation more difficult.

We performed coupled interior structure and atmospheric mass loss calculations to explore this phenomenon and provide further context for the nondetection of helium in the upper atmospheres of our targets. In this modeling, we consider a scenario where the planets consist of an Earth composition rocky core surrounded by a primordial hydrogen-helium envelope. By matching interior structure models to the observed characteristics of GJ 1214b, HD 97658b, and GJ 9827d, we can constrain the planetary envelope fractions and predict the hydrogen and helium mass loss rates as a boundary condition for the atmosphere models. Furthermore, we can show how the bulk envelope fractions and wind compositions of these planets evolve over the course of billions of years.

We performed these calculations using the Modules for Experiments in Stellar Astrophysics (MESA) code (v12778; Paxton et al. 2011 Paxton et al. , 2013 Paxton et al. , 2015 Paxton et al. , 2018 Paxton et al. , 2019 and the coupled interior structure and atmospheric mass loss models outlined in Malsky & Rogers (2020) . We instituted the regime of hydrogen and helium escape from Hu et al. (2015) to model the coupled thermal, massloss, and compositional evolution of hydrogen-helium envelopes surrounding our sub-Neptune size targets. More details of our calculations can be found in Malsky & Rogers (2020) .

When the incident EUV radiation on a planet is extremely high the escaping wind is transonic and increases in incident flux do not drive further escape. Instead, the deposited energy is converted into translational and thermal energy in the atmosphere. To account for this we included the reduction factor from Johnson et al. (2013) . This reduction factor reduces the mass loss rate when the incident EUV heating rate is above the critical heating rate of the planet. This results in relatively constant mass loss rates for highly irradiated planets even with decreasing EUV flux. Figure 10 . Hydrogen and helium mass loss rates (top row), as well as bulk envelope compositions (bottom row) for GJ 1214b, HD 97658b, and GJ 9827d. These models have initial envelope fractions of 0.022, 0.0063, and 0.003, respectively. At 10 Gyr, the models for GJ 1214b and HD 97658b had radii consistent with their currently measured values (see Table 1 ). We were unable to simulate GJ 9827d for 10 Gyr due to extreme envelope depletion. However, the model evolved to match the present-day radius of this planet by 3 Gyr. All models had homopause temperatures of 3,000 K and heating efficiency factors of η = 0.10.

We parameterized our grid of models as follows. The planetary masses, orbital separations, and host star characteristics were taken from Table 1 . For each planet, we ran a grid with varying values for the EUV heating efficiency factor (η), initial envelope mass fraction, homopause temperature, and masses. All models started with an initial solar composition (X = 0.74, Y = 0.24, Z = 0.02), and evolved from a 'hot-start' to 10 Gyr. After 6.0 Myr of evolution, hydrodynamic mass loss was turned on. All other model parameters are identical to the full model description in Malsky & Rogers (2020) .

We define the planetary radius to be at 1.0 mbar, roughly the pressure level corresponding to the observed transit radii ). MESA does not calculate atmospheric structure down to these low pressures. Therefore, we extrapolate from the outermost zone in our models, assuming an isothermal temperature profile and constant values for gravity and the mean molecular mass.

The planet envelope has constant elemental abundances throughout each zone modeled in MESA. Mass is lost via a wind at the homopause radius. We defined the homopause radius as the location where the H-He binary diffusion coefficient is equal to the eddy diffusion coefficient. In general, the homopause radius is approximately 25% to 60% larger than the transit radius, with larger differences for smaller-mass planets with lower surface gravity. Our assumption that the wind is launched at the homopause is somewhat inconsistent with the predictions of the location of the thermosphere for our targets by Salz et al. (2016) , who suggest that the τ = 1 level for high-energy photons is at radii < 1.1 R p . However, this shouldn't have a strong impact on our results because the sonic points for the planets all occur at much larger radii (> 3 R p ). We adopt a nominal homopause temperature of 3,000 K for all the planets to be consistent with the predictions of Salz et al. (2016) . Figure 10 shows the mass loss rates and bulk envelope evolution of three models, with radii matched for our individual targets. We find that the bulk atmospheric composition of GJ 1214b does not change significantly from its initial solar abundance when evolved with hydrodynamic mass loss. In order to achieve the observed radius, the GJ 1214b model requires an envelope approximately 0.81-0.96% of the total planet mass (depending on the system's uncertain age of 3 -10 Gyr, Charbonneau et al. 2009 ), which precludes the bulk composition of the envelope from becoming significantly enhanced in helium. This result is robust against variations in age, homopause temperature, and planet mass.

HD 97658b and GJ 9827d have small enough transit radii that the bulk compositions of their envelopes may change through hydrodynamic mass loss. In particular, HD 97658b may have a bulk envelope composition that is modestly enhanced in helium if it has been subject to fractionated mass loss over billions of years. For the model shown in Figure 10 , 56% of the envelope mass was lost, and the envelope helium number fraction increased from 0.075 to 0.168 over 10 Gyr. However, Henry et al. (2011) estimated the age of the HD 97658 system to be 7 Gyr using three different indicators, including their measurement of 38.5 days for the rotation period of the host star. Guo et al. (2020) found a shorter rotation period of 34 days and also an activity cycle of 9.6 yr, which both point to even younger ages. Therefore, the two times enhancement (depletion) of helium for the bulk (wind) of this planet is likely an overestimate because the fractionation is expected to ramp up toward later ages as the stellar high-energy flux and overall loss rate decrease.

Although GJ 9827d has likely lost a significant fraction of its envelope, we predict that the bulk envelope composition is at most mildly enhanced in helium. For the model shown in Figure 10 , GJ 9827d lost over 84% of its total envelope mass, and increased in helium number fraction over 3 Gyr to 0.117, corresponding to helium mass fractions going from 0.240 to 0.303. This is due to the fact that the planet is so highly irradiated that hydrogen and helium are lost in approximately equal proportion to their abundances in the bulk of the planetary envelope.

For models with similar radii to GJ 1214b, HD 97658b, and GJ 9827d, and a mass loss efficiency factor of η=0.10, we find average hydrogen mass loss rates of approximately 1.1 × 10 9 , 3.8 × 10 8 , and 4.2× 10 8 g s −1 , respectively. Helium mass loss rates for the same set are 3.2 × 10 8 , 1.8 × 10 8 and 1.8 × 10 8 g s −1 , respectively.

Numerical instabilities prevented some of the models with the lowest initial envelope mass fractions from evolving for the full 10 Gyr. In particular, the mass and radius of GJ 9827d are at the boundaries of the equation of state (EOS) tables in MESA, and are therefore at the limit of what we can simulate. A more detailed discussion of MESA's capability of modeling highly irradiated sub-Neptunes can be found in Malsky & Rogers (2020) .

The key question for this paper is, what can the nondetections of the helium feature tell us about the atmospheric properties of sub-Neptune size planets? Our coupled interior structure and atmospheric mass loss calculations indicate that the compositions of the bulk atmospheres of these planets shouldn't be significantly altered from their original state. Therefore, our targets will have helium abundances roughly corresponding to the solar value in their bulk and winds if sub-Neptune size planets accrete primary atmospheres during formation. While helium-enhanced winds would be more detectable, our observations for GJ 1214b and GJ 9827d were sensitive to upper atmospheres with solar composition for a wide range of temperatures and mass loss rates. Thus the question becomes, what temperatures and mass loss rates were expected for these planets? Salz et al. (2016) presented planet-specific simulations of atmospheric escape for GJ 1214b and HD 97658b using a coupled 1D radiative-hydrodynamical simulation code. They predict thermosphere temperatures of approximately 2,500 and 3,500 K and mass loss rates of 5 × 10 9 and 3 × 10 9 g s −1 for these two planets, respectively. Since these predictions are so similar, and because GJ 9827d is intermediate between GJ 1214b and HD 97658b in terms of irradiation, we take average values of T = 3,000 K andṀ = 4 × 10 9 g s −1 as nominal predictions for all the planets for the sake of simplicity. The models indicated by the solid orange lines and dots in Figures 5, 7 -9 , and 12 correspond to this nominal set of parameters. The red and magenta models correspond to changing the temperature or mass loss rate by an illustrative amount from this point in parameter space.

Note that the mass loss rates predicted by Salz et al. (2016) are higher than those predicted by our own calculations that were presented in §5 by about a factor of four. This is probably due to our adoption of a low efficiency factor, while Salz et al. (2016) self consistently account for the physics that is captured in this term. Higher mass loss rates would reduce the fractionation Figure 11 . Equivalent height of the absorbing atmosphere as a function of high-energy instellation for planets with helium feature detections and upper limits. The equivalent height of the absorbing atmosphere (δRp) is normalized by the scale height of the bulk atmosphere (Heq). This metric was originally propsed by Nortmann et al. (2018) . The black points with errors are measurements taken from the literature. The GJ 3470b data are from Palle et al. (2020) , and the rest are from the compilation of dos Santos et al. (2020) . The red points are measurements for our sub-Neptune targets. The arrows represent 90% confidence upper limits.

between hydrogen and helium, thus supporting our assumption that the winds of our targets should be roughly solar composition if they originally began with primary atmospheres.

As can be seen in the figures, the nominal predictions of Salz et al. (2016) and the surrounding parameter space can be ruled out for GJ 1214b and GJ 9827d at high confidence. The results for HD 97658b are less constraining due to a combination of an unfavorable planetto-star radius ratio, poor observing conditions limiting the signal-to-noise of the data, and the telluric contamination that is endemic to ground-based observations. Further observations of this target would be useful.

We place our measurements in the context of other planets with published helium detections and upper limits in Figure 11 . Following Nortmann et al. (2018) , this figure plots the "equivalent height" of the absorbing atmosphere normalized by the scale height of the bulk atmosphere as a function of the high-energy instellation of the planets. The measurements for our planets are tabulated in Table 4 . The equivalent height of the absorbing atmosphere measurements we quote are 90% confidence upper limits. The results for GJ 1214b and GJ 9827d do not follow the correlation that is seen for larger planets between the high-energy instellation and the normalized height of the helium absorption. This suggests that the upper atmospheres of our targets are fundamentally different than those of giant planets.

The lack of a detection and the limits we can set for three similar planets that orbit stars with a range of high-energy fluxes points to a problem with at least one fundamental assumption in our chain of hypotheses about the upper atmospheres of sub-Neptune planets. One possible explanation for the mismatch between the model predictions and the data is that the models are inaccurate. We previously showed that the combination of the Salz et al. (2016) and Oklopčić & Hirata (2018) models perfectly reproduced observations of helium absorption in the Neptune-sized exoplanet HAT-P-11b (Mansfield et al. 2018 ). However, the predictions of the temperature and escape rate in the Salz et al. (2016) model, and the population of the metastable level that gives rise to the helium feature in the Oklopčić (2019) model both depend on the input stellar spectra in the EUV range. The EUV flux is usually reconstructed from the observed X-ray or UV fluxes, which are not precisely established for our targets or HAT-P-11b ). Oklopčić (2019, see her Figure 7 ) has shown that reasonable variations in how the stellar EUV spectrum is reconstructed can lead to factor of three differences in the strength of the helium feature.

We tested the impact of changing the high-energy input stellar spectra in our radiative transfer calculations by almost an order of magnitude, which is the typical level of uncertainty in the reconstructed EUV spectral energy distributions (Youngblood et al. 2019; Drake et al. 2020) . Figure 12 shows the result of this test. For our fiducial model for GJ 1214b (T = 3,000 K, logṀ = 9.6), changing the flux at wavelengths shorter than 1,000Å by a factor of three relative to the nominal MUSCLES spectrum of GJ 1214, changes the depth of the predicted absorption signal by several percent. However, even in the low EUV flux case, which produces the weakest He absorption signal, the predicted signal is ruled out by the observations at more than 5σ confidence. Nevertheless, this test shows that good characterization of the high-energy fluxes of stellar host stars are essential for understanding atmospheric loss in the sub-Neptune regime.

Another reason why our models might be inaccurate is that the lower temperatures of our sub-Neptune planets could impact the composition of their thermospheres. Both the Salz et al. (2016) and Oklopčić & Hirata (2018) models assume that all the hydrogen is in atomic form. This is likely a safe assumption for HAT- b Determined by integrating our stellar models from 100 to 504Å. We conservatviely estimate that these values are accurate to within a factor of three. Figure 12 . The transmission spectrum of GJ 1214b compared to models using different stellar EUV flux levels. The data for GJ 1214b are plotted as grey points with errors. The black dashed line indicates the standard deviation of the data around the He feature. The solid orange model is created utilizing the fiducial model (T = 3,000 K, logṀ = 9.6) and the nominal stellar EUV flux. The dashed line model was created utilizing a factor of three less EUV flux. The dotted model was created utilizing a factor of three more EUV flux. All three models are rejected at more than 5σ significance. P-11b, where the thermosphere temperature was predicted to be ∼7,000 K at the radial distance from the planet probed by the helium feature. In contrast, the thermosphere temperatures of our targets (∼3,000 K) are around the point that a significant amount of hydrogen could be in its molecular form. The presence of molecules will reduce the escape rate because heavier chemical species will have lower thermal velocities for a given temperature. Salz et al. (2016) describe a test calculation for GJ 1214b where they include molecule formation. They quote H 2 abundances of 10 -25% throughout the thermosphere of the planet, and corresponding reductions in the mass loss rate of 15%. This level of uncertainty in the mass loss rate is not significant for our work. Also, it is not clear if their calculations include the impact of molecular diffusion, which will cause the heaver H 2 to preferentially settle out, thus reducing its impact on the atmospheric escape. Nevertheless, more sophisticated models are needed to understand atmospheric mass loss in the low temperature and low surface gravity regime of sub-Neptune size planets.

Finally, the possibility that these planets have metalrich atmospheres cannot be ruled out. Such atmospheres would have lower helium abundances, which reduces the strength of the helium feature (see Figure 7) . Also, having more metals in the thermosphere could increase the cooling efficiency, thus leading to lower temperatures that will reduce the mass loss rate. We are not aware of any papers that study the details of atmospheric mass loss specifically for thick, metal-rich atmospheres. Hopefully our observational results motivate work in this area.

To summarize, we set upper limits on the presence of the helium infrared triplet in high resolution transmission spectra of the sub-Neptune size planets GJ 1214b, HD 97658b, and GJ 9827d. The nominal predictions of models for GJ 1214b and GJ 9827d are inconsistent with our non-detections at high significance if we assume these planets were born with primary atmospheres. These findings are somewhat robust (more so for GJ 1214b than GJ 9827d) to uncertainties in the expected temperatures, mass loss rates, and helium abundances of the planets' winds, and the EUV fluxes of their host stars. The results for these two planets are also inconsistent with the emerging trend of the detection of large helium features in giant planets with modest to to large high-energy instellation. The data quality for HD 97658b was insufficient to detect the nominal model predictions due to the observations being obtained in poor conditions. More sophisticated models for the atmospheric evolution of cool, metal-rich atmospheres are needed to interpret these results. More precise charac-terization of the stellar EUV fluxes would also be valuable.

Atmospheric Evolution on Inhabited and Lifeless Worlds

VizieR Online Data Catalog, I/345

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

Design and development of NIRSPEC: a near-infrared echelle spectrograph for the Keck II telescope

We thank James Kirk, Jon Otegi, B.J. Fulton, and Leonardo dos Santos for sharing data from their recent papers. We thank Greg Doppmann for his assistance with the observations, and we thank Elisabeth Newton for discussions about the radial velocity of GJ 1214. We thank Allison Youngblood for helpful conversations regarding the stellar EUV reconstruction and Ilaria Carleo for sharing the results on the observed Ly α flux of GJ 9827. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Software: astropy (Astropy Collaboration et al.