© The Authors 2025

1 Introduction

Sub-Neptune exoplanets with a radius between 2 R_⊕ and 4 R_⊕ make up the most abundant planet class found by currently available surveys (e.g., Borucki et al. 2011; Batalha et al. 2013; Fulton et al. 2017). The absence of a Solar System analog highlights the importance of examining sub-Neptune exoplanets to understand this planet class. However, observing sub-Neptunes is challenging with current technology, as these small-radius planets only cover a fraction of their host stars’ light during their transit, leading to signal strengths of only a few hundred parts per million (ppm). Therefore, facilities with high sensitivities are needed to distinguish the planet signal from its star and various noise sources. The required precision is best met by space-based facilities such as the Hubble Space Telescope (HST) and the James Webb Space Telescope (JWST). But even these telescopes operate close to their limits to detect sub-Neptune planetary atmospheres.

The few sub-Neptune atmospheres characterized today hint at very diverse atmospheric conditions on these planets. Their compositions range from hydrogen-helium dominated atmospheres on GJ 3470 b (Benneke et al. 2019), HD 3167 c (Guilluy et al. 2021), K2-18 b (Madhusudhan et al. 2023), and TOI-421 b (Davenport et al. 2025) to metal-rich, miscible atmospheres as detected on TOI-270 d (Benneke et al. 2024). The molecular inventory differs across these planets: CH₄ is the dominant absorber in K2-18 b and TOI-270 d, H₂O has been seen in TOI-421 b, and SO₂ absorption has been detected in GJ 3470 b (Beatty et al. 2024) and is thought to be a consequence of disequilibrium photochemistry. Furthermore, many of those atmospheres possibly host aerosols, as detected in GJ 3470 b (Benneke et al. 2019; Beatty et al. 2024) and GJ 1214 b (Kreidberg et al. 2014a; Kempton et al. 2023).

In contrast to the well-characterized planets, many of the observed sub-Neptune spectra are featureless or have muted spectral features. Examples of these include HD 97658b (Knutson et al. 2014; Guo et al. 2020), HD 106315 b (Kreidberg et al. 2022), TOI-836c (Wallack et al. 2024), GJ 1214 b (Kreidberg et al. 2014a; Kempton et al. 2023; Schlawin et al. 2024), GJ 9827 d (Roy et al. 2023; Piaulet-Ghorayeb et al. 2024), and GJ 3090 b (Ahrer et al. 2025), where muted features are detected for the last three planets. These muted spectral features were predicted as a consequence of high mean molecular weight atmospheres (e.g., Miller-Ricci et al. 2009) or high-altitude aerosols (e.g., Benneke & Seager 2013; Morley et al. 2013). Which of these effects dominates is still unclear for most of the observed planets. Trends in the spectra of Neptune-sized to sub-Neptune-sized planets suggest that there is a correlation between planet equilibrium temperature and atmospheric feature size (Crossfield & Kreidberg 2017; Yu et al. 2021; Brande et al. 2024). Brande et al. (2024) attribute this trend to clouds, which rain out at low temperatures below 500 K. For temperatures above 800 K, a reduced formation of organic hazes may lead to clear atmospheres (e.g., Zahnle et al. 2009; Fortney et al. 2013; Crossfield & Kreidberg 2017; Gao et al. 2020; Yu et al. 2021). The characterized sample of sub-Neptunes is small, and not all planets follow this trend. More systematic observations are needed to fully understand which environment favors which atmospheric conditions.

The Sub-neptune Planetary Atmosphere Characterization Experiment (SPACE) with the HST is a multi-cycle treasury program that aims to conduct such a systematic survey (Kreidberg et al. 2022). It probes a sample of eight sub-Neptunes with near-infrared transmission spectroscopy using the Wide Field Camera 3 (WFC3¹) and combines these data with ultraviolet (UV) stellar characterization using the Space Telescope Imaging Spectrograph (STIS²). As shown in the top panel of Fig. 1, the targets span a physically motivated grid across planet radius (2–3.5 R_⊕) and temperature (300–1400 K), which allows to examine how the atmospheric conditions are altered with these parameters. Planet growth models (Zeng et al. 2019) predict the presence of atmospheres on all SPACE targets based on their radii and masses between 3 M_⊕ and 13 M_⊕ (see bottom panel of Fig. 1). The associated host stars range in spectral type between M and G, corresponding to the types with the most planet detections today. The experiment aims to reveal how sub-Neptune atmospheres are shaped by metal enrichment, disequilibrium chemistry, and aerosols. Furthermore, it guides upcoming observations with JWST and identify sub-Neptunes with clear atmospheres.

Here, we examine HD 86226 c, which is a 2.2 R_⊕ sub-Neptune that orbits its G-type host star on a 3.98 day orbit. It has a high equilibrium temperature of 1310 K (Teske et al. 2020), and its low mass of 7.25 M_⊕ suggests the presence of an atmosphere (Teske et al. 2020) on the planet (see Fig. 1). These characteristics place HD 86226 c close to both the “radius gap” (Fulton et al. 2017) and the “hot Neptune desert” (Mazeh et al. 2016). Moreover, HD 86226 c has a long-period giant planet companion, HD 86226b (Arriagada et al. 2010). These factors make HD 86226 c an interesting candidate for atmospheric characterization, as it could inform planet formation and evolution models. The G-type host star is both beneficial and hindering for the observations. While we do not expect high stellar activity based on its estimated age of 4.6 Gyr (Teske et al. 2020), the expected transit depth is only on the order of 400 ppm and therefore just on the lower edge of what can be observed with HST. However, in this sparsely populated part of the temperature-radius parameter space, HD 86226c is one of the best available targets for this study.

The hot temperatures on HD 86226 c suggest a high chance for a clear atmosphere. At these temperatures, CO is the dominant carbon reservoir, instead of CH₄. This hinders the formation of organic hazes (Zahnle et al. 2009; Fortney et al. 2013; Gao et al. 2020), which could mute the observed atmospheric features. Previous observations of HD 86226 c with the Transiting Exoplanet Survey Satellite (TESS) and the 6.5 m Magellan Telescopes (Kokori et al. 2023; Teske et al. 2020) constrained its orbital parameters and mass. Possible volatile-rich atmospheric compositions were recently predicted by Piette et al. (2023), but no spectroscopic observations of the planet’s atmosphere were conducted prior to this experiment.

Section 2 provides details on the observations of HD 86226 c and its host star. The independent data reductions of the transit spectra and stellar data are presented in Sect. 3. The resulting light curves are analyzed in Sect. 4, where we also show the planetary spectrum. Forward models and retrievals of the planet spectrum are detailed and shown in Sect. 5. The implications of these results are discussed in Sect. 6, and we summarize our findings in Sect. 7.

Fig. 1 Sub-Neptune targets of the SPACE program. HD 86226c is shown with a green square marker, and the other SPACE targets are shown in purple. Top: temperature-radius plane. Gray circles show the known population of planets with radii between 1.8 R⊕ and 4 R⊕ based on the entries of the NASA Exoplanet Archive3 in May 2024. Black crosses mark JWST Cycle 1–4 targets, except for the SPACE targets TOI-431 d and HD 86226 c, which will be observed in Cycle 4. Bottom: mass-radius plane. Color curves show models from Zeng et al. (2019) for various planetary compositions and temperatures.

2 Observations

2.1 HST observations with WFC3 and STIS

HD 86226 c was observed with HST as part of the Cycle 30 Large Treasury Program 17192 (PI: Laura Kreidberg). For the transit observations, we used WFC3 with the G141 grism, spanning the wavelength between 1.1 and 1.7 μm. Individual exposures of the time series observations used the MULTIACCUM mode with three nondestructive readouts and a total integration time of 46.696 s. The spatial scanning mode was used to allow for longer integration times before saturation by spreading the light across the detector (see e.g., Deming et al. 2013). The scan rate of 0.754 “s⁻¹ yields a scan height of 35.2” for the final readout, corresponding to approximately 290 pixels on the detector.

Roundtrip spatial scans (alternating between forward and reverse) were used to maximize the observing time per orbit, leading to 26 exposures observed per HST orbit. In addition, a direct image of HD 86226 was obtained at the start of each orbit using the F139M filter and an exposure time of 5.118 s. We achieved an observing efficiency of 42%, including overheads and target acquisition time.

Nine transits were observed, with four HST orbits per transit between May 2023 and January 2025. A guiding failure occurred on 2024 January 10, when the fine guidance sensor could only acquire one guiding star. The telescope observed orbit 29 on gyro control, which resulted in the spectrum moving by up to one pixel row on the detector. As the flux measured in this dataset is correlated with the position of the spectral trace on the detector (see Sect. 3.1), the uncertainty of the flux in orbit 29 is increased. However, the data of this orbit generally remain usable for transit spectroscopy.

The host star HD 86226 was also observed with STIS over two orbits in May 2023. These observations consist of three exposures using the G140M, G140L and E230M gratings, covering 1194–1248 Å (including the H I Lyman α 1215 Å line), 1160–1715 Å and 2274–3118 Å respectively, with exposure times of 1803.169 s, 1978.177 s and 200.020 s.

2.2 Host star monitoring

In addition to HST spectroscopy, we also conducted long-term photometric monitoring of the host star HD 86226 from the ground to verify the level of stellar activity. On more than 60 observing nights from March 2023 to April 2024, HD 86226 was observed with the automated 24-inch PlaneWave CDK telescope at Wesleyan University’s Van Vleck Observatory in Connecticut, USA. The telescope has an aperture of 0.61 m and is equipped with a FLI PL4240 camera. The field of view is approximately 23 arcmin × 23 arcmin at approximately 0.7 arcsec per pixel. Typically, twenty exposures were taken every observing night in each of the BVRI bands, with exposure times of 6.0 s, 5.0 s, 2.0 s, and 1.5 s, respectively. We present the reduction of this dataset and the photometric monitoring results in Sect. 3.5.

3 Data reduction and analysis

We reduced the HST/WFC3 data with two independent pipelines: PACMAN (Zieba & Kreidberg 2022) and Eureka! (Bell et al. 2022). Sect. 3.1 covers the data reduction with PACMAN, a well-tested (e.g., Kreidberg et al. 2014a,b, 2018a,b; Colón et al. 2020) pipeline optimized for HST/WFC3 data. Details about the Eureka! analysis can be found in Section 3.2. The stellar parameters used in these analyses are estimated with ARIADNE in Sect. 3.3. Section 3.4 describes the reduction of the STIS data of HD 86226, and the monitoring results can be found in Sect. 3.5.

3.1 Transit reduction with PACMAN

The PACMAN reduction used the _ima files provided by the Barbara A. Mikulski Archive for Space Telescopes (MAST⁴) as a starting point. The time stamps of these files were corrected for the motion of the HST with respect to the Solar System Barycenter.

The spectrum for each exposure was derived from the three non-destructive up-the-ramp readouts using an optimal extraction algorithm (Horne 1986). To extract each spectrum, PACMAN created a sequence of differenced images consisting of the difference between subsequent up-the-ramp samples. The differenced images contain the first order of the spectral trace (columns 153–334). Each differenced image was background corrected by subtracting the median value of all pixels with a value below 500 counts. The optimal extraction algorithm was applied to all rows in the differenced image within 40 rows of the edges of the spectral trace in the spatial direction. These edges were determined by accumulating the flux in the spatial direction and identifying the rows between which the flux declines the strongest. The final spectrum and its uncertainty per exposure were determined by summing the spectra and variances for all of the differenced images.

For wavelength calibration, we used the stellar model grids of Castelli & Kurucz (2003). The model assumes the stellar parameters log(Z) = $-{0.07}_{-0.14}^{+0.15}$ and log₁₀(g[cm s⁻²]) = ${4.55}_{-0.48}^{+0.40}$, obtained from fits to the stellar spectral energy distribution (SED) with ARIADNE, which are described in Sect. 3.3. To obtain the calibration, we interpolated the stellar model and fitted it to the extracted observed spectrum. In this fit, we derived the shift in wavelength direction and allowed for a linear stretch in flux and wavelength.

The raw HD 86226 c light curve shows a non-constant offset between spatial scans in the forward and reverse direction. While the known upstream-downstream effect (McCullough & MacKenty 2012) more typically results in a constant offset between scan directions, a similar non-constant offset was also seen in WFC3 observations of 55 Cancri e (Tsiaras et al. 2016). This behavior might be a consequence of the fast scan speed used for observing bright objects: For fast scans, the position of the spectrum in the spatial direction can shift by up to two pixels between orbits (see Fig. 2). Because the detector is read out row-by-row and has no shutter, a shift in spatial position effectively changes the exposure time. This introduces a systematic change in flux with the position of the trace in the spatial (row) direction on the detector.

To correct this effect, we introduced a new systematics decorrelation parameter, the “row shift”, based on the position of the spectral trace. We independently assumed a linear correlation between flux and shift of spectral trace in spatial direction (row shift) for both scan directions. The row shift was determined for each exposure by summing the flux of the first non-destructive readout along the dispersion direction (columns) and comparing this spatial profile to the profile of the first exposure of all observations. This resulted in the row shift relative to the first exposure, which was used in the further analysis of the data.

Figure 2 shows that the spatial extent of the spectral trace on the detector is correlated with the position of the trace. Analogous to the row shift measurement, the relative trace extent was derived from the spatial profiles of the last detector readouts of each exposure. The profile of each exposure was compared to the profile of the first exposure of the observations by fitting a shift and a linear stretch in the spatial direction. While a tight negative linear correlation is apparent for scans in the forward direction, the trace extents differ stronger for scans in the reverse direction and show a change in slope at −0.2 px. To simplify the correction of this systematic, we continued assuming a single linear correlation for each scan direction.

The method of determining the relative trace extent for each exposure is inaccurate, as spatial profiles of individual exposures can deviate in shape. In these cases, the fit is not only guided by the trace extent but also by other characteristics of the spatial profile. Therefore, we directly de-correlated the flux against the row shift derived from each first non-destructive exposure during the light curve fits in Sect. 4.1.

Fig. 2 Relative stretch of the light trace as a function of trace position on the WFC3 detector. Both values are relative to the first exposure of the observations. Values for forward and reverse scan directions are green and purple, respectively. Different symbols indicate the nine individual visits of the planet. The gray diamonds mark orbit 29, which was removed for the PACMAN analysis.

3.2 Data reduction with Eureka!

Eureka! (Bell et al. 2022) is a publicly available, community-developed pipeline designed to reduce and analyze exoplanet time-series observations with particular emphasis on data from JWST. It converts raw, uncalibrated FITS images into precise exoplanet transmission or emission spectra. It features a modular design of six stages, of which four (Stages 3–6) are employed for HST WFC3 transit observations.

Each stage in the pipeline is guided by “Eureka! Control Files” (ECFs) and “Eureka! Parameter Files” (EPFs), with the latter used to adjust transit model fit parameters. Eureka! currently provides template ECFs for JWST’s MIRI, NIRCam, NIRSpec instruments, and HST’s WFC3 instrument. Further details can be found on the Eureka! ReadTheDocs page.

To reduce HST WFC3 data with Eureka!, we used the _ima files provided by MAST as a starting point. Stage 3 of the pipeline performed background subtraction and optimal spectral extraction on calibrated image data (Bell et al. 2022; Ashtari et al. 2025). The extraction apertures and regions for background estimation were assigned individually for each visit to achieve the best possible spectrum extraction. As already demonstrated with the PACMAN reduction in Fig. 2, this was necessary due to the non-negligible shift of the spectral trace in scan-direction. This stage resulted in a time series of 1D spectra.

Stage 4 used this time series of 1D spectra to generate the light curve by binning it along the wavelength axis, then removed drift and jitter to produce spectroscopic light curves (Ashtari et al. 2025). Stages 5 and 6 of Eureka! conducted the light curve fitting and plotting of the results. Details about the light curve fitting can be found in Sect. 4.2.

Eureka! also detected the non-negligible shift of the spectral trace in the scan direction on the detector (see Fig. 2). For seven of the eight transits analyzed with Eureka!, the shift in scan direction is greater than the recommended 0.2 pixel (peak-to-peak) threshold for precise time series observations (Stevenson & Fowler 2019). Because each pixel has distinct sensitivity characteristics, spatial shifting of the spectrum introduces variability in the measured flux (Stevenson & Fowler 2019). Correcting for these effects requires precise modeling of positional shifts that are dynamic and do not follow simple patterns. As such, modeling is currently not implemented in Eureka!, the positional dependence limits the precision that can be achieved with this pipeline. The limitations and capabilities of Eureka! for fitting the light curves of this dataset is discussed further in Sects. 4.2 and 4.3.

3.3 Stellar parameter estimation with ARIADNE

For the analysis of the planet data, it is essential to obtain precise parameters of the star, especially those needed to determine the best-fitting stellar models. We derived the stellar parameters based on fits to the SED with ARIADNE (Vines & Jenkins 2022) code.

ARIADNE is an open-source Python package designed to automatically fit the SED of stars to different stellar atmosphere model grids. For HD 86226, we used the Phoenix v2 (Husser et al. 2013), BT-Settl (Hauschildt et al. 1999; Allard et al. 2012), Kurucz (1993), and Castelli & Kurucz (2003) model grids, which were previously convolved with several broadband filters. ARIADNE operates under a Bayesian framework, where each grid is modeled using a Nested Sampling (Skilling 2004, 2006; Higson et al. 2019) algorithm as implemented by dynesty (Speagle 2020). The algorithm calculates the evidence of each model, while also computing the posterior distributions for the effective temperature, log₁₀(g), [Fe/H], extinction in V, and radius of the star. After each grid was fitted, we combined the posterior distributions using Bayesian model averaging, which leverages the evidence of each model and accounts for model intrinsic uncertainties. We summarize the results from ARIADNE in Table 1.

3.4 STIS stellar spectral analysis

The STIS spectra were automatically calibrated with CALSTIS v.3.4.2 and retrieved from MAST. The automated extraction successfully identified the spectral trace in all three gratings, and no custom extractions were required. We detect multiple strong emission lines in the spectra and a clear continuum signal redward of ≈1350 Å. To splice the three spectra together, as well as the orders of the echelle E230M spectra, we retrieved the Hubble Advanced Spectral Products (HASP⁵) complete coadd of the spectrum.

The Lyman α line is by far the brightest far-UV (FUV) line, but is strongly absorbed by the local interstellar medium (LISM, see Redfield & Linsky 2008) and must be reconstructed in order to complete the FUV spectrum. Fortunately, the wings of the Lyman α line are clearly detected in the G140M spectrum, allowing us to reconstruct the intrinsic line profile with lyapy⁶ (Youngblood et al. 2016, 2022). We used the example start file provided with the package, changing the input spectrum, line spread function (to the G140M LSF), and the number of steps to 50 000. We also fixed the self-reversal parameter p = 2.4 (Taylor et al. 2024). The key parameters of the Lyman α reconstruction are given in Table 2. The LISM column density and velocity are consistent with the sample of LISM absorptions of nearby stars (Redfield & Linsky 2008) within 3 σ. We note that the large uncertainties in the derived values are caused by an air-glow line at the position where we expect the LISM absorption. This caused small uncertainties in the background subtraction to propagate into large uncertainties in the LISM constraints.

The gap between the G140L and G230M spectra covering 1715–2274 Å and wavelengths greater than 3118 Å were filled with a Phoenix model spectrum from the Lyon BT-Settl CIFIST 2011_2015 grid (Allard 2016; Baraffe et al. 2015) using the stellar parameters from Sect. 3.3. Uncertainties were estimated using models with that account for the uncertainty in T_eff. We find that the model is a good match to both sides of the gap in the UV spectrum, although the model diverges from the spectrum at wavelengths smaller than 1700 Å. This indicates that the continuum contribution from the stellar chromosphere begins to dominate over the photosphere flux at those wavelengths.

Wavelengths below 1160 Å were filled using the FUV to extreme UV (EUV) scaling relations from France et al. (2018). We used the scaling relationship for the Si IV 1400 Å doublet as it is detected with high S/N in the G140L spectrum. The line flux was measured by integrating the G140L spectrum over the region 1390–1410 Å then subtracting a flat continuum flux fitted to the surrounding spectrum. We find F_Si IV = (1.25 ± 0.18) × 10⁻¹⁵ erg s⁻¹ cm⁻². Propagating both the error on the flux and on the scaling relationships, we find log₁₀ (F(90 − 360 Å)[erg s⁻¹ cm⁻²]) = $-{13.6}_{-0.04}^{+1.00}$ and log₁₀(F(360 − 911 Å)[erg s⁻¹ cm⁻²]) = $-{13.8}_{-0.04}^{+1.00}$. We generated a two-step spectrum for the region 90–1160 Å with these relationships, extending the higher wavelength band to meet the blue end of the STIS spectrum.

Our completed SED is shown in Fig. 3 with the various components labeled. We compare the SED to that of the Sun. As might be expected from the stellar parameters, they are almost identical. The slight offsets between the spectra are within the uncertainties of the reconstructed spectrum.

The flux at wavelengths shorter than 1400 Å is strongly affected by stellar activity, in addition to the effective temperature. The close agreement of HD 86226 with the flux of the Sun therefore implies similar activity levels of both stars. This is also traced by the 25-day rotational period of HD 86226 (Arriagada 2011), which is similar to the solar rotation period at the equator.

Fig. 3 Full reconstructed SED of HD 86226. Various components ar described in the text labeled in the top legend. The solar spectrum was taken from Woods et al. (2009) and scaled to the distance of HD 86226. The Phoenix models of HD 86226 were binned to a lower resolution for this figure to allow for better comparison to the solar spectrum.

3.5 Host star optical monitoring

The ground-based photometric monitoring data were handled and reduced with a customized pipeline. The World Coordinate System (WCS) parameters of the images were determined by the Astrometric STAcking Program ASTAP⁷, using star databases made with Gaia data releases. We used Astropy’s ccdproc package (Craig et al. 2015) for bias, dark, and flat calibrations and photutils (Bradley et al. 2024) for aperture photometry. The aperture photometry process found centroids of stars in the field and placed circular apertures with a radius of 1.4 times the smallest full width at half maximum (FWHM) of stars in the field. Local backgrounds were determined with annular apertures between 2.4 and 3.6 times the FWHM and subtracted. We used field stars and their magnitudes obtained from SIMBAD⁸ as standards, with seven stars in the B band, five in the V band, and two in the R and I bands each.

The light curves shown in Fig. 4 retain the magnitude of HD 86226 as measured with single exposures and binned by observing nights. The figure also highlights the HST visit times. The photometric monitoring dataset covers most HST visits except for the last two WFC3 observations. We characterize the statistics of the photometry results in Table 3. The sigma-clipped mean magnitudes, averaging entire light curves, match well with literature values, indicating good overall photometry quality (see Cols. 2 and 5 in Table 3). The night-to-night standard deviation is comparable to the median exposure-to-exposure standard deviation within given nights (Cols. 3 and 4 in Table 3), indicating no night-to-night brightness fluctuation. The lack of significant photometric stellar activity matches our expectation that the host star is relatively quiet and its activity should not significantly impact the transmission spectrum.

Fig. 4 Ground-based photometric monitoring light curve of HD 86226. The large colored markers represent nightly mean magnitudes. The small gray markers represent individual exposures. The sigma-clipped mean magnitudes averaging the entire light curves and the ±1 σall ranges are represented with horizontal lines and shaded regions. Vertical lines correspond to HST observations; the last HST visit in January 2025 is not shown. We identify no evidence of stellar activity during the monitoring period.

4 Light curve fitting

4.1 Light curve fitting with PACMAN

We fitted the PACMAN-reduced data simultaneously with a batman (Kreidberg 2015) transit model and a systematics model, which are applied multiplicatively. Following common practice, we removed the first orbit in each visit due to its relatively strong exponential ramp. As the flux of the first exposure in each orbit is much lower than the later exposures, we also removed the first exposure of each orbit. Lastly, we removed all exposures from orbit 29, as the increased spectral trace instability hindered the row shift correction (see Sect. 2 and Fig. 2).

For the transit model, we fixed the planet’s semi-major axis, inclination, eccentricity, and argument of periastron to the values derived by Teske et al. (2020). The period of the planet was fixed to the value derived by Kokori et al. (2023). We used a quadratic limb darkening law and fixed the limb darkening coefficients to the values obtained from the Castelli & Kurucz (2003) stellar grid with the exotic-LD (Grant & Wakeford 2024) Python package. Attempts to fit the limb darkening coefficients with a linear limb darkening law and free parameters led to fit coefficients similar to the model values. In contrast, attempts to fit a quadratic law with free parameters resulted in strong degeneracies between the free parameters. Due to this, fixing the limb darkening coefficients to the model values for a quadratic law provided the most accurate solution for our fit. Therefore, free parameters for the astrophysical model were only the transit mid-time and the transit depth. For both parameters, we fitted one common value across all visits. The transit mid-time was hereby extrapolated from the first to the subsequent visits using the orbital period of the planet.

For the systematics model, we fitted a constant flux and offset between forward and reverse scan flux for each visit. In addition, the flux behavior with time for each individual visit was described by a linear slope across the orbits and an exponential ramp with the same shape for each orbit. The exponential ramp was parametrized by $$\exp (-r_1\times t_{\text{orbit }}-r_2),$$(1)

where t_orbit is the time since the start of the associated orbit and r₁ and r₂ are free parameters. In addition, a linear slope in flux with row shift (see Sect. 3.1) was fit to all visits simultaneously for forward and reverse scans, respectively. To account for possible neglected sources of uncertainty, we additionally fitted a scale factor that inflates the uncertainties to achieve a reduced χ² of one.

We ran a Marcov Chain Monte Carlo (MCMC) algorithm with the Python package emcee (Foreman-Mackey et al. 2013) to obtain the final parameter estimates and uncertainties. We used broad, uninformed linear priors for all parameters of the systematics model and the transit depth. We ran chains with 100 000 steps and 500 walkers, removing the first 4000 steps of each chain as burn-in. This ensured that the chains were longer than 50 times the autocorrelation time for astrophysical and systemic parameters and that the chains converged. The uncertainties on the parameters were derived from the posterior distributions and based on the Gaussian 1 σ deviations of the parameters from their median value. Minor correlations are seen in the posterior distributions of the constant flux offsets and the linear slopes of flux with time. Furthermore, the row shift slopes between forward and reverse scans are strongly correlated. Since both parameters describe the shift of the spectral trace on the detector, a tight correlation between these two systemic corrections is expected. The obtained transit depth is not strongly correlated with the other parameters.

The broadband light curve is shown in Fig. 5. The root mean square (RMS) noise of 45.2 ppm is above the expected photon noise of 25.2 ppm, but within the best achieved RMS values for HST data (see e.g., Kreidberg et al. 2014a; Tsiaras et al. 2016). Major contributors to the remaining red noise are a wave-like residual pattern with time and the systematically lower egress flux. The wave-like red noise can be removed by fitting a second-order polynomial with time to the flux of each visit instead of a linear slope. We tested if this model is a better fit for the systematics by calculating the Bayesian information criterion (BIC, Kass & Raftery 1995) for both models. Using a second-order polynomial improved the BIC of the broadband light curve fit significantly by 78.4, as compared to the linear slope with time. However, as seen in Table 4, this improvement only holds for five spectroscopic bins, while the first-order polynomial is preferred for most bins. To avoid overfitting the spectroscopic data, we opted to fit a linear slope with time to all datasets.

The observed egress flux is systematically lower than the model flux by approximately 50 ppm. Since this lower flux is not seen during ingress, we can exclude issues with the limb darkening as the cause of this lower egress flux. Most of this offset can be accounted for by fitting a second-order flux baseline for each visit. However, due to the many steps involved in fitting the light curve, we cannot unambiguously say if the decreased egress flux is astrophysical or just an artifact of the fitting process. A continuous time series observation of this planet could clarify if this egress flux drop persists.

From the broadband light curve, we derive a transit depth of HD 86226c of 410 ± 11 ppm in the full 1.1–1.7 μm band observed with HST/WFC3. This results in a planet radius of 2.313 ± 0.051 R_⊕, assuming the stellar radius given in Table 1. This planet radius is slightly larger than the value of 2.16 ± 0.08 R_⊕ derived by Teske et al. (2020). While the planet may appear smaller in the wavelength range covered by TESS, the radius difference of less than 2 σ does not allow for definite conclusions. For the analysis of the planet spectrum in Sect. 3, we assumed a planet radius of 2.313 ± 0.051 R_⊕, as derived from the HST data.

We explored different variations of the treatment of the row shift decorrelation parameter for the fitting of the spectroscopic light curves. We tested whether it is possible to fix the slopes of the row shift correction to the values of the broadband light curve. Except for three bins, the BIC strongly favors leaving these slopes as free parameters for the spectral light curve fits (see Table 4). Furthermore, we tested whether fitting the linear slopes with row shift per visit is useful, analogous to the remaining systematics. As this adds 16 additional parameters to the systematics models, the BIC disfavors this approach even for the broadband light curve. Removing the correction for the row shift in the systematics model increases the BIC in every wavelength bin between 28 and 543, demonstrating that it is necessary to consider this effect when fitting the data.

To obtain the spectrum of the planet, we divided the signal into 21 spectral bins between 1.12 μm and 1.65 μm. The transit mid-point was fixed to the value derived from fitting the broadband light curve, while transit depth and all systematics were free parameters. Limb darkening coefficients were calculated for each wavelength bin using the ExoTiC-LD (Grant & Wakeford 2024) package for Python. The spectrum is shown in Fig. 6 and Appendix A lists the transit depths of the individual spectral light curves. We do not see major contributions of red noise to the RMS of the spectroscopic light curve fits with PACMAN (corresponding figures are available via Zenodo). Possible atmospheric models that can explain this spectrum are discussed in Section 5.

Fig. 5 Broadband (top panel) and spectroscopic (bottom panel) light curves obtained with PACMAN by combining all nine observed transits of HD 86226 c. Differently colored data points in the broadband light curve mark the respective visits. The right panel shows the residuals from the fit transit model. Numbers in the bottom panel state the wavelength bin in μm (left panel) and the RMS in ppm (right panel).

Fig. 6 Near-infrared transmission spectrum of HD 86226 c obtained with PACMAN (green) and Eureka! (purple). The corresponding shading indicates the ±1 σ range of a constant transit depth fit to each spectrum.

4.2 Light curve fitting with Eureka!

The data reduced with Eureka! were also simultaneously fit with a batman (Kreidberg 2015) transit model and a systematics model that was applied multiplicatively. In contrast to PACMAN, Eureka! was not designed for analyzing multiple transits simultaneously, so we implemented a custom method to perform a joint fit.

This method shared the fit parameters of the astrophysical model across all visits, while allowing the parameters of the systematics model to be freely determined for each visit. With this setup, each iteration of the fitter determined the best-fit systematics model for each visit individually while keeping a shared astrophysical model. The best-fitting model was determined by the combination of systematic and astrophysical parameters that yielded a minimized deviation between data and model across all visits.

The two scan directions were fit separately to account for the offset between forward and reverse scans. We determined the final transit depth by calculating a weighted average between the fits to the two scan directions. If the flux offset between scan directions were constant, normalizing the respective light curves would lead to identical transit depths. As the spectrum shift on the detector leads to a non-constant offset between both scan directions, this averaging provides an intermediate solution between both directions without correcting for the shift.

The astrophysical model had the transit midpoint and the transit depth as the only free parameters. The remaining orbital parameters were fixed to the literature values derived by Teske et al. (2020) and Kokori et al. (2023), analogous to the PACMAN light curve fit. A quadratic limb darkening law with fixed coefficients from the Castelli & Kurucz (2003) stellar grid was used.

The systematics model contained a first-order polynomial per visit and the hstramp model implemented in Eureka!. This model fits an exponential ramp to each orbit according to: $$h_0\times \exp (-h_1\times t_{\text{batch }}),t_{\text{batch }}=(t_{\text{local }}-h_5)\mathrm{\%}{h}_4.$$(2)

In this equation, h₀ − h₅ are free parameters, and t_local is the time since the visit started. Eureka! furthermore provides the capability to fit a second-order polynomial per orbit with the parameters h₂ − h₃, which we did not use. An additional parameter, h₆, was varied as a fixed parameter with possible values of 0.2, 0.3, or 0.4. It fits the timing offset of the hstramp model.

A range of other correction techniques were tested, including individually and simultaneously fitted transits, positional adjustments (with the inbuilt functions xpos, ypos, ywidth), and separate consideration of forward and reverse scan directions, where spectra generated by only single scan directions were evaluated. While effective under more stable conditions, these methods proved challenging to implement for this dataset. Isolating data segments by removing integrations or entire transits did not fully resolve the issues.

The spectroscopic light curves were fit using the divide_white method (see e.g., Kreidberg et al. 2014a). In addition to fitting the systematics model of the respective wavelength bin, the data were further corrected by the residuals between the data and the broadband light curve model. This was necessary to remove the remaining red noise seen in the broadband light curve (see Appendix B and the corresponding figures available via Zenodo). Furthermore, the transit mid-time was fixed to the value derived from the broadband light curve fit.

Parameter estimates and uncertainties for the fits were obtained using a MCMC algorithm for the light curve fits with the python package emcee (Foreman-Mackey et al. 2013). The final light curve fits still show a large contribution of red noise in most of the 27 spectral bins between 1.12 μm and 1.66 μm. A major contributor to this systematic noise is the inclusion of the first orbit of every visit. The exponential ramp of this orbit is not identical to the ramp in the subsequent orbits, even after removing different amounts of exposures at the start of this orbit. Unfortunately, removing the first orbit of each visit in this reduction was impossible, as the many half-covered transits terminated the fitting with Eureka! in this case. The transit depth of the first wavelength bin between 1.12 μm and 1.14 μm deviates strongly from the remaining bins, which is why we remove it from further analysis.

Due to the remaining red noise, the uncertainties on the transit depths derived with Eureka! are underestimated. To account for this, we applied a β scaling similar to the methods explained by Pont et al. (2006) and Cubillos et al. (2017) to our results. The resulting spectrum with inflated uncertainties is shown in Fig. 6, and the inflation method is detailed further in Appendix B.

4.3 Differences between the spectra obtained with PACMAN and Eureka!

The spectra obtained with PACMAN and Eureka! are both featureless but offset by 30 ppm. In addition, the Eureka! data points show a larger scatter, especially at small wavelengths. Major differences can be found in the light curve fitting processes of both pipelines: We included different data, chose different systematics treatments, and corrected for the trace position with the PACMAN pipeline only.

The offset in transit depth is caused by the inclusion and exclusion of the first orbit of every visit by Eureka! and PACMAN, respectively. The ramp shape of this first orbit differs from the subsequent ramps, which introduces a systematic offset in transit depth when all orbits are treated equally. Including the first orbit of every visit in the PACMAN light curve fits also leads to an increased transit depth comparable to the average value in the Eureka! spectrum.

One cause for the increased scatter in the Eureka! reduction might be the use of the divide_white method: This method assumes that the fitted systematics model is wavelength independent. However, the spectroscopic fits with PACMAN show that this is not the case. Especially, the first-order polynomials of the baseline flux as a function of time are strongly dependent on the underlying stellar flux that changes continuously with wavelength.

Another difference between the systematics models is how we account for the upstream-downstream effect and the associated flux dependence on the spectral trace position. While PACMAN jointly fits all data with a constant offset between forward and reverse scans of a visit, Eureka! fits the scan directions separately. Including more data in a single fit might be beneficial for robustly determining the parameters of the complex underlying model. The row shift correction implemented in PACMAN significantly improved the light curve fits and reduced the scatter in the derived spectrum. However, the derived spectroscopic transit depths deviate from the uncorrected values by less than 1 σ. This correction therefore only partially accounts for the differences in the spectra.

To further investigate the reasons for the deviating spectra, we also inspected the raw light curves. We noticed that these already show an average systematic deviation of 600 ppm in the normalized flux. While the flux in some orbits only shows minor deviations of a few ppm, other orbits show a systematic flux offset across the full orbit. This is indicative of possible differences in the spectrum extraction. While both pipelines implemented an optimal extraction algorithm, the extraction regions differ. Eureka! uses a fixed, user-assigned region, while PACMAN determines the optimal extraction region individually for each exposure. Given the trace position instability in and between orbits (see Sect. 3.1), this individual treatment of exposures is important for our data. Furthermore, the background estimations differ: While Eureka! considers pixels outside of a certain row interval, PACMAN considers all pixels in the difference-image of the spectral traces below a threshold value.

These differences in the reduction also lead to different qualities of the light curve fits: While the RMS of the PACMAN spectroscopic light curve fits is typically only about 10% above the photon noise, the Eureka! RMS is about twice the photon noise in many bins. Based on these findings, we consider the PACMAN results more reliable than those from Eureka!. PACMAN was designed specifically for HST/WFC3, and this is a particularly challenging dataset. The main problems in the Eureka! analysis arose from the bright host star and corresponding poor trace position stability and from the partial transit observations. While these effects were accounted for with PACMAN, Eureka! does not currently provide the required capabilities. Given the featureless nature of the spectrum, implementing these capabilities to Eureka! is out of the scope of this analysis. For further analysis and modeling of the planetary atmosphere, we proceeded with the spectrum obtained with the PACMAN pipeline.

5 Transit spectrum models

Given the precision of our measurements, a visual inspection of the transmission spectrum of HD 86226c (see Fig. 6) indicates that it is likely featureless. As a simple test for the presence of spectral features, we performed a least-squares fit to a constant transit depth with the curve_fit function of the Python package scipy.optimize (Virtanen et al. 2020; Vugrin et al. 2007). For the PACMAN spectrum, this leads to a value of 418 ± 14 ppm, with a χ² deviation of 16.6. According to chi-squared distributions with 20 degrees of freedom, the data is compatible with a constant transit depth within 0.4 σ.

The featureless spectrum strongly rules out a clear solar-metallicity atmosphere for HD 86226 c, but other scenarios such as high atmospheric metallicity or high altitude clouds are still possible (Benneke & Seager 2013). To quantify these potential scenarios, we used the retrieval package petitRADTRANS (pRT; Mollière et al. 2019; Nasedkin et al. 2024) to fit atmospheric models to the spectrum of HD 86226 c (see Sect. 5.1). In addition, we ran self-consistent cloud and haze models for the planet, which are described in Sections 5.2 and 5.3.

5.1 Retrievals with petitRADTRANS

Retrieving a featureless spectrum does not lead to a precise constraint of physical or chemical properties of the atmosphere but rather to a broad range of conditions consistent with the observed spectrum. Due to this, we employed a simplified approach with a minimal number of retrieved parameters.

To simulate the spectra of the planet, we assumed an isothermal atmosphere with the equilibrium temperature T_eq = 1310 K (Teske et al. 2020). This approximation is sufficient, as we observe the planet in transmission and are only sensitive to a narrow layer of the atmosphere (e.g., Fortney et al. 2010). We simulated the atmosphere in a pressure range between log₁₀(P[bar]) = 2 and log₁₀(P[bar]) = −8. To calculate the transmission spectrum, pRT first converts the pressure structure of the atmosphere to a radius structure, assuming hydrostatic equilibrium. To specify this conversion, it is necessary to define a pair of reference pressure, P_ref, and reference radius, R_ref. Based on the transit depth of the broadband light curve (see Sect. 4.1), we fixed the reference radius to R_ref = 2.313 R_⊕. The reference pressure was left as a free parameter with a log-linear prior across the entire pressure range to account for offsets in the transit depth. We retrieved the planet mass with a Gaussian prior centered on the value 7.25 M_⊕ with the 1 σ width of 1.19 M_⊕, following the value obtained by Teske et al. (2020). Lastly, to calculate the transit depths, we fixed the stellar radius to R_star = 1.047 R_⊙, as derived from the ARIADNE models (see Sect. 3.3).

We also included a gray cloud with variable pressure in the retrieval. This was implemented such that the atmosphere is fully transparent at altitudes above the cloud level and fully opaque below. The clouds were placed at the layer with the pressure P_cloud, which was retrieved with a uniform prior across the simulated pressure range between log₁₀(P_cloud[bar]) = 2 and log₁₀(P_cloud[bar]) = −8. With this cloud implementation, the spectra do not contain absorption signals that would be imprinted at altitudes below the cloud deck. The amplitude of the corresponding absorption features can therefore be muted.

We used two approaches for the atmospheric chemistry. In the first approach, we assumed that the atmosphere has a scaled solar composition in chemical equilibrium. The abundances were obtained from the pRT subpackage chemistry.pre_calculated_chemistry, whose tables were calculated with the easyCHEM python package (Lei & Mollière 2024)⁹. Since the spectrum does not inform us about the abundances of individual species, we fixed the C/O ratio to 0.55. This value was estimated based on the C/H and O/H abundances of HD 86226 stated in the Hypatia catalogue (Hinkel et al. 2014) and is equal to the solar value (Asplund et al. 2009). We retrieved the metallicity of the atmosphere $$\left[\frac{\mathrm{M}}{\mathrm{H}}\right]={\log }_{10}{\left(\frac{N_{\mathrm{M}}}{N_{\mathrm{H}}}\right)}_{\text{planet }}-{\log }_{10}{\left(\frac{N_{\mathrm{M}}}{N_{\mathrm{H}}}\right)}_{\odot },$$(3)

where N_H and N_M are the number fractions of hydrogen and species heavier than helium, respectively. We let [M/H] vary from −1 to 3 and included opacities from H₂O (Polyansky et al. 2018), CO (Rothman et al. 2010), CO₂ (Yurchenko et al. 2020), and CH₄ (Hargreaves et al. 2020).

The second approach assumed a simplified atmosphere consisting only of H₂, He, and H₂O. We focused on H₂O because it has a strong feature in the observed wavelength range between 1.1 and 1.7 μm and is a species with high expected abundances based on planet formation models (Fortney et al. 2013; Bitsch et al. 2021; Izidoro et al. 2021; Burn et al. 2024). For a planet with such a high equilibrium temperature (T_eq = 1310 K), we do not expect high abundances of CH₄ as it is photodissociated (Zahnle et al. 2009; Gao et al. 2020). Other species that could be expected in the atmosphere of HD 86226 c are CO and CO₂ (e.g., Moses et al. 2013). However, these only have shallow features in the observed wavelength range, which is why we neglected them for this analysis. During the retrieval, we fitted for the H₂O mass fraction, X_H2O, and calculated the H and He mass fractions according to their solar ratio (see Asplund et al. 2009). We note that this ratio could differ for sub-Neptunes, as a preferential mass-loss of H over He could alter the atmospheric composition (Malsky et al. 2023). However, our data is not sensitive enough to determine a change in the H-He ratio.

In addition to molecular absorption, we also included collision-induced absorption of H₂ and He, as well as Rayleigh scattering from H₂, He, and H₂O as additional opacity in both approaches. We used the opacities for collision induced absorption between H₂ molecules from Borysow et al. (2001); Borysow (2002), and between H₂ and He from Borysow et al. (1988, 1989); Borysow & Frommhold (1989). Opacities from Rayleigh scattering were computed internally by pRT, based on the Rayleigh scattering cross-sections of H₂ (Dalgarno & Williams 1962), He (Chan & Dalgarno 1965), and H₂O (Harvey et al. 1998).

The best-fit models and the posterior distribution of the parameters for both chemistry treatments are shown in Fig. 7. From these, it is clear that we cannot distinguish whether high-altitude clouds or a high-metallicity atmosphere cause the featureless spectrum. High-altitude cloud models allow for most values of [M/H] since the cloud always hides the corresponding molecular features. Also, models with the highest [M/H] allow the opaque cloud deck to be placed at most pressure levels. For these, the spectrum is sufficiently flat due to the high mean molecular weight alone. At intermediate cloud pressures (log₁₀(P_cloud[bar]) = −3 to 0) and metallicities below [M/H]=2, feasible models cluster in a triangle shape, allowing for higher cloud pressures at lower metallicities. This is reasonable, as the respective molecules form absorption features deeper in the atmosphere for low metallicities. Therefore, low-altitude clouds can also sufficiently mute their spectral features.

The reference pressure and gray cloud pressure level have a nearly one-to-one correlation. Both reference and cloud pressure control the effective planet radius in the model and, consequently, the transit depth of the featureless spectrum. The parameters are therefore expected to have similar values. A deviation from this correlation is seen for high cloud pressures. For these, the transit depth is not set by the cloud but mainly by the assumed reference pressure.

Attempts to run the retrievals with a free atmospheric temperature were unreliable, as the models always ran into the lower prior boundaries of the temperature interval. Expanding the prior interval to between 50 K and 3000 K led to a preferred value of T = 200 K, much lower than the equilibrium temperature of the planet. These solutions fit the spectrum equally well as the high-temperature models with high metallicity or high-altitude clouds. An atmosphere with a lower temperature also has a lower scale height, decreasing the feature size. For the retrieval, this leads to a large prior volume with equally good fits of low-temperature models, which is then also seen in the posterior. To avoid biasing our results with unphysically low temperatures, we keep the planet temperature in the retrieval fixed at 1310 K.

We also tested our sensitivity to cloud-free atmospheres. To determine the metal enrichment that could produce a sufficiently flat spectrum, we ran the retrievals with scaled solar chemistry again without the gray cloud deck. These retrievals show that the featureless spectrum can also be reproduced with a metallicity of [M/H]=${2.9}_{-0.2}^{+0.1}$ and a reference pressure of log₁₀(P_ref[bar]) = $-{2.6}_{-0.4}^{+0.5}$. In the absence of clouds, we derive a 3 σ confidence lower limit on the atmospheric metallicity of HD 86226c of [M/H]=2.3, corresponding to a minimum mean molecular weight of approximately 6 u. We note that the precalculated grids of pRT only cover atmospheric metallicities of up to [M/H]=3, which is also the upper boundary of our prior range. Therefore, the actual metallicity of the planet may be much higher than the derived lower limit.

To further quantify the significance with which we can rule out a cloud-free solar-metallicity atmosphere, we retrieved the spectrum a final time without the gray cloud deck and a fixed solar [M/H]=0. To avoid artificially muting the feature size by an overly increased planet mass, we furthermore fixed the planet mass to 7.25 M_⊕ (Teske et al. 2020). The best-fitting spectrum with a χ² of 89.9 has a reference pressure of log₁₀(P_ref) = $-{1.50}_{-0.07}^{+0.07}$. As shown in Fig. 7, the cloud-free solar-metallicity model does not fit the observed spectrum well. Based on the derived χ² for one free parameter, this model is ruled out with a confidence of 6.5 σ. We note that including mass as a free parameter only excludes the cloud-free solar-metallicity model with confidence of 2.9 σ. However, the corresponding retrieval finds a planetary mass of 11.1 ± 0.8 M_⊕, much higher than the value of ${7.25}_{-1.12}^{+1.19}$ M_⊕ found by Teske et al. (2020).

Fig. 7 Results of the pRT retrievals for the atmosphere of HD 86226 c. Left: PACMAN spectrum (black data) and 1,2, and 3 σ percentiles of the spectra from the posteriors of the pRT retrieval with scaled solar chemistry (green shadings). The gray spectra show cloud-free models with 1× and 1000× solar metallicity. Right: posterior distribution for the two retrieval runs with scaled solar chemistry (green) and H-He-H2O atmosphere (purple). The pressures are provided in units of log10([bar]). Contours are drawn at the 1,2, and 3 σ levels. Vertical dashed lines mark the 16%, 50%, and 84% quantiles.

5.2 Cloud forward models

We also explored cloud forward models of metal-enriched atmospheres to investigate plausible cloud compositions. To do this, we used the Thermal Stability Cloud (TSC) model (Blecic et al. 2025). This model was already successfully applied to multiple synthetic and actual HST and JWST targets to distinguish various cloud species present in their atmospheres and constrain the overall cloud structure (e.g., Venot et al. 2020; Taylor et al. 2023; Bell et al. 2024). TSC is a complex Mie-scattering cloud model that is built on the approaches of Benneke (2015) and Ackerman & Marley (2001), allowing for additional flexibility in the location of the cloud base (determined by the overall atmospheric metallicity) and in the width of the log-normal distribution (which affects the droplet sizes and distribution). The model formulation employs multiple free parameters, enabling us to capture the complex structure of the clouds. These include the cloud’s shape, location, and extent; the droplet sizes, distribution, and volume mixing ratio; and the overall atmospheric metallicity. This parametrized model assumes thermally stable equilibrium atmospheric conditions (Marley et al. 2013) and does not account for non-equilibrium microphysical processes. For more details on the model parameters, refer to the model description in Venot et al. (2020), as well as Ackerman & Marley (2001) and Benneke (2015).

Before applying the TSC cloud model, we first generated 1D self-consistent cloud-free temperature profiles over a range of different metallicities, [M/H] = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 (1, 3, 10, 30, 100, 300, and 1000× solar), applying both radiative-convective and chemical equilibrium calculations (Cubillos, Blecic et al., in prep). Our radiative equilibrium approach followed the iterative procedure of Malik et al. (2017) and Heng et al. (2014) and utilized the system parameters provided in Table 1 and the stellar spectrum derived in Section 3.4 as input. The approach assumed the presence of all major molecular species formed from these most abundant elements: H, He, C, O, N, Na, K, S, Si, Fe, Ti, and V; and the following ions: e⁻, H⁻, H⁺, ${\mathrm{H}}_2^{+}$, He⁺, Na⁻, Na⁺, K⁻, K⁺, Fe⁺, Ti⁺, V⁺. We set the planetary interior temperature of HD 86226 c to T_int = 100 K and the energy redistribution factor to f = 2/3 (Hansen 2008). The model included the opacities from alkali metals, Na and K, and molecular species H₂O, CH₄, NH₃, HCN, CO, CO₂, C₂H₂, OH, SO₂, H₂S, and SiO (see Cubillos & Blecic 2021, for all relevant references and the line-lists sources). Fig. 8 shows these cloud-free temperature profiles color-coded by metallicity.

After obtaining the self-consistent temperature profiles for all atmospheric metallicities, we investigated which cloud species could condense within the HD 86226c temperature regimes. These include the species KCl, ZnS, Na₂S, MnS, SiO₂, Cr, MgSiO₃, Fe, and Mg₂SiO₄ (e.g., Gao et al. 2021; Morley et al. 2012). We excluded KCl, ZnS, Na₂S, Mg₂SiO₄, Cr, and SiO₂ from our analysis for the following reasons: KCl, ZnS, and Na₂S condense at lower temperatures and do not intersect the temperature-pressure (T − P) profiles for the metallicity range we considered, preventing cloud formation based on our thermalstability approach. Mg₂SiO₄ condenses at slightly higher temperatures than MgSiO₃. At the temperatures of HD 86226 c, MgSiO₃ condensates form first and be more abundant (Lodders & Fegley Jr 2006). We note that Cr is much less abundant than MnS or MgSiO₃, thus significantly reducing the likelihood of Cr cloud formation (Lodders 2003; Morley et al. 2012). In the presence of both Mg and Si elemental species, Visscher et al. (2010) suggest that MgSiO₃ forms before SiO₂, making SiO₂ less likely to form at all. Additionally, neither Cr nor SiO₂ has prominent features at the wavelengths of our observations (Blecic et al. 2025), making them indistinguishable from other, more likely, cloud species.

Figure 8 shows the vapor–pressure (V − P) curves of the species we included in our analysis, which could condense between temperatures of 1000–3000 K and pressures of 100–10⁻⁸ bar. It is apparent that the position of the V − P curves in the T − P space strongly depends on the atmospheric metallicity (e.g., Morley et al. 2012) shifting to lower temperatures for lower solar metallicity values and to higher temperatures for super-solar metallicities. We show this range as shaded regions for several cloud species. The plot demonstrates that the MnS, MgSiO₃, and Fe V − P curves intersect with the T − P curves for our explored metallicity and pressure range. Even at a metallicity of 1000× solar, Na₂S clouds could condense only at extremely low pressures, where their abundance is too low to produce visible signatures in the planetary spectra. Thus, we only proceeded with our grid exploration for the MnS, MgSiO₃, and Fe clouds. In Fig. 9, we also present the scattering and absorption coefficients for these species. While MnS and MgSiO₃ exhibit distinct features within the displayed wavelength range, they, along with Fe, are nearly indistinguishable in the range of our observations (1.1–1.7 μm). This lack of differentiation has posed a challenge for our analysis, as discussed below.

We generated a grid of cloud models across all metallicities, varying fundamental cloud parameters that can produce visible signatures in the spectra. In addition to metallicity, these include droplet radius, r_g, and droplet volume mixing ratio, q*, described by the number of cloud droplets over the number of H₂ molecules per unit volume. The presence of clouds effectively increases the radius of the planet and decreases the atmospheric scale height. A similar effect on the spectra could result from uncertainties in the planetary transit radius, potentially leading to either an underestimation or overestimation of the cloud droplet volume mixing ratio and droplet size. For this reason, we made the transit radius a free parameter in our model, and we also allowed for a vertical offset of the observed spectrum.

We only considered plausible values for our free parameters across our metallicity range. For the droplet radii, we chose a range from 0.01 to 100 μm, equally spaced in log space. These radii represent the median values of the underlying log-normal distribution (geometric mean droplet radii), which is smaller than the effective droplet radii reported by Ackerman & Marley (2001). We set the width of the log-normal distribution to 1.5 for all metallicities, as we empirically determined that this value best matches the data for all cloud species. We note that this choice naturally leads to a somewhat smaller droplet radius than if we fixed the width to the commonly used value of 1.8 (see Ackerman & Marley 2001). For the droplet volume mixing ratio, we chose a range from 10⁻⁷ to 10⁻³, also equally spaced in log space. Our cloud base was calculated at the intersection of the T−P and V−P curves for the corresponding metallicity. We assumed that the cloud extends to the top of the atmosphere and that the droplet volume mixing ratio is uniform across all cloud layers. For the planetary radius at the reference pressure of 1 mbar, we explored a grid of planetary radii, from 2.1 to 2.7 R_⊕, which also includes the radius of 2.313 R_⊕ derived from our broadband light curve (see Sect. 3.1). Our detailed analysis showed that this planetary radius best fits the data; therefore, in the following text, we focus on the results derived from it.

In addition to assessing how models with different planetary radii fit the data, we also allowed for a vertical offset of the observed spectra when calculating the goodness of fit (${\chi }_{\nu }^2$). This vertical offset ranged between ±11 ppm, corresponding to the 1σ uncertainty interval obtained from the broadband light curve fit (see Sect. 3.1). This few percent rescaling also accounts for uncertainty in the stellar radius. Significant instrumental offsets are not expected due to the large number of stacked transits and the consistency in transit depth across the epochs.

With each parameter set, our cloud model calculated the cloud’s scattering and absorption cross-sections and optical depth. The refractive indices of MnS and MgSiO₃ condensate species were taken from Wakeford & Sing (2015), while those for Fe were taken from Palik (1985). As a test, we also computed the MnS models with the refractive indices of Kitzmann & Heng (2018), which yielded qualitatively the same results. The scattering and absorption efficiencies were calculated using Mie-scattering theory (Toon & Ackerman 1981) by utilizing the Mie-scattering code provided by Mark Marley (private communication). This code generates pre-calculated tables of scattering and absorption efficiencies over a range of wavelengths and droplet radii, on which we then interpolated depending on our wavelength region and the cloud droplet radius (see also Fig. 9).

To generate transmission spectra including clouds, we used the PyratBay open-source framework, which is designed for spectral synthesis and retrieval of 1D atmospheric models of planetary temperatures, species concentrations, altitude profiles, and cloud coverage (Cubillos & Blecic 2021). This framework utilizes current knowledge of atmospheric physical and chemical processes, as well as the opacities from alkali metals, molecules, collision-induced opacities, Rayleigh scattering, and clouds. The tool can run in both forward and retrieval modes, using either self-consistent approaches or alternative parameterization methods to model atmospheric processes.

Although several molecules have spectral features in the WFC3/HST wavelength bands (MacDonald & Madhusudhan 2017), for this analysis, we only considered H₂O as the main and most dominant molecular source of opacity. We also assumed the presence of small aerosol particles by including both the Dalgarno & Williams (1962); Kurucz (1970) models for the H, He, and H₂ species and the Lecavelier Des Etangs et al. (2008) Rayleigh scattering model, for which we set the enhancement factor and the scattering cross-section opacity to fixed values of 0 and −4, respectively.

Figure 10 shows some of the lowest reduced chi-squared (${\chi }_{\nu }^2$) models and the grid exploration with the goodness of fit for MnS, MgSiO₃, and Fe clouds. We explored the full range of our model parameters, but only the ranges that best fit the data are shown in the tables to reduce clutter. The tables list the minimum ${\chi }_{\nu }^2$ value for each set of parameters. The color bar range is set between 1 and 2 to emphasize only the models that best fit the data.

The top panels of Fig. 10 explore the MnS cloud models. In the top left panel, we show several clear, no-cloud models and some representative cloud models with the lowest ${\chi }_{\nu }^2$. Only the flat model and high-metallicity models fit the data well (with reduced ${\chi }_{\nu }^2$ ~ 1), while lower metallicity models require some clouds to match the data. We note that the no-cloud model at 1× solar metallicity requires a high planetary radius of R_p = 2.65 R_⊕ to match the data level. This value is much larger than the value our broadband light curve analysis reports, which confirms our initial assessment that this planet either has high-altitude clouds, which increase the planetary radius, or a high metallicity atmosphere, which effectively does the same.

The middle and right panels present the ${\chi }_{\nu }^2$ values for various models across the grid. To investigate whether there are trends when exploring the grid of our model parameters, these panels depict how changing each crucial cloud parameter (while keeping the others fixed at plausible values) influences the goodness of fit. Our analysis explored all combinations of cloud parameter values across our grid, but we chose to display panels with the most ${\chi }_{\nu }^2$ values smaller than two (shown as yellow to green colors in the panels). Other parameter combinations also occasionally yielded ${\chi }_{\nu }^2$ < 2, but not as frequently as the selected examples.

The middle panel explores the goodness of fit when we keep the volume mixing ratio fixed at q* = 10⁻⁴ and change the droplet radius. The right panel shows models with a fixed cloud droplet radius and varying droplet volume mixing ratio, q*. In general, MnS models favor higher metallicity solutions, with a slight preference for a droplet radius of 0.1 μm and a droplet volume mixing ratio of q* = 10⁻⁶. However, the models are highly degenerate, particularly with respect to droplet radius, as different values of model parameters lead to similar values of the reduced ${\chi }_{\nu }^2$. This is independent of the fact that some models with a higher volume mixing ratio of 10⁻⁴ produce featureless spectra, while others with a lower droplet volume mixing ratio of 10⁻⁶ and smaller droplet radius of 0.1 μm produce significantly less muted spectra (see orange and purple models, respectively, shown in the first panel).

We performed a similar exploration for MgSiO₃ clouds (Fig. 10, middle row). MgSiO₃ clouds are located slightly deeper in the planetary atmosphere, around the 0.1 mbar level (see Fig. 8). The lowest ${\chi }_{\nu }^2$ MgSiO₃ models are generally somewhat flatter than MnS clouds. Based on our grid exploration shown in the middle and right panels, the lowest ${\chi }_{\nu }^2$ models correspond to lower metallicities than MnS clouds. They also fit the data marginally better than the MnS clouds, preferring a larger volume mixing ratio, q* = 10⁻⁴, and a slightly smaller droplet radius, r_g = 0.01 μm. We emphasize again that our q* defines the volume mixing ratio of the cloud droplets, where each droplet contains many individual cloud particles. To calculate the volume mixing ratio of the actual cloud condensates (particles), one must start with the droplet radius, account for the cloud particle radius, and consider the actual size of the droplet nucleus. For this analysis, we assumed a relatively small size for the radius of the core droplet nuclei, representing 15% of the total droplet volume. Although the actual nucleus size is not well known, if the nuclei occupied up to 50–70% of the droplet volume (Ackerman & Marley 2001), the droplet volume mixing ratio would decrease by up to an order of magnitude, from 10⁻⁴ to 10⁻⁵ for a droplet radius of 0.01 μm. We tested this hypothesis and observed only minor changes in the chi-squared values, with no change in the trend for any of the condensate species (as seen in Fig. 10). It is important to note that even with a nucleus size set to 15% of the total droplet volume, for the droplet radius of r_g = 0.01 μm, the droplet volume mixing ratio of q* = 10⁻⁴, and the MgSiO₃ particle radius of 2.5 Å, our calculated particle volume mixing ratio remains within the condensate particle volume mixing ratio reported by Morley et al. (2012), Table 2, for a solar atmospheric composition. Since most of the low ${\chi }_{\nu }^2$ models correspond to higher atmospheric metallicities, the reservoir of available MgSiO₃ particles would exceed that reported by Morley et al. (2012) for a solar composition. Overall, as for the MnS clouds, we reconfirm our conclusion that there is a high degeneracy between low and high-metallicity cloudy models and no-cloud high-metallicity models.

The Fe grid model exploration led to even more degeneracies. The spectra are flatter than for MnS clouds, but ${\chi }_{\nu }^2$ values are a bit larger than for MgSiO₃ clouds. The Fe models with ${\chi }_{\nu }^2$ values around one spread more uniformly across all metallicities, providing no clear indication of which atmospheric metallicities might be preferred but suggesting that a slightly larger droplet radius (on the order of 1–10 μm) and q* = 10⁻⁵ are needed to fit the data.

We also performed a full grid exploration assuming a planetary radius of 2.16 R_⊕ as given by Teske et al. (2020). These models, as models assuming other planetary radii, do not provide a good fit to the data. Teske et al. (2020)’s planetary radius is derived based on only optical observations and is smaller than our broadband light curve analysis radius of 2.313 R_⊕. Models with the transit radius set to 2.16 R_⊕ have the base level of the spectrum at a much lower transit depth. To compensate for this and attempt to match the data (by lifting the spectrum and increasing the transit depth), the cloud droplet volume mixing ratio increases to high, less plausible levels (on the order of 10⁻³ to 10⁻² for super-solar metallicities above 30× solar), strongly biasing our results toward a low-metallicity solution.

Overall, the results show that we cannot conclusively identify which clouds condense in the atmosphere of HD 86226 c, or whether we have a clear, high metallicity atmosphere or an atmosphere with clouds. None of the considered species (MnS, MgSiO₃, and Fe) exhibit distinct features within the wavelength range of our observations (1.1–1.7 μm), which complicates the identification of unique cloud signatures in the data. If clouds were present, their effect would be grayer and uniform across all wavelengths, which would make it impossible to identify the particular species that mutes the spectral features of HD 86226 c. Just marginally, our cloudy solutions match the data better than the high metallicity, no-cloud solutions. The implications of these findings are further discussed in Section 6.

Fig. 8 Radiative-convective temperature–pressure profiles and vapor-pressure curves. Solid lines with rainbow colors show radiative-convective temperature profiles for metallicities from 1 to 1000×solar. This figure also shows the vapor–pressure curves of the species we considered in our analysis that could condense on the temperature regimes of HD 86226 c. The vapor-pressure curves at 100× solar atmospheric metallicity are given as dashed lines, while the extent of the vapor-pressure curves over the range of 1 to 1000× solar for MnS, MgSiO3, and Na2S are given as shaded regions. The Fe range is not shown to avoid clutter, as it largely overlaps the MgSiO3 range. Na2S vapor-pressure curves do not cross our temperature-pressure profiles for any atmospheric metallicity. Thus, based on our thermal-stability approach, these clouds could not condense in the atmosphere of HD 86226 c and we did not include them in our analysis.

Fig. 9 Cloud scattering and absorption coefficients of MnS, MgSiO3, and Fe condensates. All species are featureless on the narrow wavelength range of our observations, 1.1–1.7 μm (marked with blue vertical lines), for the expected sizes of the cloud droplets of 0.01–100 μm. As shown, this is not the case for MnS and MgSiO3 clouds in the range of 0.2–15 μm, where their effect on the spectra should be distinguishable with observations that have sufficiently high resolution. We note that values are plotted in log-log space to highlight the narrow wavelength range of our observations.

Fig. 10 Comparison between different grid models and the goodness of fit for the MnS (top row), MgSiO3 (middle row), and Fe (bottom row) clouds calculated at the planetary radius set to 2.313 R⊕. For each species, the left panel shows our observations with uncertainties compared against some selected lowest reduced chi-squared (${\chi }_{\nu }^2$) models, with the corresponding ${\chi }_{\nu }^2$ values in the bottom legend. The parameter values for clear, no-cloud models are given in the left legend, and for cloudy models, they are given in the right legend. The flat model is drawn as a dotted line at a transit depth of 418 ppm. The tables in the middle and right panels display the models’ ${\chi }_{\nu }^2$ values on a grid of cloud droplet radius rg, droplet volume mixing ratio q*, and metallicity [M/H]. The table titles list fixed parameters. The color bar extent is set between 1 and 2 to highlight only the models that best fit the data.

Fig. 11 Comparison of the haze models to the HST spectrum of HD 86226 c. (Left) Goodness of fit for the parameter space (metallicity and haze formation efficiency) explored using the haze models. (Right) Spectrum models for hazy (maximum 100% haze formation efficiency; solid lines) and clear (dotted lines) atmospheres across the explored metallicity grid (1, 3, 10, 30, 100, 300, and 1000× solar). We note that in most metallicity cases, the solid and dotted lines overlap due to the negligible effect of haze.

5.3 Haze models

We also investigated the potential formation of haze in the atmosphere of HD 86226c by simulating transmission spectra that account for the presence of haze. The size and number density distributions of haze particles were computed using the photochemical and haze microphysical models of Kawashima & Ikoma (2018), following the approach of Kawashima & Ikoma (2019). First, we ran the photochemical model to determine the abundances of gaseous species. In this calculation, we adopted the UV spectrum of the host star that was observed and reconstructed in Sect. 3.4 (see Fig. 3). Next, we used the haze microphysical model to simulate the size and number density distributions of haze particles. We assumed that a certain fraction of the photodissociation of precursor molecules − CH₄, HCN, and C₂H₂ − produces haze monomers. We refer to this fraction as the haze formation efficiency. Finally, we computed the transmission spectrum models using the derived distributions of gaseous species and haze particles. For the optical properties of haze, we adopted the refractive index of tholin from Khare et al. (1984). We also tested the refractive index of soot from (Hess et al. 1998) and found that the difference in the resulting spectra is negligible compared to our observational uncertainties. We considered the same range of metallicities as in the cloud models, namely, 1, 3, 10, 30, 100, 300, and 1000× solar, and used the same temperature-pressure profile for each case. Other parameters, such as the eddy diffusion coefficient (K_zz), monomer radius (s₁), and internal density (ρ_p), were set to the same values as those used in Sect. 5.2 of Kreidberg et al. (2022); K_zz = 10⁷ cm² s⁻¹, s₁ = 1 nm, and ρ_p = 1.0 g cm⁻³.

The left panel of Fig. 11 shows the reduced χ² values for the parameter space we explored, namely metallicity and haze formation efficiency. The degrees of freedom are given by the number of data points minus the number of free parameters (metallicity and haze formation efficiency), resulting in 21 − 2 = 19. For almost all metallicity cases, the photodissociation rates of the precursor molecules are so low that the effect of haze on the transmission spectrum is negligible. This is true even at the maximum efficiency of 100% (= 10⁰), as shown in the right panel of Fig. 11. This can be explained by the relatively high atmospheric temperatures of HD 86226 c, where CO is stable instead of CH₄. The low abundance of CH₄ also suppresses the formation of the other two haze precursors, HCN and C₂H₂, via its photodissociation, leading to negligible photodissociation rates of haze precursor molecules. We note that the somewhat low χ² value for the 1× solar metallicity and 100% haze formation efficiency case is due to the relatively high abundance of CH₄, which results from the low temperature of the 1× solar atmosphere (see Fig. 8). We conclude the haze formation is suppressed across the entire metallicity range. When hazes are the only considered aerosol, high-metallicity cases are favored due to their high mean molecular weight, regardless of the haze formation efficiency.

Fig. 12 Trends in the 1.4 μm feature amplitude of gaseous exoplanets. The featureless spectrum of HD 86226 c (green marker) does not follow the trends seen for other planets. The purple parabolic trend was found by Brande et al. (2024) and describes Neptune to sub-Neptune-sized objects (black markers). Non-filled black markers show GJ 1214 b and the super-puffs Kepler-51 b and d, that were excluded from this trend. Gray markers show the feature amplitudes of gaseous planets calculated by Edwards et al. (2023). Here we only include their results for planets with equilibrium temperatures above 900 K, which mostly correspond to gas giants. The red region shows the feature size predictions of Gao et al. (2020) for planets with gravities of 10 m s−2 at pressures of 1 bar. This shaded region includes atmospheres with metallicities between 1× (dotted edge) and 10× (dashed edge) solar.

6 Discussion

The featureless spectrum of HD 86226 c continues a series of many muted sub-Neptune spectra observed in the past years (i.e., Kreidberg et al. 2014a; Knutson et al. 2014; Kreidberg et al. 2022; Guo et al. 2020; Wallack et al. 2024). This is somewhat surprising, as observation and theory predict clearer atmospheres on sub-Neptunes with equilibrium temperatures as high as the 1310 K of HD 86226c (see Fig. 1).

Previously characterized Neptune- to sub-Neptune sized objects show a parabolic trend in the observed 1.4 μm transmission feature size with planet equilibrium temperature (see e.g., Yu et al. 2021; Brande et al. 2024): Large features are observed at equilibrium temperatures below 400 K and above 700 K. In contrast, planets with intermediate temperatures show muted spectral features. This trend is possibly shaped by hazes (Crossfield & Kreidberg 2017; Yu et al. 2021) and clouds (Brande et al. 2024) in the planet atmospheres, which form most efficiently at intermediate temperatures around 600 K. However, the trend is based on planets with equilibrium temperatures below 1000 K. The featureless spectrum of HD 86226 c demonstrates that we cannot reliably extend the quadratic trend to equilibrium temperatures of 1310 K (see Fig. 12). Compared with the planets analyzed by Brande et al. (2024), the spectral feature amplitude of ${0.01}_{-0.01}^{+0.17}$ scale heights of HD 86226 c is among the smallest in the sample.

When giant planets are included, it is possible to extend the probed sample to equilibrium temperatures of 2500 K. This is done by Fu et al. (2017) and Edwards et al. (2023), who find an overall increasing feature size with equilibrium temperature. They also see a local maximum in feature size at approximately 1400 K (see Fig. 12), aligned with the featureless spectrum of HD 86226 c. This discrepancy might be caused by the differences between giant planet and sub-Neptune atmospheres. While giant planets in general have approximately solar-metallicity atmospheres, recent JWST observations suggest that many sub-Neptunes have metal-rich atmospheres (e.g., GJ 1214 b, GJ 9827 d, TOI-270 d, and GJ 3090 b; Kempton et al. 2023; Roy et al. 2023; Piaulet-Ghorayeb et al. 2024; Benneke et al. 2024; Ahrer et al. 2025). High-metallicity atmospheres have a smaller scale height, which decreases the spectral feature amplitude. In addition, their reservoir of potentially cloud-forming metals is larger. This may allow for an increased formation of clouds on sub-Neptunes, further muting their spectral features. For giant planets, Gao et al. (2020) explain the local feature size maximum with a sinking of silicate clouds to higher atmospheric pressures as the equilibrium temperature of a planet decreases. In sub-Neptunes, different atmospheric metallicities and temperature structures may prevent the clouds from sinking.

6.1 The hazy scenario

From a theoretical perspective, the formation of organic hazes is hindered at the high temperatures of HD 86226 c (see e.g., Zahnle et al. 2009; Morley et al. 2015; Gao et al. 2020), eliminating one possible cause for muted spectral features. Our haze models in Sect. 5.3 agree with these findings. For HD 86226 c specifically, haze formation from the photodissociation of CH₄, HCN, and C₂H₂ is suppressed for all considered atmospheric compositions between 1 and 1000× solar.

For atmospheric compositions with metallicities above 100× solar, He et al. (2020) demonstrate that CO and CO₂ can also function as haze precursors. If these species also contribute to lower-metallicity atmospheres, the haze production rate could be much higher than our calculations predict. Specifically, for a 1× solar atmosphere with 100% haze formation efficiency, the haze production rate increases by four orders of magnitude from 1.0 × 10⁻¹⁴ g cm⁻² s⁻¹ to 1.3 × 10⁻¹⁰ g cm⁻² s⁻¹. However, it remains unclear whether CO and CO₂ also act as haze precursors in low-metallicity atmospheres.

Since CO and CO₂ are possible haze precursors, our haze models based on CH₄, HCN, and C₂H₂ serve as a calculation with the lower limit of the haze production rate. In this lower limit, the haze produced is not sufficient to mute the atmospheric features of HD 86226 c. If aerosols shape the muted spectrum of HD 86226 c, these are either hazes primarily formed through the photodissociation of CO and CO₂ or cloud species.

6.2 The cloudy scenario

Many cloud species may contribute to opaque high-temperature atmospheres, as shown in Fig. 8. We particularly focus on MnS, MgSiO₃, and Fe clouds, as the corresponding condensation curves indicate that these species could condense at the atmospheric pressures that are probed by our observations, 10⁻²−10⁻⁴ bar. The presence of other silicate species such as Mg₂SiO₄ and SiO₂ (see e.g., Grant et al. 2023; Dyrek et al. 2024; Inglis et al. 2024) is also plausible for this planet, although these would condense at higher atmospheric pressures than MgSiO₃ and are likely to be less abundant (Lodders & Fegley Jr 2006; Visscher et al. 2010). Most cloud species typically predicted on smaller Neptune-sized planets, such as Na₂S, KCl, and ZnS, evaporate at the temperatures of HD 86226 c.

Our cloud models show that any of the considered species can mute the planetary spectrum of HD 86226 c. High-metallicity atmospheric solutions without clouds fit the data equally well as cloudy solutions (Section 6.3). To reproduce the featureless spectrum of HD 86226 c with low atmospheric metallicity and MgSiO₃ clouds, our models require droplet volume mixing ratios of 10⁻⁴. This relatively high value results partly from our sparse grid (which is not sensitive to solutions between 10⁻⁴ and 10⁻⁵) and partly from our choice of droplet nucleus size. To verify the robustness of our results, we conducted tests using a more finely spaced grid and larger nucleus sizes. The trend observed in Fig. 10 remained unchanged. Both cloudy and high-metallicity solutions stayed equally probable.

Another point to consider is that not all of these particles may be available for cloud droplet formation. Some may participate in other chemical processes or be inhibited by rapid vertical mixing. Vertical mixing can also replenish certain species, thereby sustaining the gas-phase reservoir available for condensation. In this case, the necessary species may be supplied from deeper atmospheric layers. In contrast, larger droplets may also sink to deeper layers as a consequence of gravitational settling. Our cloud prescription applies a simplifying assumption that clouds are homogenously distributed above altitudes where the condensation curves cross the pressure-temperature profile. In the presence of large-scale zonal flows such as the ones present in hot Jupiters (see e.g., Fortney et al. 2021), this description may be inaccurate. These flows may mix condensates below the pressure level predicted by this crossing. Furthermore, they would also preferably transport a specific size range of condensate droplets. Larger droplets would sink to lower atmospheric layers, while smaller droplets could be lofted. Depending on the position and extent of such a flow in the planet’s atmosphere, a vertically homogenous cloud droplet size distribution does not reproduce this effect.

Overall, the observed featureless spectrum of HD 86226c, along with its associated uncertainties, does not allow us to accurately constrain the cloud parameters or break the degeneracy between the cloudy and high metallicity solutions, and differentiate the actual condensate species muting the spectrum. However, these could be achieved utilizing the JWST, particularly the observations made with the MIRI instrument. Unlike Fe and MnS condensates, MgSiO₃ exhibits a prominent feature between 8 and 12 microns (see Fig. 9), which may be observable.

6.3 The metal-rich scenario

Recent JWST observations of the sub-Neptunes suggest that metal-rich atmospheres are common on sub-Neptunes (e.g., GJ 1214 b, GJ 9827 d, TOI-270 d, and GJ 3090 b; Kempton et al. 2023; Roy et al. 2023; Piaulet-Ghorayeb et al. 2024; Benneke et al. 2024; Ahrer et al. 2025). The featureless spectrum of HD 86226 c can also be explained by a high mean molecular weight atmosphere alone, without the need for aerosols.

As derived in Sect. 5.1, a cloudless atmosphere could reproduce the spectrum with a metallicity of [M/H]=${2.9}_{-0.2}^{+0.1}$. This results in a 3 σ confidence lower limit on the atmospheric metallicity of HD 86226c of [M/H]=2.3, corresponding to a minimum mean molecular weight of approximately 6. Assuming H₂O is the only metal in the atmosphere, this would require a water fraction of 25% in the planet’s atmosphere. Many scenarios could lead to such a high metallicity atmosphere. The atmosphere could have been enriched in volatiles by an outgassing from the planet’s interior (Piette et al. 2023). In this scenario, the dayside of HD 86226c is assumed to be hot enough to evaporate the surface of the planet. In contrast, HD 86226 c could also have formed volatile-rich (see e.g., Izidoro et al. 2021; Bitsch et al. 2021; Burn et al. 2024). Such a volatile-rich atmosphere would then be mostly stable against mass loss (Lopez 2017), explaining how the atmosphere of HD 86226 c survived until today.

A high mean molecular weight atmosphere of HD 86226 c in combination with the presence of the outer giant companion HD 86226b (Teske et al. 2020) would imply interesting scenarios for its formation history. The high mean molecular weight atmosphere could originate from the accretion of water vapor from an enriched inner disk (Bitsch et al. 2021), which might have formed from inward drifting and evaporating water pebbles (Schneider & Bitsch 2021). However, the outer growing giant beyond the water ice line could open a gap in the disk to block the inward-flowing pebbles. This would reduce the vapor enrichment in the inner regions of the disk and thus reduce the mean molecular weight of a planet forming in it. However, a high molecular weight atmosphere of HD 86226 c would indicate two possible scenarios: The giant planet could have formed interior to the water ice line of the system, where the water ice pebbles have already evaporated. Alternatively, HD 86226 c could have formed exterior to the ice line and migrated inward, before the outer giant planet could block the inward flow of pebbles (e.g., Bitsch et al. 2021). This system is the second system after HAT-P-11, where such a formation analysis is possible (see Chatziastros et al. 2024). Knowing the composition of the inner planets can therefore help us learn about the different formation scenarios of outer giant planets, a puzzle also present in the Solar System (e.g., Kruijer et al. 2017).

6.4 Possible interior structures of HD 86226 c

Previous works reported a planetary radius of R = 2.16 ± 0.08 R_⊕ for HD 86226 c (Teske et al. 2020). Our transmission spectrum updates the radius to R = 2.313 ± 0.051 R_⊕, implying a lower bulk density. Given this density estimate and the cloud-metallicity degeneracy, HD 86226 c could have either an H/He-dominated envelope with high-altitude condensates or a metal-rich atmosphere. In the first scenario, H/He would constitute no more than 1% in mass. H/He could form an atmospheric layer that is distinct from a rocky core (i.e., a gas dwarf interior structure; see Figure 1 in Lopez & Fortney 2014). However, this scenario is unlikely, as planetary cores and hydrogen envelopes are chemically reactive, producing H₂O, CH₄, and other gaseous species that mix within the envelope and core (Schlichting & Young 2022). For high-metallicity envelopes ([M/H] > 3), HD 86226 c could accommodate an envelope mass fraction of 20–50% (see Figure 5 in Aguichine et al. 2021).

The combination of mass, radius, and atmospheric metallicity alone presents a degenerate problem in constraining the size and the composition of the deep core (Kramm et al. 2011; Otegi et al. 2020). Therefore, the core composition could range from pure rock to a mixture of rock and iron with dissolved volatiles (Vazan et al. 2022; Luo et al. 2024). Determining the atmospheric abundances of H₂O, CO₂, and CH₄ with future JWST transmission spectra may help constrain the composition of the deep interior, particularly the water-to-hydrogen ratio (Yang & Hu 2024). However, further modeling is required to confirm whether this is also applicable to hot (T_eq > 1000 K) sub-Neptunes similar to HD 86226 c.

Fig. 13 Simulated JWST spectra for different atmospheric compositions of HD 86226c. The colored data points show JWST NIRSpec/G395H mock data of three transits, simulated with PandExo (Batalha et al. 2017). Colored contours show 1 σ percentiles drawn from the posterior distributions of the pRT retrievals on the simulated data. Simulated spectra are based on three edge-case atmospheric scenarios that are feasible based on the pRT fits to the HST data in Sect. 5.1: high-metallicity (M/H] > 2.5) solutions with a cloud-deck below −1 bar (green), high-metallicity solutions with a cloud-deck above −3 bar (purple), and low-metallicity ([M/H] < 1) solutions with a cloud-deck above −3 bar (gray, data points removed to avoid clutter).

6.5 How to solve the cloud-metallicity degeneracy

The only atmospheric scenario that can be unequivocally ruled out with the HST spectrum of HD 86226 c is a cloud-free, solar metallicity atmosphere. Based on our atmospheric retrievals in Sect. 5.1, we can exclude this model with a significance of 6.5 σ. Nevertheless, we cannot conclusively identify if clouds condense in the atmosphere of HD 86226 c, or whether it has a clear, high-metallicity atmosphere. The atmosphere of HD 86226 c could also combine high-altitude clouds and an enhanced metallicity.

This degeneracy could be broken by observing HD 86226 c with the NIRSpec/G395H instrument of JWST, as CO₂ has a prominent feature above 3 μm that is highly sensitive to the atmospheric metallicity. Figure 13 displays the predicted spectra for different possible atmospheric edge-cases on HD 86226 c, as derived from the HST data: high-metallicity ([M/H] > 2.5) atmospheres with a gray cloud-deck below −1 bar or above −3 bar, and a low-metallicity ([M/H] < 1) atmosphere with a cloud-deck above −3 bar. For each of these scenarios, we used pRT to generate spectra from all retrieval posteriors (see Sect. 5.1) that match the respective description. For these spectra, we included the opacities of H₂O, CO, and CO₂, as these are the dominant species predicted for planets similar to HD 86226 c by Moses et al. (2013). We used the median of these spectra to obtain a model spectrum for each scenario, respectively. For each of these three model spectra, we used pandexo (Batalha et al. 2017) to simulate NIRSpec/G395H data of three transit observations. On the simulated data, we additionally added a 1σ random noise.

To verify that the three scenarios can be distinguished, we retrieved the simulated spectrum with pRT. Assuming scaled solar abundances, high-metallicity solutions show a detectable CO₂ spectral feature. This species does not contribute to the observable opacity of a low-metallicity atmosphere, which allows us to differentiate between the scenarios. While we only included CO, CO₂, and H₂O in our predictions, OCS could be another opacity contributor that traces metal-rich atmospheres Moses et al. (2013). Depending on the exact composition of the atmosphere, the amplitude of the 4.8 μm OCS feature could be similar to that of CO₂. The cloud opacity usually diminishes in the wavelength region between 3 to 5 μm (e.g., Benneke et al. 2024; Beatty et al. 2024; Schlawin et al. 2024), facilitating a measurement of the atmospheric metallicity.

Observations with JWST/MIRI may enable a better distinction of the cloud species that mute the spectra of HD 86226 c, by probing the clouds directly. If silicate clouds are present in HD 86226 c, it may be possible to observe the silicate features around 8–10 μm (see Fig. 9 and Marley & Robinson 2015). However, directly detecting silicate features in the possibly metal-enriched atmosphere of HD 86226 c may require much observing time, as previously simulated by Piette et al. (2023). At shorter wavelengths, Rayleigh scattering on small cloud particles might produce a spectral slope. High sensitivity observations at visual wavelengths could therefore provide a further probe of the cloud coverage on HD 86226 c. Finally, another possibility for detecting and characterizing clouds on the planet is to measure the occultation depth, or phase curve, of the planet (see e.g., Kempton et al. 2023; Coulombe et al. 2025). In the presence of silicate clouds, the planet is expected to have a high albedo, reflecting a large fraction of the incoming radiation (e.g., Garcia Munoz & Isaak 2015).

7 Summary

To characterize the atmosphere of the hot sub-Neptune HD 86226c, we observed nine transits of the planet with HST/WFC3 over the wavelength range 1.1–1.7 μm. We also measured the UV spectrum of the host star with HST/STIS and reconstructed the panchromatic stellar spectrum with scaling relationships and stellar models. In addition, we performed photometric monitoring with the 0.6 m telescope at the Van Vleck Observatory.

For the HST/WFC3 data reduction, we used the open-source PACMAN and Eureka! packages. In contrast to typical WFC3 spatial scanning observations, the high scan rate for this bright star resulted in shifts in the position of the spectral trace on the detector. We observed shifts of up to two pixels, resulting in a flux measurement that was correlated with the position of the trace. To account for this effect, we implemented a function in PACMAN that decorrelates flux and trace position for both forward and reverse scan directions. For the data reduction with Eureka!, we modified the pipeline to allow for the simultaneous fitting of multiple HST transits. We find that PACMAN and Eureka! result in qualitatively similar spectra, but we note that Eureka! does not correct the trace position instabilities, resulting in a higher level of correlated noise in the light curves.

The optical photometric monitoring revealed HD 86226 as a quiet host star with a standard deviation in nightly variability of 31 mmag in the V band. Our nine transit observations do not show evidence of stellar flares or spot crossings of the planet.

The transmission spectrum of HD 86266c is featureless and consistent within 0.4 σ with a constant transit depth of 418 ± 14 ppm. This corresponds to an amplitude of just 0.01 scale heights for an H/He-dominated atmosphere. Based on an atmospheric retrieval analysis with petitRADTRANS, we ruled out a cloud-free solar metallicity atmosphere with a confidence of 6.5 σ. However, the existing data cannot unambiguously distinguish whether aerosols or a high-metallicity atmosphere are responsible for the muted spectral features. In the absence of clouds, the observed spectrum could be produced by a metal-enriched atmosphere with [M/H] > 2.3 (3 σ confidence lower limit). In addition to the atmospheric retrieval analysis, we also compared the spectrum to predictions from forward models of clouds and hazes. In contrast to most sub-Neptunes, we find that the high temperature of HD 86226 c prohibits the formation of organic hazes with CH₄, HCN, and C₂H₂ as haze precursors. On the other hand, it allows for the formation of refractory cloud species such as MnS, MgSiO₃, and Fe clouds, which are similar to those expected for hot Jupiters. All of these cloud species could plausibly match the observations.

The featureless spectrum of the hot and small sub-Neptune HD 86226 c is surprising in the context of previously observed trends for sub-Neptune atmospheres, namely that hotter planets should have larger amplitude spectral features (e.g., Brande et al. 2024). Follow-up observations with JWST may clarify the physical processes and atmospheric chemistry that shape the observed spectrum. Determining this planet’s exact atmospheric composition will show if it hosts a cloud species atypical for sub-Neptunes or if it was formed with a metal-enriched atmosphere.

Data availability

The HST data used in this paper are associated with the program GO 17192 (P.I., L. Kreidberg) and are publicly available via the Mikulski Archive for Space Telescope at https://mast.stsci.edu. Additional data products derived in this work are available via Zenodo at https://zenodo.org/records/16035649.

Acknowledgements

This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program 17192. We wish to thank Peter Gao for sharing the aerosol model spectra with us. Furthermore, we thank the anonymous referee for the detailed report that led to the investigation of many important details and a better overall clarity of the text. Q.C.T. and S.R. wish to thank Roy Kilgard for his assistance in the operation of the Wesleyan 24-inch telescope and the data storage system. K.A.K. gratefully acknowledges support from the DLR via project P.S.ASTR1508. Y.K. acknowledges support from JSPS KAKENHI Grant Numbers 21K13984, 22H05150, and 23H01224. J.S.J. gratefully acknowledges support by FONDECYT grant 1240738 and from the ANID BASAL project FB210003. P.E.C. acknowledges financial support by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship, program J4595-N. Portions of this research were carried out on the High Performance Computing resources at New York University Abu Dhabi. T.D. acknowledges support from the McDonnell Center for the Space Sciences at Washington University in St. Louis.

References

Appendix A PACMAN light curve fits

Table A.1 lists the transit depths derived in Sect. 4 from the light curve fits with pacman. Figures showing the behavior of the RMS when the individual HST exposures are time-sorted and binned (see e.g., Cubillos et al. 2017) are available via Zenodo. The broadband light curve shows residual red noise, which could be removed by fitting a second-order baseline in time to the flux of each visit. As argued in Sect. 3.1, this was not done to avoid overfitting the individual wavelength bins of the spectroscopic light curves and to have a common systematics model for broadband and spectroscopic light curves. To estimate the uncertainty of the broadband transit depth, we applied an uncertainty scaling with a constant β-factor of 2, as described by Pont et al. (2006) and Cubillos et al. (2017). This factor was derived from the fraction between fit RMS and expected white-noise RMS (see e.g., Cubillos et al. 2017) at a bin size of 25 exposures. This bin size corresponds to the duration of one HST orbit, the longest duration of continuous observations in our dataset.

For the spectroscopic bins, the increased RMS makes the deviation from a linear flux baseline undetectable, which is why the RMS behaves similar to white noise in most wavelength bins. Minor deviations from this are seen in the bin at 1.31 μm, as there is an artifact on the WFC3 detector. Furthermore, bins at the long-wavelength edge show slight signs of red noise. Since this contribution is minor, we did not account for it further.

Appendix B Eureka! light curve fits

Table B.1 lists the transit depths derived in Sect. 4 from the light curve fits with Eureka!. Plots of the corresponding light curves are available via Zenodo. With Eureka!, it was not possible to accurately describe the exact HST ramp behavior. This issue was accounted for in the fits to the spectroscopic bins by using the divide-white technique (e.g., Kreidberg et al. 2014a).

The spectroscopic light curve fits with Eureka! show red noise in most wavelength bins and the broadband light curve (the corresponding figures are available at Zenodo). As explained further in Sects. 3.2 and 4.3, this red noise mainly emerges from the use of the first orbit of every visit and the poor stability of the trace position on the detector during the observations. We accounted for this remaining red noise according to the method described by Pont et al. (2006) and Cubillos et al. (2017): The uncertainty on the transit depth of each wavelength bin was inflated by a bin-specific β-factor. This factor was derived from the fraction between fit RMS and expected white noise RMS (see e.g., Cubillos et al. 2017) at a bin size of 25 exposures. This bin size corresponds to the duration of one HST orbit, the longest duration of continuous observations in our dataset.

All Tables

All Figures

Fig. 1 Sub-Neptune targets of the SPACE program. HD 86226c is shown with a green square marker, and the other SPACE targets are shown in purple. Top: temperature-radius plane. Gray circles show the known population of planets with radii between 1.8 R⊕ and 4 R⊕ based on the entries of the NASA Exoplanet Archive3 in May 2024. Black crosses mark JWST Cycle 1–4 targets, except for the SPACE targets TOI-431 d and HD 86226 c, which will be observed in Cycle 4. Bottom: mass-radius plane. Color curves show models from Zeng et al. (2019) for various planetary compositions and temperatures.

In the text

	Fig. 2 Relative stretch of the light trace as a function of trace position on the WFC3 detector. Both values are relative to the first exposure of the observations. Values for forward and reverse scan directions are green and purple, respectively. Different symbols indicate the nine individual visits of the planet. The gray diamonds mark orbit 29, which was removed for the PACMAN analysis.
In the text

	Fig. 3 Full reconstructed SED of HD 86226. Various components ar described in the text labeled in the top legend. The solar spectrum was taken from Woods et al. (2009) and scaled to the distance of HD 86226. The Phoenix models of HD 86226 were binned to a lower resolution for this figure to allow for better comparison to the solar spectrum.
In the text

Fig. 4 Ground-based photometric monitoring light curve of HD 86226. The large colored markers represent nightly mean magnitudes. The small gray markers represent individual exposures. The sigma-clipped mean magnitudes averaging the entire light curves and the ±1 σall ranges are represented with horizontal lines and shaded regions. Vertical lines correspond to HST observations; the last HST visit in January 2025 is not shown. We identify no evidence of stellar activity during the monitoring period.

In the text

Fig. 5 Broadband (top panel) and spectroscopic (bottom panel) light curves obtained with PACMAN by combining all nine observed transits of HD 86226 c. Differently colored data points in the broadband light curve mark the respective visits. The right panel shows the residuals from the fit transit model. Numbers in the bottom panel state the wavelength bin in μm (left panel) and the RMS in ppm (right panel).

In the text

	Fig. 6 Near-infrared transmission spectrum of HD 86226 c obtained with PACMAN (green) and Eureka! (purple). The corresponding shading indicates the ±1 σ range of a constant transit depth fit to each spectrum.
In the text

Fig. 7 Results of the pRT retrievals for the atmosphere of HD 86226 c. Left: PACMAN spectrum (black data) and 1,2, and 3 σ percentiles of the spectra from the posteriors of the pRT retrieval with scaled solar chemistry (green shadings). The gray spectra show cloud-free models with 1× and 1000× solar metallicity. Right: posterior distribution for the two retrieval runs with scaled solar chemistry (green) and H-He-H2O atmosphere (purple). The pressures are provided in units of log10([bar]). Contours are drawn at the 1,2, and 3 σ levels. Vertical dashed lines mark the 16%, 50%, and 84% quantiles.

In the text

Fig. 8 Radiative-convective temperature–pressure profiles and vapor-pressure curves. Solid lines with rainbow colors show radiative-convective temperature profiles for metallicities from 1 to 1000×solar. This figure also shows the vapor–pressure curves of the species we considered in our analysis that could condense on the temperature regimes of HD 86226 c. The vapor-pressure curves at 100× solar atmospheric metallicity are given as dashed lines, while the extent of the vapor-pressure curves over the range of 1 to 1000× solar for MnS, MgSiO3, and Na2S are given as shaded regions. The Fe range is not shown to avoid clutter, as it largely overlaps the MgSiO3 range. Na2S vapor-pressure curves do not cross our temperature-pressure profiles for any atmospheric metallicity. Thus, based on our thermal-stability approach, these clouds could not condense in the atmosphere of HD 86226 c and we did not include them in our analysis.

In the text

Fig. 9 Cloud scattering and absorption coefficients of MnS, MgSiO3, and Fe condensates. All species are featureless on the narrow wavelength range of our observations, 1.1–1.7 μm (marked with blue vertical lines), for the expected sizes of the cloud droplets of 0.01–100 μm. As shown, this is not the case for MnS and MgSiO3 clouds in the range of 0.2–15 μm, where their effect on the spectra should be distinguishable with observations that have sufficiently high resolution. We note that values are plotted in log-log space to highlight the narrow wavelength range of our observations.

In the text

Fig. 10 Comparison between different grid models and the goodness of fit for the MnS (top row), MgSiO3 (middle row), and Fe (bottom row) clouds calculated at the planetary radius set to 2.313 R⊕. For each species, the left panel shows our observations with uncertainties compared against some selected lowest reduced chi-squared (${\chi }_{\nu }^2$) models, with the corresponding ${\chi }_{\nu }^2$ values in the bottom legend. The parameter values for clear, no-cloud models are given in the left legend, and for cloudy models, they are given in the right legend. The flat model is drawn as a dotted line at a transit depth of 418 ppm. The tables in the middle and right panels display the models’ ${\chi }_{\nu }^2$ values on a grid of cloud droplet radius rg, droplet volume mixing ratio q*, and metallicity [M/H]. The table titles list fixed parameters. The color bar extent is set between 1 and 2 to highlight only the models that best fit the data.

In the text

Fig. 11 Comparison of the haze models to the HST spectrum of HD 86226 c. (Left) Goodness of fit for the parameter space (metallicity and haze formation efficiency) explored using the haze models. (Right) Spectrum models for hazy (maximum 100% haze formation efficiency; solid lines) and clear (dotted lines) atmospheres across the explored metallicity grid (1, 3, 10, 30, 100, 300, and 1000× solar). We note that in most metallicity cases, the solid and dotted lines overlap due to the negligible effect of haze.

In the text

Fig. 12 Trends in the 1.4 μm feature amplitude of gaseous exoplanets. The featureless spectrum of HD 86226 c (green marker) does not follow the trends seen for other planets. The purple parabolic trend was found by Brande et al. (2024) and describes Neptune to sub-Neptune-sized objects (black markers). Non-filled black markers show GJ 1214 b and the super-puffs Kepler-51 b and d, that were excluded from this trend. Gray markers show the feature amplitudes of gaseous planets calculated by Edwards et al. (2023). Here we only include their results for planets with equilibrium temperatures above 900 K, which mostly correspond to gas giants. The red region shows the feature size predictions of Gao et al. (2020) for planets with gravities of 10 m s−2 at pressures of 1 bar. This shaded region includes atmospheres with metallicities between 1× (dotted edge) and 10× (dashed edge) solar.

In the text

Fig. 13 Simulated JWST spectra for different atmospheric compositions of HD 86226c. The colored data points show JWST NIRSpec/G395H mock data of three transits, simulated with PandExo (Batalha et al. 2017). Colored contours show 1 σ percentiles drawn from the posterior distributions of the pRT retrievals on the simulated data. Simulated spectra are based on three edge-case atmospheric scenarios that are feasible based on the pRT fits to the HST data in Sect. 5.1: high-metallicity (M/H] > 2.5) solutions with a cloud-deck below −1 bar (green), high-metallicity solutions with a cloud-deck above −3 bar (purple), and low-metallicity ([M/H] < 1) solutions with a cloud-deck above −3 bar (gray, data points removed to avoid clutter).

In the text