A Multiplanet System's Sole Super-puff: Exploring Allowable Physical Parameters for the Cold Super-puff HIP 41378 f

Despite the large number of exoplanets that have been physically characterized, the physics of planet formation remains shrouded in mystery. In general, the only observable piece of the processes that shape planet formation is the final product. In order to understand these processes, we must first understand the range of planet properties that can be produced by those processes. Of particular interest are the planets residing at the extremes of allowable physical properties, for example, those planets at the physical inner or outer edges of systems, those planets with the largest or smallest masses, or those planets with the highest or lowest densities.

Understanding the planets that form at the extremes of measured densities will provide an important constraint on the processes of planet formation. Planetary densities are determined using measurements of planetary radii and masses. When inferred from observations, the radius depends on more detailed physical parameters: small envelope masses can greatly inflate a planet's radius (Lopez & Fortney 2013; Lissauer et al. 2014), the bulk metallicity affects the atmospheric opacity for larger planets (Burrows et al. 2007), and the temperature (Owen & Alvarez 2016) affects the atmospheric dynamics. In these ways, otherwise similar planets can have substantially different densities, and a single planet can have different values of observed radius at different times.

The least dense exoplanets observed so far appear to have bulk densities less than 0.1 g cm⁻³, which can be compared with Styrofoam (0.05 g cm⁻³). Some of these planets are referred to as super-puffs (Lee & Chiang 2016). In Figure 1, we plot the exoplanets for which bulk densities have been measured. Of the planets known to have the lowest bulk densities, some have theoretical explanations for their densities (Jontof-Hutter 2019). For example, low-density planets with short orbital periods and the largest masses (such as WASP-107 b, Anderson et al. 2017; Piaulet et al. 2021; WASP-127 b, Lam et al. 2017; HAT-P-67 b, Zhou et al. 2017; HAT-S-62 b, Hartman et al. 2019) are characterized by their large equilibrium temperatures, suggesting that their large radii are impacted by their significant stellar insolation (Guillot & Showman 2002; Laughlin et al. 2011; Thorngren & Fortney 2018) or alternatively by tidal inflation (Millholland 2019).

Zoom In Zoom Out Reset image size

A separate super-puff archetype is seen in Kepler-51, a system containing three super-puff planets of relatively small masses (Steffen et al. 2013; Steffen 2013). Low planetary densities can be observed when heat drives strongly outflowing atmospheres with small dust grains (Wang & Dai 2019), which may be produced by a photochemical haze (Gao & Zhang 2020; Ohno & Tanaka 2021). This occurs more readily for smaller, younger planets. In the case of Kepler-51, the star is young (Masuda 2014), and the transmission spectra of the inner and outer planets in the Kepler-51 system appear to support the existence of a haze layer or dusty outflow (Libby-Roberts et al. 2020).

Finally, a more perplexing system architecture is found in the Kepler-79 system (Jontof-Hutter et al. 2014) and in the the HIP 41378 system (also known as K2-93; Vanderburg et al. 2016; Becker et al. 2019; Berardo et al. 2019; Lund et al. 2019), where only one planet in the multiplanet system is a super-puff. Remarkably, these systems both subvert the "peas in a pod" archetype seen in the Kepler sample (Weiss et al. 2018; Weiss & Petigura 2020), where the ratios of adjacent planetary radii are roughly uniform within a system. This architecture, common in the exoplanet sample, is proposed to be a favorable energetic state imposed during formation (Adams 2019); as the planets in the Kepler-79 and HIP 41378 systems are both under the mass limit where this architecture is expected to develop (Adams et al. 2020), the fact that both systems containing super-puffs subvert this trend may indicate alternative formation pathways as compared to the general multiplanet sample (Jontof-Hutter 2019).

For both Kepler-79 and HIP 41378, the old age of both these host stars, combined with the fact that the super-puff is not the innermost planet in both systems, challenges theories that work for other super-puff planets. Some other possible explanations for the low observed densities of super-puff planets in general include large internal heat fluxes from ohmic dissipation (Pu & Valencia 2017) or tidal heating (Millholland 2019); planet formation that takes place farther out in the disk past the ice line (Chachan et al. 2021) in a cold, dust-free environment (Lee & Chiang 2016); or the presence of rings in a geometry that mimics the transit of a larger planet (Akinsanmi et al. 2020; Piro & Vissapragada 2020). In the case of HIP 41378 f, the combination of its long orbital period (∼542 days; Santerne et al. 2019) and subsequently low stellar insolation level, relatively old age (2.1 ± 0.3 Gyr; Lund et al. 2019), large radius (9.2 ± 0.1 R_⊕), and low mass (12 ± 3 M_⊕; Santerne et al. 2019) paints a puzzle that challenges our understanding of the least dense cold planets.

Previous work studying feasible equations of state (EOSs) for low-mass planets has indicated that very low planetary densities are possible for particular combinations of core mass and envelope properties. Batygin & Stevenson (2013) and Chen & Rogers (2016) model exemplar planets and find physically consistent solutions that can allow very different radii for a planet of a given mass, depending on the core mass and other physical properties (such as metallicity). The inclusion of a wider range of variable parameters and atmospheric properties allows for the reproduction of a larger range of mass–radius combinations, even those that initially seem to have infeasibly low densities (Lopez et al. 2012).

In this paper, we discuss where HIP 41378 f resides within the range of allowable physical planetary parameters and assess the necessity of the various hypotheses for HIP 41378 f's measured density. We do this by first, in Section 2, describing feasible current states of the system that result in a planetary mass and radius consistent with the observed values. In Section 3, we extrapolate backward and consider how planetary mass loss due to irradiation from an evolving star would affect all planets in the HIP 41378 system. In Section 4, we discuss the degeneracies on the structure of HIP 41378 f and ways that future observations might improve the limits we provide. We then conclude in Section 5 with a statement of how HIP 41378 f contributes to our understanding of the lowest-density planets.

The HIP 41378 system has a range of planetary densities, which appear to roughly decrease with orbital radius (Santerne et al. 2019). The inner two planets (HIP 41378 b and c, which have densities of 2.2 and 1.2 g cm⁻³, respectively) have densities that appear typical for the exoplanet sample. The next two planets in the system (HIP 41378 g and d) do not have precisely measured densities. The next planet, HIP 41378 e, has a slightly low density of 0.6 g cm⁻³. The outermost planet in the system, HIP 41378 f, has a bulk density of 0.09 g cm⁻³, among the lowest in the observed exoplanet sample. In this work, we do not consider HIP 41378 g, which is a low signal-to-noise ratio detection and is only detected in the radial velocities of Santerne et al. (2019). As HIP 41378 g does not have a measured radius or density, we cannot place meaningful constraints on its physical structure. A summary of the best estimates of the orbital parameters for the remaining five planets in the HIP 41378 system can be found in Table 1.

2.1. Estimated Envelope Fractions for HIP 41378 b, c, d, and e

Except for HIP 41378 f, the other planets in the HIP 41378 system have bulk densities that are easily explained by typical mass–radius relations such as Lopez et al. (2012). As the densities of HIP 41378 b, c, and e are fairly typical, it is straightforward to compute ranges of possible envelope fractions that match their observed masses and radii. HIP 41378 d has a measured radius, but only an upper limit on its mass, and so only a lower limit on the envelope fraction can be computed.

To determine the expected envelope fraction assuming rocky cores and H/He envelopes, we use the models from Lopez et al. (2012) to construct interpolated tables that allow us to model planets up to envelope fractions of 50%.

Using these grids, we find solutions for the current envelope fraction that match the observed mass and radius for each planet by interpolating the grids. The results are shown in Figure 2 overlaid with contour lines that represent the mass–radius solution for several envelope masses at 1 F_⊕ (which is roughly correct for the outer two planets and low for the inner two planets). The contours are produced for a solar metallicity and an age of ∼2 Gyr. The full interpolated solutions for the present-day envelope fractions based on the best estimates of the planet parameters are presented in Table 1, with errors that correspond to solutions that match the observational uncertainties in mass and radius.

Zoom In Zoom Out Reset image size

We note that while the mass and radius measurements reported for HIP 41378 b and c are fairly secure, the measured mass for HIP 41378 e is based on an estimate rather than a direct detection of the orbital period (Santerne et al. 2019), as the planet transited only once in the two campaigns of available K2 data (Becker et al. 2019). This value may need to be revised if an additional transit observation is made.

2.2. A Hydrostatic Solution for HIP 41378 f

HIP 41378 f's density is less consistent with the range of accepted planetary densities, as HIP 41378 f has a radius of 9.2 ± 0.1 R_⊕ (Santerne et al. 2019) and a mass of 12 ± 3 M_⊕ (Santerne et al. 2019). The combination of the mass of a mini-Neptune and the radius of a Jupiter-like planet creates a compositional puzzle. Upon a glance at Figure 2, it is clear that HIP 41378 f has a much larger envelope fraction compared to the other four modeled planets. Indeed, this planet cannot be well reproduced by the relations used above, as its measured radius is larger than any of the models from Lopez et al. (2012).

To determine whether it is necessary to invoke a nonstandard explanation for HIP 41378 f's density, we attempted to solve for feasible parameters that create a hydrostatic planet while matching the present-day radius and mass observed in Santerne et al. (2019) for HIP 41378 f. To explore what these parameters might be, most importantly its core mass, we solved for planetary hydrostatic equilibrium using the code and method of Howe et al. (2014) and Howe & Burrows (2015). We fix the planetary mass to be 12 M_⊕, consistent with observations, and then numerically compute radii for that mass that correspond to various interior structures. We use an updated version of the code of Howe et al. (2014), ⁵ which allows us to compute mass–radius relations for an extensive range of potential planet parameter combinations, thus offering a way to test the extreme values that a super-puff planet's parameters typically have.

The details of the numerical method are given in the appendices of Howe et al. (2014). In short, we solve the equations of hydrostatic equilibrium using a fourth-order Runge–Kutta method with a combination of two EOSs: an EOS for one of three possible core compositions (a pure Fe core, a perovskite core, or an ice VII core) and the Saumon et al. (1995) tabular EOS for an H₂–He envelope. The Vinet EOS used in the core is defined as

$$\begin{eqnarray}\begin{array}{rcl}P & = & \ 3{K}_{0}{\left(\rho /{\rho }_{0}\right)}^{2/3}\left(1-{\left({\rho }_{0}/\rho \right)}^{1/3}\right)\\ & & \times \exp \left[1.5({K}_{0}^{{\prime} }-1)\left(1-{\left({\rho }_{0}/\rho \right)}^{1/3}\right)\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where P is the pressure, and the numerical coefficients are a reference density at zero pressure ρ₀, bulk modulus K₀, and the pressure derivative ${K}_{0}^{{\prime} }$. For the iron core, ρ₀ = 8.267 g cm⁻³, K₀ = 1.634 Mbar, and ${K}_{0}^{{\prime} }\,$ = 5.38 Mbar (Dewaele et al. 2006); for the perovskite core, ρ₀ = 4.064 g cm⁻³, K₀ = 2.48 Mbar, and ${K}_{0}^{{\prime} }\,$ = 3.91 Mbar (Tsuchiya et al. 2004); and for the ice core, ρ₀ = 1.4876 g cm⁻³, K₀ = 1.49 Mbar, and ${K}_{0}^{{\prime} }\,$ = 6.2 Mbar (Wolanin et al. 1997). For a range of planet envelope entropy values and envelope metallicities, we solve the equations of hydrostatic equilibrium with the above EOSs to determine whether feasible solutions exist and what the radii of those solutions should be.

Howe et al. (2014) found that using an Fe/rock core as opposed to an ice core does not generally significantly change the allowable envelope fraction for planets of a given mass and radius (see their Table 2—for planets with large envelope masses, the envelope mass may change by 1 M_⊕ dependent on the core type). In this work, we test three possible core compositions for HIP 41378 f under the assumption that there is no metallicity gradient in the envelope.

Table 2. Feasible Values for HIP 41378 f's Core Mass Given Varying Combinations of Entropy and Metallicity That Result in Planets with the Observed Present Physical Radii within 1σ of the Observationally Determined Value (9.2 ± 0.1 R_⊕; Santerne et al. 2019)

Model	Entropy	Metallicity	Allowable Core Mass	Allowable Core Mass	Allowable Core Mass
	(kb /baryon)	(Zsolar)	(Iron) (M⊕, 1σ)	(Perovskite) (M⊕, 1σ)	(Ice) (M⊕, 1σ)
1	6.0	1.0	[0.1, 0.4]	[0.1, 0.4]	[0.1, 0.5]
2	6.5	10.0	[1.7, 2.0]	[1.8, 2.1]	[2.0, 2.3]
3	6.5	1.0	[2.4, 2.7]	[2.6, 2.9]	[2.7, 3.1]

Download table as:  ASCIITypeset image

To explore the potential mass–radius relations for a planet of HIP 41378 f's mass, we computed six main models for each core type under the above prescription, each with a varying combination of entropy (either 5.5, 6.0, or 6.5 k_b /baryon, the range considered in Howe et al. 2014) and metallicity (either solar or 10 times solar). For each of these six models, we evaluated the interior structure of a planet with 1 of 14 discrete core mass fractions. For each combination of entropy, metallicity, core composition, and core mass, we solve for the hydrostatically stable solution and derive the planetary radius that corresponds to that solution if it exists. In Figure 3, we show interpolated versions of the core mass—planetary radius relationship for the case of a planet with an iron core. All points were computed for a planet with a mass of 12 M_⊕. This result shows a wide range of possible values for the radius of the planet. The primary driver of the variability in allowed radii is the core mass, with larger core masses corresponding to planets with overall smaller radii. Secondary factors that affect the planetary radius are the entropy and metallicity of the envelope, with higher envelope entropy corresponding to puffier planets and higher metallicity corresponding to slightly smaller radius planets.

Zoom In Zoom Out Reset image size

While these relationships are intuitive, the important result is the point at which the observed radius can be reproduced and the physical consistency of those solutions with the core accretion paradigm. By evaluating the range of HIP 41378 f's core mass within its given radius 9.2 ± 0.1 R_⊕, we found the range of core masses of HIP 41378 f that reproduce the observed planetary density.

The allowable parameters for each model and each core type are summarized in Table 2. The iron and perovskite core compositions result in very similar ranges of allowable core masses that result in planets with radii consistent with the observationally derived value. For each model, the ice composition allows a slightly more massive core. Among the three models, Model 3 (k = 6.5 k_b /baryon and Z = 1 × solar) allows for the largest possible core masses, so we consider this model moving forward. The 1σ upper limit on core mass is 2.7, 2.9, and 3.1 M_⊕ for the iron, perovskite, and ice cores, respectively. Moving forward, we take a 25% core mass, 75% envelope mass to be representative of these solutions; however, it is possible that the core mass is much less massive (as would be the case in Models 1 and 2).

In Section 2, we found plausible planetary structures for all five planets with measured masses and radii that are consistent with the present-day observed values. Due to the fact that the inner two planets have the highest densities in the system and the importance of photoevaporation in sculpting the radius distribution of close-in planets (Sanz-Forcada et al. 2011; Owen & Wu 2017; Uzsoy et al. 2021), we next consider whether the system once consisted of planets with higher envelope mass fractions that diminished over time owing to the radiation from their host star, HIP 41378, as it evolved. In this scenario, the innermost planets would be most strongly affected by the stellar X-ray and ultraviolet (XUV) radiation and experience the largest fractional loss of the envelopes, potentially resulting in the density gradient that we see today.

3.1. General Analysis of Mass Loss for HIP 41378 b, c, d, e

We used the program VPLanet ⁶ (Barnes et al. 2020) to model the evolution of the five planets with measured masses and radii. VPLanet allows us to solve how the envelope mass of each planet in the system changes as the star evolves. The mass of a planet will decrease over time owing to XUV radiation from the host star. XUV radiation from the host star will heat a planetary atmosphere via the ionization of atmospheric H/He, and the heat will drive atmospheric outflows (Lammer et al. 2003; Vidal-Madjar 2003; Owen & Jackson 2012; Murray-Clay et al. 2009). This process will occur at a higher rate with a higher incident XUV flux. In Figure 4, we show the evolution of the luminosities as the star ages, as computed by VPLanet for a 1.16 M_⊙ star (which is the mass of the host star HIP 41378; Santerne et al. 2019) and based on the stellar models of Baraffe et al. (2015). The XUV luminosity, which drives mass loss, is strongest when the star is young.

Zoom In Zoom Out Reset image size

We compute the mass loss due to XUV radiation over two billion years using the atmesc and stellar modules along with the parameters set in Table 1. We use variable time stepping with an accuracy coefficient of 0.01 and integrate between 5 × 10⁶ yr and 2 × 10⁹ yr using the Baraffe et al. (2015) stellar model, and we assume a saturated XUV luminosity fraction of 10⁻³, which persists for the first 0.1 Gyr of the star's lifetime (Lammer et al. 2014). The planetary semimajor axis and mass values are set to the best-fit values from Santerne et al. (2019) and summarized in Table 1. We use VPLanet's automatic atmospheric escape prescription, which switches between radiation/recombination-limited (Luger et al. 2015) and energy-limited escape according to the incident XUV flux. In contrast to energy-limited escape, which occurs at lower incident FUV fluxes, some of the inputted energy is lost via radiative cooling in recombination-limited escape (Murray-Clay et al. 2009). The VPLanet analysis works in terms of the planetary mass, with the mass–radius models of Lopez et al. (2012) being subsequently used to convert the mass into a radius. As discussed in Section 2, the Lopez et al. (2012) mass–radius models describe well the structures of all planets except for HIP 41378 f. We note that the radius cannot accurately be computed for HIP 41378 f using this model, an issue we will return to in Section 3.2. We also note two additional caveats: (1) we assume that the stellar X-ray and EUV fluxes will evolve on the same timescale, but this may not always be true (King & Wheatley 2021); and (2) the efficiency factor for photoevaporative mass loss will depend on the specific level of XUV radiation and the planet gravity (Caldiroli et al. 2021), so further refinement of the planet masses or stellar luminosity could alter our results.

While we have the final present-day values for the mass, radius, and expected envelope fraction for each planet, we want to test what feasible envelope fractions could have been at formation. To test this and determine which HIP 41378 planets may have had significantly higher envelope mass fractions at formation, we perform a parameter space survey where we vary the estimated initial envelope masses for each planet in the system and then compute the mass loss over 2 Gyr, attempting to match the current value. Once we found which initial envelope mass fractions led to the integration output matching present-day observed values, we could then estimate by how much each planet's envelope mass fraction likely diminished over two billion years, the estimated current age of the system (Lund et al. 2019). Once the parameter space survey yielded a result consistent with the currently measured system parameters, we computed a full system evolution with VPLanet. The resultant evolutions of the planetary radii, envelope masses, and densities are shown in Figure 5.

Zoom In Zoom Out Reset image size

We also show a zoomed-in view of the planetary envelope mass fraction over time in Figure 6. We observe that the two inner planets HIP 41378 b and HIP 41378 c had initial envelope masses that substantially decreased owing to the XUV luminosity from the system's host star. HIP 41378 d, e, and f lose very minimal amounts of mass from their envelopes as the star evolves, due to their larger orbital separations from the host star. Our VPLanet analysis supports our initial hypothesis that the two inner planets in the HIP 41378 system likely had low densities in the past and lost some of their envelopes owing to the XUV radiation from the host star. However, the amount of XUV radiation on the atmospheres of HIP 41378 b and c is not expected to have been enough to drive the evolution from an initial super-puff density to the present-day value. Instead, it is likely that HIP 41378 b and c had at formation envelope fractions a factor of two larger than their present-day values. Similarly, the middle planets in the system (HIP 41378 d and HIP 41378 e) likely did not experience much XUV-driven photoevaporation, and their present-day observed envelope fractions are likely very similar to the values at formation.

Zoom In Zoom Out Reset image size

3.2. Detailed Analysis of HIP 41378 f's Evolution

To further study HIP 41378 f's nonstandard interior structure, we perform additional explorations of its time evolution using coupled interior structure and atmospheric modeling of planets with Earth composition rocky cores and primordial hydrogen−helium envelopes (Rogers et al. 2011). We performed these calculations using the Modules for Experiments in Stellar Astrophysics (MESA) code (v12778; Paxton et al. 2011, 2013, 2015, 2018, 2019) and with the additions outlined in Malsky & Rogers (2020). The framework from Malsky & Rogers (2020) specifically allows for the modeling of sub-Neptune-mass exoplanets over a range of envelope mass fractions, planet core masses, and mass-loss parameterizations. By matching the observed properties of HIP 41378 f to these 1D structure models, we can characterize the interior structure of the planet and how it evolves over the course of billions of years.

We parameterized the 1D evolution model in MESA as follows. Based on the results of the previous section, which showed that the mass of HIP 41378 f did not significantly change over the system age, we choose the initial mass and envelope fraction to be equal to the values found above. As such, we assume an initial planet mass of 12.0 M_⊕, an initial envelope mass fraction of 0.75, an orbital separation of 1.37 au, and initial abundances of X = 0.74, Y = 0.24, and Z = 0.02. We set the outer atmospheric boundary layer within MESA at τ = 2/3 and implement a gray Eddington T(τ) relation with the atm_T_tau_relation option. This is the outermost zone for which MESA solves the set of coupled structure and composition equations. We then extrapolate the atmosphere up to 1.0 mbar assuming an isothermal temperature profile and a constant atmospheric mean molecular mass. The pressure level of 1.0 mbar corresponds to the radii observed by transit surveys (Miller et al. 2009). All other model parameters are identical to the full model description in Malsky & Rogers (2020), except for the removal of the "hot start" evolution step. This step was required in Malsky & Rogers (2020) in order to standardize the starting entropy of highly irradiated super-Earths, and it is not required for HIP 41378 f's evolution.

Our 1D models of HIP 41378 f are consistent with a bulk density of 0.087 and match observational results. By 2.1 Gyr, the model had a transit radius of 9.1 R_⊕, consistent to 1σ with the measured value from Santerne et al. (2019), and as expected from the previous analysis, the planet did not exhibit significant mass loss over its lifetime. The approximate lower bound 75% envelope fraction derived from the hydrostatic solutions does reproduce the observed planet radius through the MESA modeling.

Another concern for the atmospheric retention of HIP 41378 f is the possibility that the planet would have lost its atmosphere early owing to Parker winds (Parker 1958). For some planets, the Parker wind causes rapid boil-off of their atmospheres on kiloyear timescales (Owen & Wu 2016). For super-puff planets, Parker winds can potentially prevent the retention of a planetary atmosphere.

To test whether Parker winds could have caused HIP 41378 f to have lost its atmosphere in a problematically short timescale, we compute the mass-loss rate using Equation (5) of Chachan et al. (2020):

$$\begin{eqnarray}&&{dM}/{dt}=4\pi {r}_{s}^{2}{c}_{s}{\rho }_{p}\exp (3/2-2{r}_{s}/{R}_{p}),\end{eqnarray} \tag{ 2 }$$

where ${r}_{s}={{Gm}}_{p}/2{c}_{s}^{2}$ is the sonic radius, m_p the planetary mass, ${c}_{s}=\sqrt{{k}_{b}T/\mu }$ the isothermal sound speed, μ ≈ 2.2 amu the assumed atmospheric mean molecular weight, ρ_p the atmospheric density at the radius of the planet, and R_p the planetary radius. For HIP 41378 f, the equilibrium temperature is T ≈ 300 K (Santerne et al. 2019), and we compute the speed of sound to be c_s = 1 km s⁻¹, sonic radius to be r_s = 0.016 au, and ρ_p = 10⁻³ g cm⁻³. Using these values, we find that the Parker wind mass-loss rate is ∼10⁻⁷ g s⁻¹, corresponding to an atmospheric mass-loss timescale longer than the Hubble time.

HIP 41378 f's large mass (12 M_⊕) causes the Parker wind mass-loss rate to be small enough that this planet would not have been expected to lose a substantial amount of its atmosphere via early boil-off. However, it is important to note that while the Parker wind mass loss is not important for HIP 41378 f's evolution, the same may not be true for other similar planets; if HIP 41378 f's mass had been 25% of its current value, for example, its mass-loss rate would been 10¹⁸ g s⁻¹ and the atmosphere would have been lost on kiloyear timescales. In this way, HIP 41378 f's relatively large mass allowed it to retain its atmosphere (similarly to other higher-mass, larger gas-to-core ratio planets as in Vissapragada et al. 2020).

In this work, we have presented self-consistent solutions for HIP 41378 f's interior structure (including core mass and envelope entropy) both today and after formation (2 Gyr ago). We have also derived estimates of the present-day envelope fractions for HIP 41378 b, c, d, and e and used VPLanet to simulate their evolution while accounting for the evolving XUV flux, and thus we estimated the initial envelope fractions. We find that while HIP 41378 b and c (with periods of roughly 15 and 30 days, respectively) likely had significantly (by a factor of ∼2) higher envelope fractions at formation, none of the planets except HIP 41378 f formed with anomalously low densities.

Planets with high H/He envelope mass fractions may experience mass loss due to large XUV fluxes (Lammer et al. 2003). This is most relevant for planets residing close to their host stars (Owen & Alvarez 2016), in which case substantial initial envelope fractions may prevent the destruction of a low-density planet's planetary atmosphere (Hallatt & Lee 2022), particularly for short-period planets. Short-period planets will also be affected by mass loss due to Parker winds early in their lifetimes; for some super-puff planets, this could lead the entire atmosphere to be lost on ∼kiloyear timescales (Chachan et al. 2020). For HIP 41378 f, the predicted mass-loss rate is predicted to be very small, meaning that mass loss due to isothermal Parker winds will be negligible. While this rate could be an underestimate if HIP 41378 f has significant internal heat flux, that is unlikely for a planet with such a large orbital radius (Pu & Valencia 2017). In cases like HIP 41378 f, where the planet's orbital separation is large enough that mass loss will be minimal, the primary factor driving the change in physical radius from formation to today is the change in the envelope's entropy as the planet gravitationally contracts. At formation, cold planets with masses around 10 M_⊕ can have entropies ranging between 9 and 9.5 k_b /baryon (Mordasini et al. 2017). Our investigations with PlanetSolver and the MESA simulations show that at 2 Gyr HIP 41378 f can have a radius and density consistent with the observed values if the planet has an envelope entropy of 6.5 k_b /baryon.

Our results show that a high envelope fraction (around 75%) and envelope entropy of 6.5 k_b /baryon are also consistent with HIP 41378 f's measured planetary mass and radius. A transient enhancement in planetary envelope entropy (which increases the effective radius of a planet) could be recent collisions from small rocky planetesimals or even from other planets (Ketchum et al. 2011; Anderson & Adams 2012). However, the necessary mass flux of solid materials required to enhance the planetary entropy sufficiently at present day (2 Gyr) would be several orders of magnitude larger than the current mass of the planet, and so this explanation is not feasible. Similarly, we can assess the feasibility of some of the alternative hypotheses for super-puff planetary densities. Due to HIP 41378 f's relatively large (1.3 au) orbital distance, theories that depend on stellar insolation such as tidal inflation (Millholland 2019) or ohmic dissipation (Pu & Valencia 2017) can be dismissed for this particular planet.

Transmission spectroscopy can further differentiate between the remaining feasible theories, one of which is an outflowing atmosphere with a large population of small grains (Wang & Dai 2019) or photochemical hazes (Gao & Zhang 2020), which may produce dust grains (Ohno & Tanaka 2021). This explanation, although more likely to affect younger planets, will be testable with transmission spectroscopy of HIP 41378 f in the future, as such dusty atmospheres are expected to have flat transmission spectra (Kawashima & Ikoma 2019). However, we compute very low expected mass-loss rates for HIP 41378 f via both photoevaporation and Parker winds, which makes it unlikely that dust would exist at altitudes sufficient to obscure transmission spectra. In the case of Kepler-79, water features are not detected in the transmission spectrum (Chachan et al. 2020), supporting the haze hypothesis. In contrast, if transmission spectra show evidence for water features in HIP 41378 f's atmosphere, that could be evidence in favor of the hypothesis of Lee & Chiang (2016; see also Lee et al. 2018), in which super-puffs form dust-free in a cold formation environment far from the host star. Indeed, the largest core mass is possible with an ice core, and if HIP 41378 f did have an ice core, this would place additional constraints on its formation process. An ice core requires very particular formation conditions: the planet likely would have formed past the ice line in the disk and then migrated inward—an origin possibly consistent with the scenario outlined in Lee & Chiang (2016). However, recent Hubble Space Telescope (HST) results from Alam et al. (2022) suggest no evidence for water features in HIP 41378 f's transmission spectrum.

A core mass fraction of 25% or less for a planet of HIP 41378 f's mass also appears immediately inconsistent with the paradigm of giant planet formation. Planets this small are unlikely to form via disk instability (Boss 1997; Dodson-Robinson et al. 2009; Boley 2009), so the most likely pathway to a planet like HIP 41378 f is via core accretion. In the core accretion model, cores form via accretion of planetesimals. Following the formation of the core, the planet will accrete nearby gas until the mass of the core and the mass of the envelope are roughly equal, at which point runaway accretion can occur and the envelope may become far more massive. Barring the absence of a sufficient reservoir of gas, the outcome of runaway accretion is generally a Jupiter-mass planet (Pollack et al. 1996). All of our hydrostatic solutions require envelope fractions for HIP 41378 f of 75% or greater, and for a planet with this core mass, the gas fraction will be set by how much the core can accrete during the disk lifetime (Lee 2019). To create a planet with a mass of only 12 M_⊕ would subsequently require the termination of accretion slightly after the start of runaway accretion (before HIP 41378 f would have had the chance to grow to a Jupiter mass). This issue could be solved with particularly serendipitous (and low-probability) timing for the disk dissipation, a solution that has also been invoked for the paucity of close-in Jupiter-mass planets (Lee et al. 2014). However, evidence suggests that the lifetime of the protoplanetary disk is not generally the primary driver of the end of planetary mass accretion (Adams et al. 2021)

There is a second, more fatal flaw with the core accretion picture for HIP 41378 f: the allowable core mass we compute is smaller than the limit required to accrete a significant (>50%) envelope fraction, no matter where in the disk the planet formed. The critical core mass for runaway accretion is estimated to be upward of 5–10 M_⊕ (Papaloizou & Terquem 1999; Rafikov 2006; Coleman et al. 2017), but we compute a maximum core mass of 3 M_⊕, with the value being only 0.4 M_⊕ under more standard assumptions about the planetary entropy. Since this is well below the threshold required to trigger runaway accretion, there is no easy way to explain the 75%+ envelope fraction in the core accretion paradigm.

In this paper, we attempted to find a feasible structure for HIP 41378 f that does not require invoking the existence of planetary rings. Although hydrostatic solutions do exist, these solutions are not consistent with the paradigm of core accretion due to the combination of core mass and envelope fraction. Because of this, there is no easy explanation for how HIP 41378 f could have formed with its measured parameters. Our analysis provides support for the hypothesis of Piro & Vissapragada (2020) and Akinsanmi et al. (2020) that HIP 41378 f's anomalously low density could be due to the presence of planetary rings. As discussed in Piro & Vissapragada (2020), HIP 41378 f is the only known super-puff planet for which the rings hypothesis could be assessed using ground-based telescopes. Simultaneously with the review process of this manuscript, Alam et al. (2022) announced HST results showing a flat transmission spectrum for HIP 41378 f. Ohno & Fortney (2022) also demonstrated that flat transmission spectra are expected for planets with ringed geometries, supporting the potential for rings in the case of HIP 41378 f. It is important to continue to monitor HIP 41378 f for transit timing variations and refine its transit ephemeris (as in Bryant et al. 2021) so that future James Webb Space Telescope (JWST) observations can be made.

This work considers both the current envelope fractions and feasible past evolutions of envelope fractions for the planets in the HIP 41378 system, with a particular focus on the super-puff HIP 41378 f. We find estimated present-day envelope fractions by solving the equations of hydrostatic equilibrium and then use VPLanet to deduce their envelope fractions at formation. We find that the inner two planets (HIP 41378 b and c) did not have anomalously low densities at formation but HIP 41378 f did.

Although HIP 41378 f's structure can be explained by a combination of larger envelope entropy values and high envelope fractions, there remain problems with the theoretical interpretation of those values. Namely, as discussed in this paper, consistent solutions for HIP 41378 f's interior structure require present-day envelope fractions upward of 75%, which provides a tension with the core accretion paradigm. Because of this, ways in which the planetary radius may have been overestimated by the transit data should be considered, including the theory advanced in Piro & Vissapragada (2020) and Akinsanmi et al. (2020) that the anomalously large radius of HIP 41378 f could be due to planetary rings. This particular theory appears consistent with current data (Alam et al. 2022; Ohno & Fortney 2022) and will be directly testable with additional observations with JWST.

M.B. thanks the University of Michigan Undergraduate Research Opportunity and Research Scholars programs for their support of this project. J.C.B. has been supported by the Heising-Simons 51 Pegasi b postdoctoral fellowship. We thank Rodrigo Luger and Rory Barnes for useful conversations and assistance with VPLanet. We also thank the anonymous referee for their invaluable insight into our paper, which led to substantial improvements. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Facility: Exoplanet Archive. -

Software: pandas (McKinney 2010), IPython (Pérez & Granger 2007), matplotlib (Hunter 2007), scipy (Jones et al. 2001), numpy (Oliphant 2006), Jupyter (Kluyver et al. 2016), VPLanet (Barnes et al. 2020), PlanetSolver (Howe & Burrows 2015).