Helium in the Extended Atmosphere of the Warm Superpuff TOI-1420b

Superpuffs are planets with exceptionally low densities (ρ ≲ 0.1 g cm−3) and core masses (M c ≲ 5M ⊕). Many lower-mass (M p ≲ 10M ⊕) superpuffs are expected to be unstable to catastrophic mass loss via photoevaporation and/or boil-off, whereas the larger gravitational potentials of higher-mass (M p ≳ 10M ⊕) superpuffs should make them more stable to these processes. We test this expectation by studying atmospheric loss in the warm, higher-mass superpuff TOI-1420b (M = 25.1M ⊕, R = 11.9R ⊕, ρ = 0.08 g cm−3, T eq = 960 K). We observed one full transit and one partial transit of this planet using the metastable helium filter on Palomar/WIRC and found that the helium transits were 0.671% ± 0.079% (8.5σ) deeper than the TESS transits, indicating an outflowing atmosphere. We modeled the excess helium absorption using a self-consistent 1D hydrodynamics code to constrain the thermal structure of the outflow given different assumptions for the stellar XUV spectrum. These calculations then informed a 3D simulation, which provided a good match to the observations with a modest planetary mass-loss rate of 1010.82 g s−1 ( Mp/Ṁ≈70 Gyr). Superpuffs with M p ≳ 10M ⊕, like TOI-1420b and WASP-107b, appear perfectly capable of retaining atmospheres over long timescales; therefore, these planets may have formed with the unusually large envelope mass fractions they appear to possess today. Alternatively, tidal circularization could have plausibly heated and inflated these planets, which would bring their envelope mass fractions into better agreement with expectations from core-nucleated accretion.


INTRODUCTION
Thirty years of exoplanet discovery have revealed planetary archetypes for which we have no direct counterparts in our solar system.Among these myriad surprises are the super-puffs, a class of low-mass planets with densities ρ ≲ 0.1 g cm −3 (Lee & Chiang 2016;Jontof-Hutter 2019).For planets with M p ≲ 50M ⊕ , such low densities imply core masses smaller than 5M ⊕ , which are challenging to reconcile with the core accretion theory for planet formation unless the envelopes are Corresponding author: Shreyas Vissapragada shreyas.vissapragada@cfa.harvard.eduaccreted "dust-free" (Lee & Chiang 2016) or atop icerich cores (Stevenson 1984;Venturini et al. 2015) beyond the water ice line.Moreover, if super-puffs are indeed as gas-rich as their densities imply, then many should be prone to catastrophic mass loss via photoevaporation and/or boil-off, with atmospheric lifetimes shorter than or comparable to the system ages (i.e.M p / Ṁ ≲ 1 Gyr for catastrophic loss; Cubillos et al. 2017;Wang & Dai 2019;Chachan et al. 2020).A number of mechanisms have been proposed to help ease these tensions and bring super-puff envelope mass fractions into closer agreement with expectations from core accretion, including inflation by tides (Millholland 2019;Millholland et al. 2020) or heat from formation (Masuda 2014;Libby-Roberts et al. 2020), planetary rings (Akinsanmi et al. 2020;Piro & Vissapragada 2020;Ohno & Fortney 2022), and high-altitude aerosols (Kawashima et al. 2019;Wang & Dai 2019;Gao & Zhang 2020;Ohno & Tanaka 2021).
We show the super-puffs in context on a M p − T eq diagram in Figure 1.All transiting planets on the NASA Exoplanet Archive (Akeson et al. 2013) with densities ρ < 0.3 g cm −3 and masses measured to better than 30% are shown, and planets on grazing orbits have been removed.Most puffy planets hotter than 1000 K are typical gas giants, with large H 2 /He envelopes that are readily explained by core accretion models (e.g. Lee 2019).Gas giants above this temperature threshold have long been known to have inflated radii driven by the large incident stellar flux (e.g.Bodenheimer et al. 2003;Thorngren & Fortney 2018).Planets usually considered "super-puffs" (i.e.those with inferred core masses too small to trigger classical runaway gas accretion) are found in the lower left of this figure.The super-puff and gas giant regimes appear to be separated by a deficit of planets, which may result from sensitivity biases as a function of orbital period: most gas giants have masses measured with radial velocities, and most super-puffs with transit timing variations (Steffen 2016;Mills & Mazeh 2017).
The population is shown in more detail in the right panel of Figure 1.We focus on the planets on the upper right corner of the plot, WASP-107b (Anderson et al. 2017) and TOI-1420b (Yoshida et al. 2023).Both have inferred core masses < 5M ⊕ (Piaulet et al. 2021;Yoshida et al. 2023), too small to trigger runaway gas accretion in the traditional core accretion paradigm for gas giants (which requires core masses ≳ 10M ⊕ ; e.g.Bodenheimer & Pollack 1986).However, these planets stand out as warmer and higher-mass than the other superpuffs.They also differ from the rest of the super-puff population in system architecture: WASP-107b is accompanied only by a distant non-resonant companion (Piaulet et al. 2021), and TOI-1420b also lacks nearby companions (Yoshida et al. 2023).All other super-puffs are found in multi-planet systems near mean-motion resonances, which facilitates their characterization via transit timing variations (Jontof-Hutter 2019), although a few have now also been detected with radial velocities (Petigura et al. 2018;Santerne et al. 2019).
Although their short orbital periods (P < 10 days) make them more prone to XUV-driven photoevaporation, TOI-1420b and WASP-107b are not expected to be susceptible to catastrophic atmospheric loss, as their gravitational potentials are relatively large despite their low densities (Wang & Dai 2019;Chachan et al. 2020;Belkovski et al. 2022).Indeed, He 1083 nm measurements of WASP-107b suggest an outflow rate of about one Earth mass per gigayear (M/ Ṁ ≈ 30 Gyr; Spake et al. 2018;Allart et al. 2019;Kirk et al. 2020;Spake et al. 2021;Wang & Dai 2021), in line with predictions from photoevaporation models (Kubyshkina et al. 2018;Caldiroli et al. 2022).This planet has not lost much envelope mass over its lifetime, so it remains plausible that WASP-107b formed with an unusually large envelope (Lee & Chiang 2016;Piaulet et al. 2021).This is a major difference from the lower-mass super-puffs, which would be expected to lose large envelopes on short timescales necessitating many of the aforementioned alternative mechanisms (rings, hazes, etc.).
In this work, we study atmospheric escape in TOI-1420b (Yoshida et al. 2023).This planet orbits a late G dwarf every 6.96 days, and with M p = 25.1±3.8M⊕ and R p = 11.9 ± 0.3R ⊕ , TOI-1420b has an exceptionally low density of just 0.08±0.02g cm −3 .Moreover, at J = 10.6 the system is well-positioned for narrowband photometric observations with the Wide-field InfraRed Camera (WIRC; Wilson et al. 2003) on the Hale 200-inch Telescope at Palomar Observatory.We recently used Palomar/WIRC to constrain helium absorption in three systems with similar J magnitudes: HAT-P-18 (J = 10.8),WASP-52 (J = 10.6), and HAT-P-26 (J = 10.1;Vissapragada et al. 2020a;Paragas et al. 2021;Vissapragada et al. 2022b).The former two measurements have recently been confirmed by JWST NIRISS/SOSS and Keck/NIRPSEC, demonstrating the utility of narrowband photometry for this science case (Fu et al. 2022;Kirk et al. 2022).In Section 2, we describe our narrowband photometric observations of TOI-1420b.We model our narrowband helium light curves in Section 3, interpret the results using outflow models in Section 4, discuss the implications for warm, high-mass super-puffs in Section 5, and conclude in Section 6.

HELIUM LIGHT CURVES
We obtained transit photometry of TOI-1420b with Palomar/WIRC.All events were observed with the beam-shaping diffuser (Stefansson et al. 2017;Vissapragada et al. 2020b), which spreads the light into a tophat point-spread function with a 3 ′′ full-width at halfmaximum.This allows us to control time-correlated systematics and to integrate for longer on bright sources (improving the duty cycle of the observations).We used the custom metastable helium filter (centered on 1083.3 nm with a bandwidth of 0.64 nm) presented in Vissapragada et al. (2020a).
We observed a full transit of TOI-1420b on UT 2022 June 30 from 05:17 (airmass 1.8) to 12:02 (airmass 1.2) and a partial transit on UT 2022 July 21 from 05:28 (airmass 1.5) to 11:08 (airmass 1.2).The transit duration is T 14 = 3.37±0.15hr (Yoshida et al. 2023).We attempted a partial transit observation on UT 2022 July 28, but poor weather rendered the data unusable for our analysis.All observations were taken with 90 s exposures.To correct the structured background imposed by the narrowband helium filter (further details in Vissapragada et al. 2020a), we obtained 8 dither frames immediately before both transit observations.Additionally, we maintained roughly the same detector positioning for both nights, so we observed the same six reference stars each night (which were selected to be the brightest sources in the field besides TOI-1420, ranging from J = 11.3 to J = 8.8).None of the selected reference stars are known variables, and none exhibited detectable shortterm variability in our photometry.
We dark-subtracted, flat-fielded, and backgroundcorrected the data using the methods presented previously in Vissapragada et al. (2020aVissapragada et al. ( ,b, 2022b) ) and Greklek-McKeon et al. (2023).We then performed aperture photometry on the target and reference stars using photutils (Bradley et al. 2023).We tested apertures from 5 pixels (1.′′ 25) to 20 pixels (5.′′ 00) in radius.The centers of the apertures were allowed to shift to track variations in telescope pointing.Our guiding was stable on both nights, with centroid variations not exceeding 1.5 pixels.We estimated and corrected any residual local background using annular apertures around each source with inner radii of 25 pixels (6.′′ 25) and outer radii of 50 pixels (12.′′ 50).The optimal aperture was selected to be the one that minimized the per-point rms scatter in the target star photometry with the average comparison star divided out (removing common-mode systematics) and a moving-median filter applied (removing the transit event).We adopted optimal apertures of 8 pixels (2.′′ 00) for the first night and 7 pixels (1.′′ 75) for the second night.

LIGHT CURVE MODELING
We first fit each helium curve separately using the exoplanet code (Foreman-Mackey et al. 2021a,b).We fit the target light curve f as the product of a transit light curve model T and a systematics model S: (1) To model the transit, we included a uniform prior on the planet-to-star radius ratio R p /R ⋆ , scaled semimajor axis a/R ⋆ , and impact parameter b, uninformative priors (Kipping 2013) on the quadratic limb darkening coefficients (u 1 , u 2 ), and normal priors on the planetary orbital period P , transit epoch t 0 , and stellar radius R ⋆ .Normal priors were taken from Yoshida et al. (2023), and   are reproduced in Table 1.The planet does not appear to be on an eccentric orbit (e < 0.17 Yoshida et al. 2023) so we fixed the eccentricity to 0. We tracked the depth of the WIRC light curve at mid-transit δ WIRC numerically rather than using (R p /R ⋆ ) 2 so that limb darkening was taken into account.
We included the two minute cadence TESS light curve of TOI-1420b from Yoshida et al. (2023) in all our fits to provide a reference R p /R ⋆ for the WIRC light curves.The TESS data also helped us to constrain the degeneracy between b and a/R ⋆ with the precise transit shape.All transit parameters were shared between the WIRC and TESS fits except for R p /R ⋆ and the quadratic limb darkening coefficients (u 1 , u 2 ).When we tracked the comparison transit depth δ comp , we recalculated the TESS light curve with the WIRC limb darkening coefficients.This excluded the possibility that perceived depth differences were caused by differences in limb darkening between the two bandpasses.The excess depth at mid-transit, which describes the extra absorption due to metastable helium, was then calculated as We modeled the systematics S as a linear combination of detrending vectors s i including the six comparison star light curves and the mean-subtracted times (defining a linear trend in time): The weight w i for each detrending vector s i in the linear combination was a free parameter with uniform prior U(−2, 2).We also tested including the airmass, the centroid offsets, and the water absorption proxy described previously in Paragas et al. (2021) and Vissapragada et al. (2022b) as additional detrending vectors, and calculated the Bayesian Information Criterion (BIC; Schwarz 1978) to determine whether the increase in log-likelihood justified the addition of these parameters to the model.However, none of these additional detrending vectors improved the BIC, i.e. the marginal increase in final log-likelihood did not justify the addition of any of the parameters, so we left these out of the final models.Finally, we included a photometric jitter term describing the excess white noise in the WIRC data (added in quadrature to the photon noise for each point) and a multiplicative scaling factor for the uncertainties on the TESS data.In all of the fits, we first optimized towards the maximum a posteriori (MAP) solution with pymc3 (Salvatier et al. 2016).We identified and removed 5σ outliers from the residuals to the MAP fit, iterating the MAP fit and outlier rejection procedure until a stable MAP solution was obtained.Then, we used the No U-Turn Sampler (Hoffman & Gelman 2011) in pymc3 to sample the posterior distributions for the model parameters.
Every fit was run with 10,000 tuning steps before taking 10,000 posterior draws (distributed across two chains).To assess convergence, we verified that the Gelman-Rubin statistic (Gelman & Rubin 1992) R ≪ 1.01 for all sampled parameters, and we also visually inspected the chains to ensure they were well-mixed.The posteriors of our single-planet fits are given in Table 1.
The transit depths on the two nights were consistent at about the 1σ level.The WIRC transit depth was less precisely constrained on the second night than the first because of the partial phase coverage, but the second night also offered more post-egress baseline than the first night.Given the reasonable agreement in R p /R ⋆ between both datasets and the complimentary baseline coverage, we decided to fit both WIRC light curves jointly with the TESS light curve.The posteriors of the joint fit are included in Table 1. Figure 2 shows the joint WIRC light curve alongside the phase-folded TESS data.The final WIRC transit was 0.671±0.079%deeper than the TESS transit, indicating substantial extra absorption by metastable helium in the upper atmosphere of TOI-1420b at a confidence of 8.5σ.
The rms scatter of the unbinned data was 0.50% on the first night and 0.63% on the second night.These were 1.2× and 1.3× the photon noise on each night, respectively, indicating good photometric performance.This is the closest to photon noise we have been able to achieve with Palomar/WIRC using the helium filter (Vissapragada et al. 2022b), likely because six good comparison stars were available in the field.However, when reducing the second night of data we noticed a large correlated noise component at about a 30 min timescale in the Allan deviation plots (Figure 3).A correlated feature with this timescale is clearly visible at the end of the WIRC light curve in Figure 2, though the uncertainties in this region of the light curve suggest it is not very significant.This feature was not clearly associated with rapid changes in water absorption proxy, centroid offset, or total flux.Additionally, the second night did not contribute much information to the final excess depth measurement (δ mid = 0.608 +0.085 −0.085 % for the first night only and δ mid = 0.671 +0.077 −0.079 % for the joint fit), so we proceeded without treating this feature in detail.1D modeling, we calculated a 3D model which treats the thermodynamics approximately, but captures the interaction of the outflow with the tidal and stellar wind environment.

Characteristic temperatures
TOI-1420b's low density carries some important implications for the modeling of its outflow.
The one-dimensional, transonic Parker wind models often adopted for modeling metastable helium absorption (e.g.Oklopčić & Hirata 2018;Lampón et al. 2020; Dos Santos 2022) begin with subsonic flow near the planet that accelerates to supersonic flow at larger radii.A transonic wind requires R s > R p , where: (the tidal gravity correction is small for TOI-1420b; see Appendix E of Vissapragada et al. 2022b).The sound speed of the isothermal outflow is: where k B is Boltzmann's constant, T 0 is the isothermal sound speed, μ is the mean molecular weight and m H is the mass of a hydrogen atom.The condition R s > R p can equivalently be written Φ > 2c 2 s , where Φ = GM p /R p is the gravitational potential.
TOI-1420b has a small gravitational potential of Φ ≈ 10 12 erg g −1 , so only isothermal outflows with c s ≲ 8 km s −1 yield transonic velocity structures (i.e.R s > R p ).In terms of an isothermal temperature, this corresponds to T 0 ≲ 7700 K for an assumed mean molecular weight of µ = 1 amu.This condition overlaps significantly with the range of isothermal outflow temperatures usually assumed when modeling outflows from hot giant planets (5000K ≲ T 0 ≲ 15000K; e.g.Linssen et al. 2022;Vissapragada et al. 2022a).To avoid the issue of applying the isothermal Parker wind outside the bounds of its validity, we decided to consider one-dimensional models with self-consistent (non-isothermal) thermodynamics and the nature of the transonic ouflows that develop.

One-dimensional outflow models
To self-consistently treat the wind-launching physics, we used the 1D photoionization hydrodynamics code ATES to compute outflow density, velocity, and temperature profiles (Caldiroli et al. 2021).Contrary to an isothermal Parker wind profile, ATES calculates the temperature structure based on the stellar irradiation assuming a pure H/He composition for the wind.With the equilibrium temperature as a boundary condition at R p , the wind starts hydrostatically and is indeed transonic.We then passed the outflow structure to the NLTE photoionization code Cloudy (Ferland et al. 1998(Ferland et al. , 2017) ) to obtain the metastable helium density profile.We finally took Cloudy's metastable helium density profiles and used a custom radiative transfer module (based on Oklopčić & Hirata 2018; Linssen et al. 2022) to calculate the mid-transit spectrum in the stellar rest frame, using the limb-darkening coefficients and transit impact parameter of Table 1.We assumed a hydrogen/helium abundance ratio of 90/10 by number.Finally, we Doppler The XUV flux of the host star is an important ingredient in the simulations, as it governs both the thermal structure in ATES and the metastable helium population in Cloudy.However, the high-energy spectrum of TOI-1420 is unknown (and challenging to measure at d = 200 pc), so we used spectra of similarly active stars as proxies.TOI-1420 is a late G dwarf (T eff = 5493 K) with log(R ′ HK ) = −4.87(Yoshida et al. 2023), indicating similar chromospheric activity to the mean solar value of −4.91 (Mamajek & Hillenbrand 2008).To cover the uncertainty in the stellar XUV spectrum, we repeated the modeling using three different proxies: the MUSCLES spectrum of HD 40307 (an early K star representing "high activity" albeit with a modest log(R ′ HK ) = −4.99;Mayor et al. 2009;France et al. 2016;Loyd et al. 2016;Youngblood et al. 2016), a quiet solar spectrum (representing "medium activity"; Oklopčić & Hirata 2018; Dos Santos 2022), and the solar spectrum with EUV (< 1000 Å) scaled down by a factor of three (representing "low activity").
The high and medium activity models predict massloss rates of 10 11.25 g s −1 and 10 10.93 g s −1 , respectively, but both predict much larger helium signals than we observed (Figure 4).The low activity spectrum gives the best match to the data at a mass-loss rate of 10 10.58 g s −1 , and we therefore adopt this as our preferred XUV model for the star.At the same time, the high and medium activity solutions should not be thought of as "ruled out" on the basis of the 1D modeling alone.Some of the absorption in these models comes from high altitudes > 4R p where shadowing, tidal effects, and interaction with the stellar wind becomes important for predicting the metastable helium signature (Wang & Dai 2021;MacLeod & Oklopčić 2022).These effects are not treated in ATES, so rather than reducing the activity in the 1D models further, we turned to 3D modeling informed by the results of the self-consistent 1D analysis.

Three-dimensional outflow model
In order to assess the impact of the stellar environment on the outflow of TOI-1420b -including the stellar gravity and stellar wind -we perform a three-dimensional hydrodynamic calculation of the escape of a planetary outflow similar to that of the TOI-1420 system.Our model is based on the Athena++ Eulerian hydrodynamic code (Stone et al. 2008(Stone et al. , 2020) ) and the planetary outflow modeling and radiative transfer post processing software developed by MacLeod & Oklopčić (2022).This simulation does not include magnetic fields, which can suppress the planetary helium signature (Schreyer et al. 2023).
We simulate the interaction of planetary and stellar wind outflows in the corotating reference frame.We assume that the planet and star rotate with the orbital frequency, and adopt the basic system parameters of Table 1.In our model the planetary and stellar winds are simulated with a nearly-isothermal adiabatic index of γ = 1.0001, but have different characteristic temperatures.The planetary outflow is launched from a uniform boundary condition representing the planet, and is defined by the hydrodynamic escape parameter λ p = GM p /R p c 2 s , where c s is the sound speed of the isothermal planetary outflow.We adopt λ p = 5, which corresponds to T ≈ 3200 K for µ = 1 amu.The value of the hydrodynamic escape parameter is chosen to roughly correspond to the 1D ATES model with the 1/3× Solar EUV spectral energy distribution in Figure 4. Similarly, on the basis of the ATES modeling, we adopt the 1/3× Solar EUV as our spectrum for the threedimensional model spectra.
With these basic parameters of the system set by either the known system properties (e.g.Table 1), or by comparison to the 1D ATES models, we are free to vary the relative strength of the planetary and stellar outflows in our hydrodynamic calculation by specifying the stellar and planetary boundary conditions that generate these flows.This freedom allows us to scale our calculation to reproduce the mid-transit excess absorption.We chose the density of our planetary boundary condition such that the measured emergent outflow rate in steady state is 6.7 × 10 10 g s −1 (10 10.82 g s −1 ).The stellar boundary condition is set similarly, but with a hydrodynamic escape parameter of λ * = 10, and a mass-loss rate of 1.6 × 10 11 g s −1 , approximately 3.2× the planetary mass loss rate, or ≈ 3.5 × 10 −15 M ⊙ yr −1 .We note that this last choice is somewhat arbitrary, as only much (e.g.order-of-magnitude) stronger or weaker stellar winds are ruled out.Finally, we evolve the system for roughly 6 orbital periods (41.7 d), and find that the flow reaches a steady-state after ∼ 10 d.
The upper panel of Figure 5 shows the emergent flow in a slice through the orbital plane.The star lies at the origin, and the planet is located in the −x-direction.An inset shows the flow in a region of ±20R p (roughly ±0.1 au).As the planetary outflow expands into the stellar environment, it is shaped by the tidal gravity and stellar wind.Near the planet, outflowing material is initially roughly homogeneous.Outside the planetary Hill sphere of radius ≈ 4.1R p , however, the stellar gravity dominates and the nearly-radial outflow from the planet is truncated by shocks and sheared into tidal tails leading and trailing the planet.Interaction with the stellar wind also partially shapes these tails and their kinemat-ics in our model, though we note that the value of the stellar wind mass loss rate is not measured in TOI-1420, so the influence of the stellar wind could be larger or smaller than our nominal assumption (e.g. as described by MacLeod & Oklopčić 2022).
We post-process metastable helium level populations and transit spectra in the snapshot shown in Figure 5.For the stellar SED, we adopt the low-activity model described in the preceding section, with 1/3 of the solar EUV flux.We adopt an observer in the general −xdirection but at different viewing angles in order to represent the passage of time as the planet transits.The mid-transit model spectrum is similar to the low-activity model of the preceding section.We then compute the excess absorption in the WIRC filter passband following Equation (3) of Vissapragada et al. (2020a).
The resulting excess absorption light curve is shown in the lower panel of Figure 5.This model achieves similar excess depth at mid-transit to that observed in Figure 2, and shows that the vast majority of absorption is concentrated during the optical transit duration.The extended tails of planetary atmosphere material are effectively hidden by their low optical depth.Thus the region that forms the metastable helium absorption lies within a few planetary radii -near the peak of metastable helium density in the center panel of Figure 4.At these scales the planetary outflow is nearly spherical and our models thus predict spectra similar to those of the 1D models despite the complexity of the largescale flow.This optical depth effect can be contrasted with outflows like those observed for HAT-P-32b and HAT-P-67b, where the larger mass-loss rates give the tails greater optical depths, making them readily observable even when compared to the mid-transit absorption (Gully-Santiago et al. 2023;Zhang et al. 2023).
Atmospheric retention has been a major puzzle for super-puff planets.Assuming their measured radii correspond to typical photospheric pressures on the order of 10 mbar, many lower-mass super-puffs with M p ≲ 10M ⊕ (including e.g. ) are expected to be unstable to catastrophic atmospheric escape via photoevaporation or boil-off.This has led to suggestions that the measured R p must correspond to lower pressures in these planets, which could be a result of high-altitude clouds, dust, hazes, or rings (Lammer et al. 2016;Cubillos et al. 2017;Wang & Dai 2019;Chachan et al. 2020 Figure 5. Large-scale model of the planetary outflow from TOI-1420b.The upper panel shows the logarithm of the mass density log 10 (ρ) of a slice through the orbital plane of the three-dimensional hydrodynamic model.A small bubble of unshocked mass loss around the planet is truncated by the star's tidal gravitational field, after which the outflow extends into leading and trailing tails ahead and behind the planet in its orbital path, respectively.The lower panel shows the corresponding excess absorption in the WIRC helium bandpass as a function of time during a transit, and shaded regions denote the 1σ (dark gray) and 2σ (light gray) uncertainties on the measured mid-transit excess absorption.Despite the extended tails, most of the absorption comes during the optical transit (approximately ±0.07 d, indicated by vertical lines) with a gradual return to very low levels of excess absorption outside of ±0.1 d.This indicates that most of the optical depth contributing to the excess absorption lies in dense regions at the base of the planetary outflow near the planet -not from the large-scale structure, which is mostly transparent.
Our helium observations of TOI-1420b agree well with the theoretical prediction that high-mass super-puffs should be stable against atmospheric loss.The threedimensional outflow model in Section 4 matched the observed helium signature at a mass-loss rate of Ṁ = 10 10.82 g s −1 , corresponding to a mass-loss timescale M/ Ṁ ≈ 100 Gyr.The situation is similar for this planet's doppelgänger WASP-107b, for which Wang & Dai (2021) similarly found a good fit to helium observations with M/ Ṁ ≈ 30 Gyr; moreover their model provided a remarkably good match to the observed posttransit tail for WASP-107b (Spake et al. 2021).The inferred outflow rates for both these planets are not expected to result in significant mass loss over the planetary lifetime.Therefore, despite their remarkably low densities, super-puffs with M p ≳ 10M ⊕ indeed appear stable against atmospheric loss.

Is TOI-1420b tidally-inflated?
While atmospheric retention may not be a problem for high-mass super-puffs, it remains challenging to form planets with such large envelope mass fractions (≳ 80% Piaulet et al. 2021;Yoshida et al. 2023) within an in situ core accretion framework (Lee 2019).One possible resolution is that these planets accreted their envelopes far from their stars (≳ 1 au), where dust-free atmospheres could cool and accrete rapidly (Lee & Chiang 2016;Piaulet et al. 2021).Migration would then be required to bring the planets to their observed positions, and indeed many super-puffs are observed to be parts of resonant chains that suggest such a migration history (Jontof-Hutter 2019).However, TOI-1420b and its doppelgänger WASP-107b do not reside in resonant configurations.WASP-107b is accompanied by a distant companion that might have driven past migration onto the current high-obliquity orbit (Dai & Winn 2017;Piaulet et al. 2021;Rubenzahl et al. 2021), but no additional companions are yet known in the TOI-1420 system (Yoshida et al. 2023).
We therefore turn to tidal inflation as a possible alternative explanation for the low apparent density of TOI-1420b.Although tidal inflation is not required to match the outflow data, the ATES+Cloudy approach used in Section 4 is agnostic to this mechanism.The underlying envelope mass fraction is immaterial to the model so long as the density and temperature at the lower boundary (R p ) are correctly described.Moreover, the overall escape rate of the flow is relatively insensitive to the boundary density (Salz et al. 2016;Caldiroli et al. 2022).For WASP-107b, the radial velocity data permit a small eccentricity (e < 0.14 at 2σ) and the resulting heating from eccentricity tides could bring the true en-velope mass fraction below 50% (Millholland et al. 2020;Piaulet et al. 2021), though this conclusion depends sensitively on the unknown planetary obliquity and tidal quality factor.The radial velocity curve for TOI-1420b permits a similarly small eccentricity (e < 0.17 at 2σ; Yoshida et al. 2023), so the same mechanism may apply.
Figure 6.The radius, eccentricity, and semimajor axis evolution for several possible cases of TOI-1420b's tidal evolution.The colored points indicate parts of each evolution track that match the planet's present-day radius (1.061±0.029RJ),eccentricity (< 0.17), and semi-major axis (0.0710 ± 0.0012 au).The blue curves indicate a planet that began on a circular orbit, where the large radius is due entirely to a low bulk metallicity (in this example, Zp = 0.13).In red is the "young circularizing" case, in which the planet is relatively young (∼ 1 Gyr) and at the final stages of fairly rapid circularization (log(Qp) = 5.5), allowing for a hotter planet with larger bulk metallicity Zp = 0.37.Finally, the yellow curves indicate a slower circularization track (log(Qp) = 6, Zp = 0.33), which is a potential compromise between the young age of the red line and the very low metallicity blue line.These are only three illustrative examples of a larger solution space (see Section 5.2).
To investigate the possibility that tidal circularization has heated and inflated the planet, we constructed a model which combines the giant planet thermal evolution models of Thorngren & Fortney (2018) with the equilibrium tides model (e.g.Goldreich & Soter 1966;Mardling & Lin 2002;Jackson et al. 2008).Similarly to Lopez & Fortney (2016), we are not including tidal inspiral or obliquity tides in this model, but we do calculate the amount of heating resulting from circularization to add to the energy balance of the interior.Thus for a constant planetary M p , bulk metallicity Z p , and tidal quality factor Q p , we evolve the eccentricity, semimajor axis, and interior specific entropy over time.
There is a broad space of possible solutions for TOI-1420b, as the model has a number of free parameters (e 0 , a 0 , Q p , Z p and age) but only strong constraints on a and R p .Therefore, we provide three illustrative evolution tracks for TOI-1420b in Figure 6.For reference we show the case in which the planet formed on a circular orbit, which requires an unusually low metallicity of Z p = 0.13 for a planet of this mass to match the radius (for more details on this nominal model, see Section 4 of Yoshida et al. 2023).This model corresponds to a maximum core mass (assuming all metals are confined to a rocky core) of Z p M p = 3.3M ⊕ .If instead the planet started with an eccentric orbit, the semi-major axes and eccentricity follow a sigmoid curve with the logarithm of the age.During this process of circularization a substantial portion of the orbital energy is injected into the planet, resulting in an increase in the planet's radius which decays as the excess heat is radiated away.The red and yellow curves of Figure 6 show how we may be observing the planet towards the tail end of this process, after it has lost most of its initial eccentricity (consistent with the upper limit from Yoshida et al. 2023) but while it still has a moderately hotter interior than it would otherwise.This extra heat allows us to reproduce the observed radius with higher bulk metallicities Z p = 0.3 − 0.4, implying core masses that are more in line with core-nucleated accretion (similarly to the case for a tidally-heated WASP-107b; Millholland et al. 2020;Piaulet et al. 2021).Different values for the tidal Q p are correlated with present-day ages in the solution space: having an age over a gigayear (which seems likely for this system; Yoshida et al. 2023) would indicate that the planetary Q p is probably not much less than 10 6 , lest the planet have long since finished circularizing.
Our calculations demonstrate that tidal inflation is a plausible way to explain the planet's large radius.However, they also demonstrate the strong degeneracies inherent to the problem: all available data for TOI-1420b can be matched equally well in models where the planet is young and circularizing or old and very gasrich.With relatively uninformative prior information on TOI-1420b's eccentricity (< 0.17) and age (< 10.7 Gyr, Yoshida et al. 2023), these degeneracies cannot be broken.Future efforts to more precisely constrain the eccentricity and/or age of the planet would therefore be highly valuable for understanding TOI-1420b's evolutionary state.5.3.Can high-altitude hazes "deflate" TOI-1420b?
The flat transmission spectra of lower-mass superpuffs are well-modeled by high-altitude aerosols, which may also be responsible for their large apparent radii (Kawashima et al. 2019;Gao & Zhang 2020;Ohno & Tanaka 2021).It is worth considering this hypothesis in more detail for the higher-mass super-puff TOI-1420b, as the model requires strong outflows (similar to the one we have detected) to carry aerosols to high altitudes.While high-altitude aerosols are ruled out as a possibility for WASP-107b (molecular features are detected lower in its atmosphere; Spake et al. 2018;Kreidberg et al. 2018;Dyrek et al. 2023), transmission spectroscopy sensitive to the lower atmosphere of TOI-1420b is not yet available.
Hazes can persist at µbar to nbar pressures in atmospheres undergoing outflows (Wang & Dai 2019;Gao & Zhang 2020;Ohno & Tanaka 2021), and so the 1 nbar radius should represent the maximum extent to which hazes can increase a planet's radius from a clear atmosphere case.As such, to test whether TOI-1420b could be a smaller planet with optically thick high-altitude hazes, we compute the 1 nbar radius of this planet using the atmospheric and planetary structure model of Gao & Zhang (2020).Briefly, the model assumes a core with Earth-like composition and a mass-radius relation adopted from Zeng et al. (2019), which lies beneath a gas envelope of solar composition and a simplified thermal structure consisting of a deep convective zone and an upper radiative zone that is connected to the convective zone at the radiative-convective boundary (RCB).The convective zone structure is solved following Owen & Wu (2017) assuming hydrostatic equilibrium and a polytrope equation of state with a power of 7/5 to represent a largely H 2 composition.At the RCB, the temperaturepressure profile transitions smoothly to a radiative profile following Piso & Youdin (2014), with the temperature of the RCB set to the equilibrium temperature of TOI-1420b.For the opacity in the radiative zone we use the tables from Freedman et al. (2014).We refer the reader to Gao & Zhang (2020) for more details regarding the model.Predicted planet radii at a pressure level of 1 nbar in the atmosphere for a range of atmospheric mass fractions and intrinsic temperatures Tint of 25 K (dark red), 50 K (red), 75 K (orange), and 100 K (yellow).The span in the shaded regions represents the impact of the uncertainty in the planet mass on the model.The 1σ range of the observed radius is shown in the gray region.For atmospheric mass fractions <0.3, no model is able to match the observed radius for a haze that is optically thick in the slant geometry at 1 nbar.
In Figure 7, we show that the 1 nbar radius of TOI-1420b does not reach the observed radius for any atmospheric mass fractions <30% and intrinsic temperatures <100 K, even when the uncertainty in the mass is taken into account.These limits roughly correspond to those of Neptune-mass planets with the age and insolation of TOI-1420b (Helled et al. 2011;Lopez & Fortney 2014).Therefore, it is unlikely TOI-1420b is a smaller planet that has been made to look larger with a high-altitude haze layer.

CONCLUSION
Many lower-mass super-puffs (M p < 10M ⊕ ) are expected to have atmospheric lifetimes shorter than their ages, while higher-mass (M p > 10M ⊕ ) super-puffs are expected to be stable to catastrophic atmospheric loss even when strongly irradiated.To test this expectation, we searched for an outflow in the atmosphere of TOI-1420b, a warm (T eq = 960 K), higher-mass (M p = 25.1M⊕ ) super-puff with a density of just ρ = 0.08 g cm −3 (Yoshida et al. 2023).We observed one full transit and one partial transit with the narrowband helium filter on Palomar/WIRC and compared the resulting light curve to the TESS light curve.We found that the transit in the helium bandpass was consistently 0.671±0.079%deeper than the TESS transit, indicating excess absorption by outflowing metastable helium at a confidence of 8.5σ.
We first interpreted our observations using the 1D photoionization hydrodynamics code ATES coupled to Cloudy, which allowed for a self-consistent treatment of the outflow thermodynamics.The model that best matched our data used a low activity proxy for the host star's XUV spectrum (1/3× solar XUV) with a hydrodynamic escape parameter λ p ≈ 5. Using this XUV spectrum and escape parameter, we then calculated a threedimensional isothermal outflow model using Athena++.This model gave a good match to the observed excess absorption at mid-transit with a mass-loss rate of Ṁ = 4.9 × 10 10 g s −1 .The corresponding mass-loss timescale is M p / Ṁ ≈ 70 Gyr; despite its remarkably low density, the mass-loss rate of TOI-1420b is quite modest.The situation is similar for the other warm high-mass super-puff WASP-107b, which has M p / Ṁ ≈ 30 Gyr (Wang & Dai 2021).Our observations confirm that super-puffs with M p > 10M ⊕ are easily able to retain their massive atmospheres over long timescales.These planets may therefore have formed with unusually massive envelopes, likely accreting gas outside the water ice line before migrating to their present positions (Stevenson 1984;Venturini et al. 2015;Lee & Chiang 2016;Piaulet et al. 2021).
We also explored a number of alternative hypotheses for TOI-1420b's large radius.We presented some illustrative evolution calculations showing that all available data for TOI-1420b could be matched at a larger bulk metallicity (more consistent with models of core accretion in the inner disk; Lee 2019) provided the planet recently underwent tidal circularization.On the other hand, we found that high-altitude hazes are disfavored as an explanation for TOI-1420b's large radius.The planet is simply too massive, resulting in a relatively small scale height among the super-puffs: even opticallythick hazes at nbar pressure for a Neptune-like planet would fall well short of explaining TOI-1420b's large radius.This bodes well for future atmospheric characterization efforts: we do not anticipate that atmospheric features will be completely obscured by hazes, similar to the case for WASP-107b (Kreidberg et al. 2018;Spake et al. 2018;Dyrek et al. 2023).
TOI-1420b's helium absorption is the largest signal we have observed at high confidence with Palomar/WIRC thus far (Vissapragada et al. 2022b).We encourage other investigators to follow up the helium absorption for this planet with high-resolution spectroscopy, which would provide additional constraints on the atmospheric modeling.With such a strong signal this should also be a good candidate to search for a Doppler-shifted line core and/or extended egress, which is indicative of an outflow morphology sculpted by stellar winds in WASP-69b and WASP-107b (Nortmann et al. 2018;Spake et al. 2021;Wang & Dai 2021;MacLeod & Oklopčić 2022;Tyler et al. 2024).Our light curve of TOI-1420b appears symmetric, and is thus well-fit by a 3D model where the stellar wind is not strong enough to shepherd the outflowing material into a comet-like tail, but tail features can be challenging to capture at low effective resolving power from the ground (Vissapragada et al. 2020a;Paragas et al. 2021;Fu et al. 2022).High-resolution studies of TOI-1420b would therefore enable more direct morphological comparisons to planets like WASP-107b, and may even enable comparative studies of the stellar wind strengths in these systems.
Lower-mass super-puffs are expected to be even more susceptible to atmospheric escape than TOI-1420b and WASP-107b.In particular, there are a number of planets in Figure 1 that are similar in equilibrium temperature to WASP-107b, including TOI-216b, Kepler-79d, Kepler-18d, and Kepler-223e.If their host stars are active, these planets might also exhibit triplet helium absorption.While many super-puffs reside in systems that are too faint for current facilities, a number of suitable targets have been discovered in bright systems, including K2-24c, HIP 41378 f, and TOI-216b (Petigura et al. 2018;Santerne et al. 2019;McKee & Montet 2023).Although the transit durations for these planets tend to be long, helium absorption could still be constrained with multi-night ground-based campaigns (Gully-Santiago et al. 2023;Zhang et al. 2023) or spacebased observations with JWST (Fu et al. 2022;Dos Santos et al. 2023).By characterizing the upper atmospheres of lower-mass super-puffs, we can gain greater insight into the nature of these remarkably low-density worlds.
We thank Paul Nied and Diana Roderick for assistance with telescope operations.We additionally thank Fei Dai, Ryan Rubenzahl, and Jeremy Drake for helpful conversations.This work was supported by telescope time allocated to NASA-NSF Exoplanet Observational Research (NN-EXPLORE) partnership through the scientific partnership of the National Aeronautics and Space Administration, the National Science Foundation, and the NOIRLab.This research made use of photutils, an astropy package for detection and photometry of astronomical sources (Bradley et al. 2023).This research made use of exoplanet (Foreman-Mackey et al. 2021a,b) and its dependencies (Agol et al. 2020;Kumar et al. 2019;Astropy Collaboration et al. 2013, 2018;Kipping 2013;Luger et al. 2019;Salvatier et al. 2016;Theano Development Team 2016).laboration et al. 2013, 2018, 2022), p-winds (Dos Santos 2022), ATES (Caldiroli et al. 2021), Cloudy (Ferland et al. 1998(Ferland et al. , 2017)), Athena++ (Stone et al. 2008(Stone et al. , 2020) ) Table 2. Priors and posteriors for the TESS error scaling parameter, jitters, and detrending weights.
Figure 1.(Left) The masses and equilibrium temperatures of low-density (ρ < 0.3 g cm −3 ) planets on the NASA Exoplanet Archive as of 2023 September 10 (Akeson et al. 2013) with masses measured to better than 30%.Planets on grazing orbits have been removed.TOI-1420b is highlighted with the red star.The shading indicates planetary density, and the dotted lines demarcate the bounds of the plot shown on the right.(Right) Same as left, but focused on the lower-mass super-puff population with individual systems labeled.

Figure 2 .
Figure 2. WIRC (left) and phase-folded TESS (right) light curves of TOI-1420b.The data are shown unbinned in gray (with circles and triangles denoting points from the first and second nights of WIRC observations, respectively) and binned to 10 min cadence in black.Maximum a posteriori (MAP) light curve models are overplotted in the solid curves for WIRC (red) and TESS (blue), with the shaded regions denoting the 68% credible intervals.The TESS light curve model is recalculated using the helium bandpass limb darkening coefficients on the WIRC plot, showing that differences in limb darkening cannot explain the increased depth at mid-transit.

Figure 3 .
Figure 3. Allan deviation plots for the two WIRC light curves.The black curves indicate the rms scatter as a function of bin size on each night, the solid red curves indicate expectations from photon noise, and the dashed red curves are the solid red curves scaled up to match the first point of the black curve, i.e. the expectation if our data binned down exactly like white noise.

Figure 4 .
Figure 4. Temperature (top panel, solid curves), mean molecular weight (top panel, dashed curves), and metastable helium number density (middle panel) profiles for the high activity (blue), medium activity (red), and low activity (green) ATES+Cloudy models.The solid curves in the bottom panel show the predicted metastable helium features at high dispersion in air wavelengths in the observatory rest frame.The dotted curve shows the filter transmission function.The text annotations give the corresponding excess transit depth values (after integrating through the filter transmission function) compared to the observed absorption.
;Libby-Roberts et al. 2020;Piro & Vissapragada 2020; Figure7.Predicted planet radii at a pressure level of 1 nbar in the atmosphere for a range of atmospheric mass fractions and intrinsic temperatures Tint of 25 K (dark red), 50 K (red), 75 K (orange), and 100 K (yellow).The span in the shaded regions represents the impact of the uncertainty in the planet mass on the model.The 1σ range of the observed radius is shown in the gray region.For atmospheric mass fractions <0.3, no model is able to match the observed radius for a haze that is optically thick in the slant geometry at 1 nbar.

Table 1 .
Priors and posteriors for the light-curve fits.The WIRC depth δWIRC, comparison depth δcomp, and mid-transit excess depth δ mid = δWIRC − δcomp are all derived parameters.Posteriors for the detrending weights, WIRC jitter parameters, and TESS error scaling parameter are included in Table2.
Note-N (a, b) indicates a normal distribution with mean a and standard deviation b.U(a, b) indicates a uniform distribution between a and b.