Near-infrared transmission spectroscopy of HAT-P-18 b with NIRISS: Disentangling planetary and stellar features in the era of JWST

ABSTRACT

The JWST Early Release Observations (ERO) included a NIRISS/SOSS (0.6–2.8 μm) transit of the ∼ 850 K Saturn-mass exoplanet HAT-P-18 b. Initial analysis of these data reported detections of water, escaping helium and haze. However, active K dwarfs like HAT-P-18 possess surface heterogeneities – star-spots and faculae – that can complicate the interpretation of transmission spectra, and indeed, a spot-crossing event is present in HAT-P-18 b’s NIRISS/SOSS light curves. Here, we present an extensive reanalysis and interpretation of the JWST ERO transmission spectrum of HAT-P-18 b, as well as HST/WFC3 and Spitzer/IRAC transit observations. We detect H₂O (12.5σ), CO₂ (7.3σ), a cloud deck (7.4σ), and unocculted star-spots (5.8σ), alongside hints of Na (2.7σ). We do not detect the previously reported CH₄ (log  CH₄ < −6 to 2σ). We obtain excellent agreement between three independent retrieval codes, which find a sub-solar H₂O abundance (log  H₂O ≈ −4.4 ± 0.3). However, the inferred CO₂ abundance (log  CO₂ ≈ −4.8 ± 0.4) is significantly super-solar and requires further investigation into its origin. We also introduce new stellar heterogeneity considerations by fitting for the active regions’ surface gravities – a proxy for the effects of magnetic pressure. Finally, we compare our JWST inferences to those from HST/WFC3 and Spitzer/IRAC. Our results highlight the exceptional promise of simultaneous planetary atmosphere and stellar heterogeneity constraints in the era of JWST and demonstrate that JWST transmission spectra may warrant more complex treatments of the transit light source effect.

1 INTRODUCTION

In the works for over two decades, the JWST is finally operational. Astronomers can now count on space instruments with modes designed to study exoplanetary atmospheres. Transmission spectroscopy is a commonly used method to reveal the composition and structure of an atmosphere and, therefore, to enable inferences about a planet’s formation and evolution history. During a transit, part of the starlight is blocked by the planet, whereas some light passes through the planetary atmosphere and can introduce measurable absorption features (Seager & Sasselov 2000; Brown 2001). The resulting transmission spectrum can reveal key information regarding the abundance of molecular and atomic species and the presence of clouds and hazes (e.g. Charbonneau et al. 2002; Tinetti et al. 2007; Wakeford & Sing 2015; Sing et al. 2016).

Astronomers have faced many challenges when conducting atmospheric studies with the Hubble Space Telescope (HST) and Spitzer Space Telescopes since neither observatory was designed for exoplanet observations. Numerous technical difficulties, instrument systematics, as well as narrow wavelength range and spectral resolution, have limited atmospheric inferences. For example, observations with the Wide Field Camera 3 (WFC3) instrument aboard HST are complicated by systematic trends, such as ‘HST breathing’ effects, visit-long slopes, and the ‘ramp’ effect (Wakeford et al. 2016). Nevertheless, the efforts and ingenuity of scientists have led to astonishing discoveries, such as the detection of several chemical species, including water vapour on hot Jupiters (e.g. Tinetti et al. 2007; Tsiaras et al. 2018) and even on a sub-Neptune (e.g. Benneke et al. 2019b), the inference of clouds and hazes in several gas giants (e.g. Wakeford & Sing 2015; Sing et al. 2016), and the observation of atmospheric escape on hot Neptune-like exoplanets (e.g. Ehrenreich et al. 2015; Spake et al. 2018).

HAT-P-18 b was discovered in 2010 by Hartman et al. (2011) using the Hungarian-made Automated Telescope Network. It is of approximately Saturn mass (⁠|$M = 0.197\, M\_ {\rm J}$|⁠), but with an inflated radius (⁠|$R = 0.995\, R\_ {\rm J}$|⁠), due to its higher temperature (⁠|$T_ {\rm eq}~= 852$| K) relative to the Solar system giants, which is a consequence of the planet’s comparatively short orbital period (P = 5.5 d). Ground- and space-based transmission spectroscopy has been performed on this target. Kirk et al. (2017) suggested a high-altitude haze consistent with the detection of Rayleigh scattering and the absence of the sodium absorption feature using the Auxiliary-port CAMera (ACAM) instrument on the William Herschel Telescope (WHT). Tsiaras et al. (2018) detected the presence of water vapour and a grey, opaque cloud deck in HAT-P-18 b’s atmosphere using HST/WFC3, reporting a water abundance of log H₂O = −2.63 ± 1.18 and a cloud-top pressure of log P_cloud [Pa] = 2.82 ± 0.91. Fu et al. (2022) presented an analysis of the transit observed in the Single Object Slitless Spectroscopy (SOSS) mode (Albert et al. 2023) of the Near-Infrared Imager and Slitless Spectrometer (NIRISS) instrument (Doyon et al. 2023) onboard the JWST and detected water (with an abundance of |$\log {\rm {H_2O}} = -3.03_{-0.25}^{+0.31}$|⁠), hints of methane (Δlog Z = 3.79, or a 3.2σ confidence), as well as excess helium absorption and tail in an otherwise very hazy atmosphere.

HAT-P-18 is an active K dwarf with an effective temperature of 4803 K and a slightly super-solar metallicity ([Fe/H] = 0.10 ± 0.08; Hartman et al. 2011). Stellar active regions, such as star-spots and faculae, can introduce spectral features in transmission spectra that overlap those of exoplanetary atmospheres. Occulted active regions were often masked when fitting transit models to spectroscopic light curves; however, the impact of the occulted spot on the transmission spectrum is still present despite that (e.g. Oshagh et al. 2014; Bixel et al. 2019). This can also lead to a biased transmission spectrum by impacting not only the transit depth but possibly the mid-transit time, the scaled semimajor axis, and the impact parameter (e.g. Barros et al. 2013; Alexoudi et al. 2020). Recent studies have moved towards joint inferences of transit and active region properties (e.g. Bixel et al. 2019; Espinoza et al. 2019a). The NASA Study Analysis Group on the effect of stellar contamination on space-based transmission spectroscopy (SAG 21) of the Exoplanet Exploration Program Analysis Group (ExoPAG) recommends performing these joint inferences with future observations instead of masking active region occultations (Rackham et al.2023).

Unocculted stellar active regions have long been recognized as a significant obstacle to exoplanet transmission spectroscopy. Early HST studies recognized that unocculted cool star-spots can cause strong transit depth slopes towards short-visible wavelengths (e.g. Pont et al. 2007, 2013; McCullough et al. 2014). Similarly, unocculted hot faculae can imprint a negative slope in transmission spectra (e.g. Rackham, Apai & Giampapa 2018; Kirk et al. 2021). The physical origin of this ‘stellar contamination’ is a mismatch between the intensity of the stellar surface sampled by the planet during transit and the average spectrum of the star (including the photosphere, star-spots, and faculae). Since the contrast ratio between the spectra of stellar regions at different temperatures increases as shorter wavelengths, this ‘transit light source effect’ (TLSE; Rackham, Apai & Giampapa 2018; Barclay et al. 2021) is wavelength-dependent and more significant at visible wavelengths. The most common approach to deal with the TLSE in early studies was to correct the transmission spectrum based on activity monitoring or occulted star-spot properties (e.g. Pont et al. 2008; Berta et al. 2011; Sing et al. 2011). Barstow et al. (2015) demonstrated that not accounting for star-spots when modelling transmission spectra of giant planets can bias retrieved molecular abundances, while Rackham, Apai & Giampapa (2018) further showed that the TLSE can dominate over absorption features for terrestrial planets. Moran et al. (2023) provide a recent example of this prediction, finding degenerate interpretations between unocculted star-spots and atmospheric H₂O for JWST observation of a terrestrial exoplanet.

Recent years have seen a renewed focus on incorporating unocculted stellar regions into atmospheric retrieval codes, allowing simultaneous inferences of stellar and atmospheric properties. Pinhas et al. (2018) developed a transmission spectrum retrieval framework to jointly fit a single population of unocculted stellar heterogeneities and a planetary atmosphere. Subsequent retrieval studies have explored the fidelity of TLSE retrieval assumptions (e.g. Iyer & Line 2020; Rackham & de Wit 2023; Thompson et al. 2023) and applied these joint retrievals to interpret observations from HST and Spitzer (e.g. Bruno et al. 2020; Rathcke et al. 2021), ground-based telescopes (e.g. Bixel et al. 2019; Jiang et al. 2021; Kirk et al. 2021), and JWST (Moran et al. 2023). However, the SAG 21 report (Rackham et al. 2023) notes that there is considerable scope to improve the realism and complexity of retrieval prescriptions for unocculted stellar active regions.

In this work, we aim to disentangle stellar and planetary atmosphere signals by including stellar heterogeneities in transit fits and atmospheric retrievals. We present and compare two independent atmospheric reanalyses of HAT-P-18 b, one using JWST NIRISS/SOSS Early Release Observations (ERO) transit observation and another combining transit observations from HST/WFC3 and the Infrared Array Camera (IRAC) of Spitzer. We describe the observations and data-reduction approach in Section 2 and detail our JWST NIRISS/SOSS light-curve fitting and occulted star-spot analysis in Section 3. Section 4 describes our joint stellar heterogeneity and atmospheric retrieval method and presents results from three independent retrieval codes. We summarize and discuss our results in Section 5 and conclude in Section 6.

2 OBSERVATIONS AND DATA REDUCTION

The scientific legacy of HST and Spitzer is considerable, particularly for exoplanet studies; JWST will, in many ways, build on this legacy. In this work, we present transit observations with JWST and its predecessors, HST and Spitzer, to show the potential of NIRISS/SOSS to characterize exoplanetary atmospheres and to cope with challenges such as stellar contamination.

2.1 Observations

2.1.1 JWST NIRISS/SOSS

HAT-P-18 b was observed in transit using the SOSS mode of the NIRISS instrument (Albert et al. 2023; Doyon et al. 2023) onboard JWST as part of the ERO programme (PI: Klaus M. Pontoppidan; Pontoppidan et al. 2022). The time series observation (TSO) started on 2022 June 13th at 04:36:50.861 UTC and lasted 7.15 h, which covered the 2.7 h transit as well as some baseline before and after the transit. The GR700XD grism and the CLEAR filter were used, along with the SUBSTRIP256 subarray, which captures both diffraction orders 1 and 2. There are 469 integrations, each consisting of 9 groups and lasting 54.94 s. In addition, a second exposure with the GR700XD grism and F277W filter was taken directly after the end of the main CLEAR filter TSO. This exposure had the exact same readout parameters as the CLEAR exposure, except it consisted of only 11 integrations lasting 10.07 min. The TSO was previously published by Fu et al. (2022), the major findings of which are summarized in Section 1.

2.1.2 HST/WFC3+Spitzer/IRAC

Transit observations of HAT-P-18 b were obtained using HST/WFC3 with the G141 grism in the spatial scan mode. These come from two visits, consisting of five orbits each, made as part of the Cycle 23 HST General Observer campaign (PI: Drake Deming; Deming et al. 2015) in February 2016 and January 2017. These observations spanned the wavelength range from 1.1 to 1.7 μm and were downloaded from the Mikulski Archive for Space Telescopes (MAST). The raw light curves show characteristic, non-linear systematics because of thermal effects induced by the telescope’s 96-min orbit (Foster et al. 1995). This leads us to follow standard practice (e.g. Benneke et al. 2019b) and discard from our analysis the first orbit as well as the first two measurements in each subsequent orbit. This leaves one orbit on either side of the transit and two orbits corresponding to the transit itself, including significant portions of the ingress and egress. These HST transits were previously published by Tsiaras et al. (2018).

Transit observations with Spitzer come from two visits in 2013 using the IRAC instrument in the 3.6 and 4.5 μm channels (PI: Jean-Michel Désert; Desert et al. 2012). The data were obtained from the Spitzer Heritage Archive.

2.2 Data reduction

2.2.1 JWST NIRISS/SOSS data reduction

We reduced the NIRISS/SOSS TSO using the supreme-SPOON pipeline (Coulombe et al. 2023; Feinstein et al. 2023; Lim et al. 2023; Radica et al. 2023), starting from the raw, uncalibrated files downloaded from the MAST archives. Some steps in the supreme-SPOON pipeline are handled by the official JWST pipeline. We follow a nearly identical procedure for the reduction as was presented in Feinstein et al. (2023) and Radica et al. (2023): we first process the TSO through supreme-SPOON Stage 1, which performs the detector level calibrations. This includes, in particular, the subtraction of the column-correlated 1/f noise, during which we ensure to properly mask several undispersed (zeroth order) as well as the sole dispersed (likely a second order, slightly above the target third order) contaminants of background field stars, and the ramp fitting to convert from 3D ‘ramp’ to 2D ‘slope’ images.

Stage 2 of supreme-SPOON performs further reductions on the slope-level products, including the background subtraction, for which we use the SOSS SUBSTRIP256 background model provided by the Space Telescope Science Institute (STScI),¹ and tracing of the target orders to define the extraction boxes as well as the stability of the target trace (changes in x and/or y position, changes in morphology, etc.) over the course of the TSO. For further details on these steps, please see Radica et al. (2023). We note, though, that using a constant scaling of the STScI background model, such as in Radica et al. (2023), did not perfectly remove the background step. We, therefore, separately scaled the STScI model blueward and redward of the step (e.g. Lim et al. 2023). We find values for the best-fitting background scaling to be 0.4384 and 0.4080 for pre- and post-jump, respectively.

We extract the stellar spectra using the ATOCA algorithm (Darveau-Bernier et al. 2022) using a specprofile reference file created specifically for this observation using the APPLESOSS code (Radica et al. 2022). ATOCA explicitly takes into account the overlap between the first and second orders of the target spectrum on the detector during the extraction, though the expected dilution introduced from this overlap is predicted to be near-negligible for relative measurements such as exoplanet transmission spectroscopy (Darveau-Bernier et al. 2022; Radica et al. 2022).

We use the updated APPLESOSS v2.0.0 in this work. The initial version of APPLESOSS presented in Radica et al. (2022) used the WebbPSF package (Perrin et al. 2014) to simulate the extended wings of the SOSS PSF (Point Spread Function). However, due to concerns that the simulated PSFs underpredict the SOSS wings,² we update the APPLESOSS framework to use the wings of order 1 profiles bluewards of ∼ 1.1 μm where there is no contribution from order 2 on the detector. We note, however, that the transmission spectrum is unchanged whether simulated or empirical wings are used, as the strength of the self-dilution is significantly smaller than the transit depth precision. Finally, after the extraction, we clip any data points which deviate by more than 5σ from their neighbours in the time direction; < 2 per cent of points are rejected in this way. We find no significant deviation in the target trace position with the positions included in the spectrace reference file,³ and we, therefore, use the default SOSS wavelength solution. A summary of the major reduction steps is shown in Fig. 1.

HAT-P-18 has a comoving white dwarf companion separated by only ∼ 2.66″ at a position angle of ∼ 186°, as revealed by Gaia astrometry (Mugrauer & Michel 2021).⁴ Given the aperture position angle of the SOSS observation (252.09°), the spectral trace of this companion is offset by (−9, +40) pixels relative to the trace of HAT-P-18, that is, mostly out of the extraction aperture throughout order 1 and for all wavelengths |$\lt 0.85\, \mu$|m at order 2. Combined with its measured contrast of 8 mag at |$1.4\, \mu$|m (Mugrauer & Michel 2021), it has a negligible effect on the flux extracted for HAT-P-18. No particular action was thus taken to deal with it.

2.2.2 HST/WFC3+Spitzer/IRAC data reduction and light-curve fitting

Following standard procedure for HST/WFC3 observations (e.g. Deming et al. 2013; Benneke et al. 2019b), we perform the necessary data reduction and light-curve fitting by using the modular Exoplanet Transit Eclipses and Phase curves (ExoTEP) framework (Benneke et al. 2017, 2019a), which employs the batman transit light-curve model (Kreidberg 2015). Using a Markov chain Monte Carlo (MCMC) method with the emcee package (Foreman-Mackey et al. 2013), ExoTEP jointly fits the transit and systematics models for the two visits along with photometric noise parameters. We construct spectrophotometric light curves by spectrally binning at 30 nm intervals (Kreidberg et al. 2014) the extracted flux for each visit. This binning was chosen because it offered the best trade-off between the number of data versus error bar size, as compared to 10 and 60 nm bin widths. We follow Benneke et al. (2019b) for the procedure of the white and spectrophotometric light-curve fits. The latter for both visits are shown in Appendix Fig. B1, whereas the best-fitting values from the white light-curve fit are quoted in Appendix Table A1.

We follow standard procedure for Spitzer/IRAC image processing (e.g. Benneke et al. 2019b). As was done for the HST data, we employ ExoTEP to reduce the data further and fit the Spitzer light curves. After correcting for the Spitzer-specific systematics, we obtain a final light curve for each of the two channels as displayed in Appendix Fig. B2. We used 80-s bins for the time axis due to the relatively high-cadence of the Spitzer data. Compared to the HST spectrophotometric light curves, we find that the scatter is higher, resulting in larger errors for the Spitzer data points in the combined transmission spectrum. We kept the transit depths from Spitzer/IRAC in the HST transmission spectrum in order to better constrain the atmospheric retrievals.

3 JWST NIRISS/SOSS LIGHT-CURVE FITTING AND OCCULTED STAR-SPOT ANALYSIS

3.1 Light-curve fitting with spot-crossing masked

For the SOSS data, we construct separate white light curves for orders 1 and 2 by summing the flux from all wavelengths in order 1 (0.85–2.8 μm) and from 0.6 to 0.85 μm in order 2. We then fit a transit model to each white light curve using the juliet package (Espinoza, Kossakowski & Brahm 2019b). We fix the orbital period to 5.508023 d, the eccentricity to 0.084, and the argument of periastron to 120° (Hartman et al. 2011; all their values used in this paper are listed in Table 1), and put wide, uninformative priors on the time from mid-transit, t₀, the impact parameter, b, the scaled orbital semimajor axis, a/R_*, and the scaled planet radius, R_p/R_*. We fit two parameters of the quadratic limb darkening law following the parametrization of Kipping (2013), a scalar jitter term, σ, which replaces the flux errors reported by the reduction pipeline, and finally, a term to fix the zero point of the transit baseline. Unlike previous SOSS TSOs (Feinstein et al. 2023; Radica et al. 2023), we find that no detrending is necessary, as the white light curves for both orders are best fit by a transit model with no additional systematics models: Δlog Z = 26.1, or a >5σ preference by the Benneke & Seager (2013) scale.

There is a spot-crossing event clearly visible just after the time from mid-transit. Here we do not include a model of this occulted star-spot and instead mask the integrations associated with the spot-crossing (integrations 240–270). Therefore, we fit eight parameters for each order. The reduced Chi-squared statistic for the fits is |$\chi ^2_\nu$| = 1.08 for order 1 and 1.06 for order 2. The best-fitting transit models for orders 1 and 2 are overplotted in Fig. 2.

We then proceed to fit the spectrophotometric light curves at the pixel level – that is, we fit one light curve per pixel column on the detector. This results in 2038 light curves for order 1 and 567 for order 2. For these fits, we fix the orbital parameters and the transit baseline’s zero point to their best-fitting values from the order 1 white light-curve fit because of the better signal-to-noise (S/N). We leave only the scaled planet radius, limb-darkening parameters, and jitter term free. We put Gaussian priors on the limb-darkening parameters based on calculations from the ExoTiC-LD package (Wakeford & Grant 2022) using the 3D stagger grid (Magic et al. 2015). The widths of the Gaussian priors are set to 0.1 (e.g. Espinoza et al. 2022; Pontoppidan et al. 2022). We tested larger widths of 0.2 following Patel & Espinoza (2022) and the retrieved transmission spectrum is consistent within less than 1σ.⁵ We also experimented with freely fitting the limb-darkening coefficients and found that this results in a slight (∼20 ppm) offset to the spectrum while preserving the relative amplitudes of the spectral features. We again mask the integrations associated with the spot-crossing in each fit. The resulting transmission spectrum is displayed in Fig. 3 with the HST/WFC3 one from Section 2.2.2.

Although we do not directly use the F277W exposure for science purposes in this study, as in Radica et al. (2023), it is useful to pinpoint undispersed (order 0) contaminants of field stars. Several bright contaminants are clearly visible in the F277W data frame, many of which are also readily visible in the CLEAR frame. However, unlike in Radica et al. (2023), we are unable to post-process our transmission spectrum to correct for the diluting effects of field star contamination. The contaminants are too bright preventing good reconstruction of the target trace. Moreover, as pointed out by Fu et al. (2022), one contaminant is time varying. We, therefore, mask the affected regions of the transmission spectrum as was done in Feinstein et al. (2023) and Fu et al. (2022). The wavelength regions masked in this way are as follows: 0.714–0.724 and 0.841–0.85 μm in order 2, and 0.853–0.870, 1.048–1.061, 1.366–1.384, and 1.972–2.011 μm in order 1. A frame of the F277W exposure, highlighting the masked regions is shown in Appendix Fig. B3.

We also produce a transmission spectrum using an independent reduction pipeline, NAMELESS (Coulombe et al. 2023; Feinstein et al. 2023; Radica et al. 2023), as shown in Appendix Fig. A1, to verify the self-consistency of our results (see Appendix  A). The best-fitting transit parameters from NIRISS/SOSS white light-curve fits are shown in Appendix Table A1.

3.2 Occulted star-spot method

We also perform a second light-curve fit, enabling a joint inference of the star-spot and planet properties to model the spot and study its impact on the transmission spectrum. To this end, we first compute a single broad-band light curve by summing the flux from wavelengths 0.65–0.85 μm of order 2 together with wavelengths 0.85–1.5 μm of order 1. We exclude wavelengths bluewards of 0.65 μm to maximize the S/N and wavelengths redwards of 1.5 μm because the effect of spot crossings is stronger at shorter wavelengths where the spot contrast with respect to the stellar photosphere is larger. We fit a transit model with a spot-crossing event using spotrod (Béky, Kipping & Holman 2014), which we have implemented into the juliet package (Espinoza, Kossakowski & Brahm 2019b). We fix the period, the eccentricity and the argument of periastron to the Hartman et al. (2011) values and fit the mid-transit time, the impact parameter, the scaled semimajor axis, the scaled planet radius, the two parameters of the quadratic limb-darkening law, and the jitter term with wide, uninformative priors. We also fit the zero point of the transit baseline, though the best-fitting value is again very close to zero. We fit four additional parameters to model the star-spot: the x- and y-positions of the spot, its radius, and its spot-to-stellar flux contrast. The positions and the radius of the spot are in stellar radii units. The centre of the star is at (0, 0). We use uniform priors ranging from 0 to 1. We employ dynesty (Speagle 2020) to sample the parameter space with 5000 live points. The reduced Chi-squared statistic for the fit with the highest likelihood is χ_ν = 1.08. There is strong evidence (Δlog Z = 95.7, or a > 5σ preference) for a transit model with an occulted spot. The best-fitting parameters for the broad-band light-curve fit with the spot-crossing modelled are displayed in Appendix Table A1. The choice of wavelength range (0.65–1.5 μm) results in best-fitting values that are more precise and overall consistent at 1σ to the best-fitting values with the entire wavelength ranges of order 1 and 2. We note that the y-position of the spot is not well constrained, as shown in the corner plot in Appendix Fig. B4, so we use the set of parameters with the highest likelihood instead of the medians of the posteriors from the broad-band light-curve fit, to ensure that the parameters are self-consistent. This best-fitting transit model is overplotted in the top panel of Fig. 4, and a cartoon of the spot-crossing event is shown in the bottom panel.

We then fit the spectrophotometric light curves at a resolving power of R = 100, fixing the orbital parameters (T₀, b, a/R_*), as well as the spot positions and radius from the above broad-band light-curve fit. We fit the scaled planet radius, the limb-darkening parameters, the spot contrast and the jitter coefficient with the priors described above. The spectrophotometric light curves for 14 bins with their best-fitting transit models are shown in Fig. 5 for the highest likelihood model. The spot contrast is correlated with the y-position (see Appendix Fig. B4) because there is a degeneracy between the radius and the contrast of the spot given the strong correlation between the spot temperature and filling factor⁶ (Bruno et al. 2022). To capture the uncertainties coming from the degeneracies in the spot y-position, size and contrast, we repeat the fit of the spectrophotometric light curves 50 times by fixing the same parameters (T₀, b, a/R_*, spot positions and radius) to 50 random sets taken from the posterior samples of the broad-band light-curve fit. We draw the 50 sets from the 68 per cent of the weighted samples with the highest likelihoods.

We next constrain the temperature of the occulted spot by fitting PHOENIX stellar atmosphere model spectra (Husser et al. 2013) to the spot contrast spectrum resulting from the spectrophotometric light-curves fits. Star-spots’ spectra are thought to be represented by stellar models with lower temperatures and 0.5–1 dex lower log g than the host star. The increased magnetic pressure in active features decreases the gas pressure (Solanki 2003; Bruno et al. 2022), which is akin to a stellar model with a lower surface gravity. In fitting for the spot temperature, we thus also left the spot surface gravity to vary as a free parameter. For the spot, we use stellar model spectra with temperatures from 4000 to 5000 K, logarithmic surface gravities from 1 to 5 dex and a fixed metallicity of 0.1 from Hartman et al. (2011). We interpolated these spectra in temperature and log g as needed during the fitting process. For the star, we use a stellar model spectrum with the temperature, log g and metallicity fixed to the corresponding values from Hartman et al. (2011). The model spot contrast spectrum is then simply the ratio of the model spot spectrum to the model star spectrum. We fit the temperature of the spot and the log g with wide, uninformative priors using the emcee MCMC package (Foreman-Mackey et al.2013).

3.3 Inferred occulted star-spot properties on HAT-P-18

From the aforementioned process, we obtained 51 star-spot contrast spectra from the 51 spectrophotometric fits performed: one from the highest likelihood values and 50 with random draws from the broad-band light-curve fit. Each model has its own set of parameters (i.e. T₀, b, a/R_*, spot positions and radius). The resulting density of spot contrast spectra is shown in Fig. 6. By inspecting the individual contrast spectra, we see that 28 of the 50 random fitting models (56 per cent) have a retrieved contrast spectrum within 1σ of the one retrieved with the highest likelihood set of parameters. We lump these solutions together as the ‘most likely’ one (blue model in Fig. 6). The mean and standard deviation of the spot x- and y-position, radius, temperature, and the model log g are given in Table 2 for this group of solutions. This corresponds to a spot radius of 0.116 ± 0.014 stellar radii, a temperature ΔT = −93 ± 15 K colder than the star (or T_spot = 4710 ± 15 K) and a lower log g by 1.16 ± 0.19 dex. The retrieved log g for the star-spot model is consistent with Bruno et al. (2022) expectation.

Among the remaining contrast spectra, those not within 1σ of the one with the highest likelihood, we can identify three other broad families of solutions corresponding to a colder spot because of a smaller filling factor; see Table 2. The transmission spectra for these four families of solutions do not significantly change (they are consistent at the 1σ level) from the transmission spectrum obtained with the spot occultation features masked (as done in subsection 3.1). This is expected for a spot of that size and temperature, as the TLSE from occulting such a spot results only in a small slope towards bluer wavelengths in the transmission spectrum, as shown in Appendix Fig. B5. This effect, estimated to be of the order of 10–20 ppm from equation (4) in Section 4, is smaller than our transit depth uncertainties.

4 RETRIEVAL ANALYSIS

We now present the inferences from our retrieval analysis of HAT-P-18 b’s transmission spectrum with NIRISS/SOSS (from subsection 3.1). Our retrievals jointly model the influence of the planetary atmosphere and unocculted stellar heterogeneities on the transmission spectrum and thus yield simultaneous constraints. We conducted this analysis with three independent retrieval codes to ensure robust results.

In what follows, we first describe our approach for joint retrievals of a planetary atmosphere and stellar contamination. We then describe the setup of our three retrieval codes before presenting our resulting constraints on HAT-P-18 b’s atmosphere and unocculted stellar active regions. Finally, we compare retrieval results from our HST+Spitzer transmission spectrum to those yielded from NIRISS/SOSS.

4.1 A joint planetary atmosphere and stellar contamination retrieval method

The influence of active stellar regions on exoplanet transmission spectra – the TLSE (e.g. Rackham, Apai & Giampapa 2018) – can be modelled simultaneously with the planetary atmosphere in retrievals. This has the advantage of yielding joint constraints on the planetary atmosphere and unocculted stellar heterogeneities while also propagating uncertainties from the star into the derived atmospheric constraints (see Rackham et al. 2023, for a review). Several retrieval studies have included prescriptions for unocculted star-spots and/or faculae (e.g. Pinhas et al. 2018; Bixel et al. 2019; Iyer & Line 2020; Jiang et al. 2021; Rathcke et al. 2021; Rackham & de Wit 2023; Thompson et al. 2023), with the most common being a three-parameter model that fits for the temperature and covering fraction of a single heterogeneity (Pinhas et al. 2018). While these treatments have generally sufficed for HST and ground-based data, the exceptional data quality and wavelength coverage of JWST motivates the consideration of more complex stellar contamination models.

We investigate a range of unocculted stellar heterogeneity prescriptions during our HAT-P-18 b retrievals. Our goal is to determine the appropriate level of star-spot and/or faculae model complexity required to interpret the JWST/NIRISS transmission spectrum of HAT-P-18 b while also informing future JWST retrieval studies of hot Jupiters. We begin by summarizing the equations underlying the TLSE

4.1.1 The TLSE

The transmission spectrum of an exoplanet transiting a star with a heterogeneous stellar surface can be written as (MacDonald & Lewis 2022)

where Δ_λ is the observed transmission spectrum, |$\delta _{\lambda , \, \rm {atm}}$| is the transmission spectrum from the planetary atmosphere, |$\epsilon _{\lambda , \, \rm {het}}$| is the wavelength-dependent ‘stellar contamination factor’ from a heterogeneous stellar surface, and |$\psi _{\lambda , \, \rm {night}}$| accounts for thermal emission from the planetary nightside (we assume this final term negligibly deviates from unity for HAT-P-18 b). For a planetary atmosphere with properties varying only with altitude (the 1D assumption), |$\delta _{\lambda , \, \rm {atm}}$| is given by

where |$R_{\mathrm{p, \, top}}$| is the planetary radius at the top of the modelled atmosphere, b is the ray impact parameter, R_* is the stellar radius, and |$\tau _{\lambda , \, \rm {slant}}$| is the slant optical depth. For a stellar surface with N_het unocculted heterogeneous active regions (e.g. spots or faculae), the stellar contamination factor, |$\epsilon _{\lambda , \, \rm {het}}$|⁠, can be expressed as

where |$f_{\mathrm{het}, \, i}$| is the fractional stellar disc coverage of the i^th heterogeneous region, with corresponding specific intensity |$I_{\lambda , \, \mathrm{het}, \, i}$|⁠, while |$I_{\lambda , \, \rm {phot}}$| is the specific intensity of the stellar photosphere. Equation 3 demonstrates that the stellar contamination factor deviates from unity when a heterogeneity possesses a different intensity from the photosphere. In general, |$I_{\lambda , \, \mathrm{het}, \, i}$| is shorthand for |$I_{\lambda , \, \mathrm{het}, \, i} \, (\mathrm{[Fe/H]}_i, \mathrm{log} \, g_i, T_i)$|⁠, where [Fe/H]_i is the local metallicity of the heterogeneity, |$\mathrm{log} \, g_i$| is the local log surface gravity, and T_i is the heterogeneity temperature.

Most retrieval studies consider a single population of unocculted heterogeneities with a temperature different from that of the photosphere (but with the same metallicity and surface gravity). Given these assumptions, the stellar contamination factor is given by (e.g. Pinhas et al. 2018; Rackham, Apai & Giampapa 2018)

where the three free parameters are the heterogeneity coverage fraction, f_het, the heterogeneity temperature, T_het, and the photosphere temperature, T_phot.

We also consider here a more general treatment of stellar contamination. Motivated by several strategic gaps identified in the NASA ExoPAG SAG 21 community report (Rackham et al. 2023), we also consider a two-heterogeneity model with the assumption of a common surface gravity relaxed. By defining these heterogeneities as star-spots and faculae (T_spot < T_phot < T_fac), we can write the stellar contamination factor as

The addition of a second heterogeneity population and separate surface gravities increases the number of free parameters for this stellar contamination model to eight.

4.1.2 Unocculted stellar heterogeneity models for HAT-P-18 b

We consider five treatments for unocculted stellar heterogeneities when retrieving HAT-P-18 b’s transmission spectrum:

• No stellar contamination.

• One heterogeneity, defined by three free parameters: f_het, T_het, and T_phot.

• Two heterogeneities (spots+faculae), defined by five free parameters: f_spot, T_spot, f_fac, T_fac, and T_phot.

• One heterogeneity with free surface gravity, defined by five free parameters: f_het, T_het, |$\mathrm{log} \, g_{\mathrm{het}}$|⁠, T_phot, and |$\mathrm{log} \, g_{\mathrm{phot}}$|⁠.

• Two heterogeneities with free surface gravity, defined by eight free parameters: f_spot, T_spot, |$\mathrm{log} \, g_{\mathrm{spot}}$|⁠, f_fac, T_fac, |$\mathrm{log} \, g_{\mathrm{fac}}$|⁠, T_phot, and |$\mathrm{log} \, g_{\mathrm{phot}}$|⁠.

4.2 Retrieval configuration

We applied three exoplanet retrieval codes to HAT-P-18 b’s transmission spectrum: Poseidon (MacDonald & Madhusudhan 2017; MacDonald 2023), Aurora (Welbanks & Madhusudhan 2021), and SCARLET (Benneke & Seager 2012; Benneke 2015). We initially conducted independent analyses to establish which atmospheric and stellar properties provide the necessary complexity to describe HAT-P-18 b’s transmission spectrum. After comparing our results, we chose a common set of model assumptions for the three codes: chemical opacity from Na, K, H₂O, CO, CO₂, CH₄, HCN, and NH₃ (the most prominent opacity sources expected at HAT-P-18 b’s equilibrium temperature under equilibrium chemistry and/or with chemical quenching, see for example Madhusudhan et al. 2016); an inhomogeneous grey cloud-deck combined with a power-law haze; and one unocculted stellar heterogeneity (i.e. the one heterogeneity model). Our retrievals were conducted on the R = 100 binned variant of HAT-P-18 b’s transmission spectrum with the spot-crossing masked (see Fig. 3), but we found consistent results when retrieving the full-resolution NIRISS data. We describe the specific configuration used by each retrieval code below and list the priors used in Table 3.

4.2.1 POSEIDON

We conducted a series of Poseidon retrievals to assess the model complexity required to fit HAT-P-18 b’s JWST/NIRISS transmission spectrum. This includes the five stellar contamination models listed in subsection 4.1.2, nested models to compute detection significances for each chemical species, and additional robustness tests (e.g. a retrieval of the unbinned pixel-resolution spectrum).

Our Poseidon retrievals model HAT-P-18 b has a H₂+He-dominated atmosphere with a standard configuration used for hot Jupiter retrieval studies (e.g. MacDonald & Madhusudhan 2017; Kirk et al. 2021; Taylor et al. 2023). We prescribe 100 layers, spaced uniformly in log-pressure from 10⁻⁸–100 bar, where the layers follow a variant of the pressure–temperature (P-T) profile parametrization from Madhusudhan & Seager (2009) but with the reference temperature located at 10 mbar. The inhomogeneous cloud and haze aerosol model follows MacDonald & Madhusudhan (2017). We fit for the log₁₀ volume mixing ratios of the 8 gases in the reference model, using absorption cross-sections (see MacDonald & Lewis 2022) computed from the following line lists: Na and K (Ryabchikova et al. 2015); H₂O (Polyansky et al. 2018); CO (Li et al. 2015); CO₂ (Tashkun & Perevalov 2011); CH₄ (Yurchenko et al. 2017); HCN (Barber et al. 2014); and NH₃ (Coles, Yurchenko & Tennyson 2019). We include H₂–H₂ collision-induced absorption from HITRAN (Karman et al. 2019). We also fit for the planetary radius at a reference pressure of 10 bar. We calculate model transmission spectra from 0.55–2.9 μm at R = 20 000. We calculate the stellar contamination factor (equation (3)) by interpolating PHOENIX models (Husser et al. 2013) using the PyMSG package (Townsend & Lopez 2023). Our Poseidon retrievals have up to 27 free parameters (depending on the chosen stellar contamination model from subsection 4.1.2), with 1000 nested sampling live points used by the pymultinest (Feroz, Hobson & Bridges 2009; Buchner et al. 2014) package to chart the parameter space.

4.2.2 SCARLET

We performed retrievals using the SCARLET retrieval framework (Benneke 2015; Benneke et al. 2019a). We improved the SCARLET framework from previous work on low resolution observations (Benneke et al. 2019a, b; Piaulet et al. 2023) with the addition of non-uniform cloud coverage and of the contribution of stellar contamination to the transmission spectrum.

In the SCARLET retrieval, the atmosphere is parametrized by the abundances of the spectrally-active chemical species of interest (Na, K, H₂O, CO, CO₂, CH₄, HCN, NH₃) which are assumed to be well mixed and constant with pressure, and an isothermal temperature structure. The ratio of He and H₂, which make up the remainder of the gas, is set to Jupiter’s He/H₂ of 0.157.

The retrieval includes three parameters describing aerosols (log p_cloud, f_cloud, log c_haze). The cloud parametrization consists of a grey cloud, opaque across all wavelengths, with a cloud top pressure p_cloud. Potential fractional cloud coverage is captured by the f_cloud parameter and corresponds to a weighted average of a cloud-free and a cloudy model. The contribution of haze to the transmission spectrum via a slope towards shorter wavelengths is parametrized using the c_haze parameter, which is a factor that multiplies the Rayleigh contribution to the scattering coefficient (unity corresponds to standard Rayleigh scattering).

Finally, two parameters describe the stellar heterogeneity: ΔT_spot and f_spot. Stellar heterogeneity is implemented following the standard approach (see subsection 4.1.1) assuming that the contribution of stellar contamination to the spectrum can be represented by that of spots with a temperature of T_spot = T_phot+ΔT_spot (where T_phot is the photosphere temperature) covering a fraction f_spot of the star. Stellar models are queried from the PHOENIX grid of stellar atmosphere models (Husser et al. 2013).

The forward models calculated by the retrieval have a resolving power of 15 625. We use opacity sampling from cross-section tables at a resolving power of 250 000 to build the opacity tables used, and the opacities are taken from the ExoClimes Simulation Platform data base. We choose HITEMP opacities for H₂O, CO₂, CH₄, NH₃ (Rothman et al. 2010), and ExoMol opacities for all other species (Tennyson et al. 2016). We use Nested Sampling to sample the parameter space and its python implementation in the nestle module (Skilling 2004; Mukherjee, Parkinson & Liddle 2006; Shaw, Bridges & Hobson 2007).

4.2.3 AURORA

We also perform atmospheric retrievals with Aurora (Welbanks & Madhusudhan 2021) to interpret the transmission spectrum of HAT-P-18 b. The application of Aurora to transmission spectra of exoplanets is described in Welbanks & Madhusudhan (2022) following the general setup presented in Welbanks & Madhusudhan (2019) and Welbanks et al. (2019). We follow the methods of Pinhas et al. (2018) to consider the impact of stellar heterogeneities as implemented in previous works (see e.g. Ahrer et al. 2023).

Aurora computes line-by-line radiative transfer in transmission geometry for a parallel-plane atmosphere assuming hydrostatic equilibrium. Our atmospheric model is computed using a 100 layer pressure grid uniformly spaced in log-pressure between 10⁻⁸ and 100 bar, and a wavelength grid from 0.55 μm to 2.9 μm, covering this NIRISS observation, at a constant resolution of R = 20 000. We parametrize the P–T profile of the atmosphere using the prescription from Madhusudhan & Seager (2009) and the presence of inhomogeneous clouds (e.g. Line & Parmentier 2016) and hazes (e.g. Lecavelier Des Etangs et al. 2008), using a single sector for their combined spectroscopic imprint following the description in Welbanks & Madhusudhan (2021). The volume mixing ratios for the eight chemical species considered in our fiducial model are assumed to be constant with height and use cross-sections for H₂–H₂ and H₂–He collision induced absorption (CIA; Richard et al. 2012), CH₄ (Yurchenko & Tennyson 2014), CO (Rothman et al. 2010), CO₂ (Rothman et al. 2010), H₂O (Rothman et al. 2010), HCN (Barber et al. 2014), K (Allard, Spiegelman & Kielkopf 2016), Na (Allard et al. 2019), and NH₃ (Yurchenko, Barber & Tennyson 2011), computed as described in Gandhi & Madhusudhan (2017, 2018) and Welbanks et al. (2019). Following Pinhas et al. (2018), we model the heterogeneous stellar photosphere by interpolating stellar models from the PHOENIX data base (Husser et al. 2013). The parameter estimation is performed using the nested sampling algorithm (Skilling 2004) with multinest (Feroz, Hobson & Bridges 2009) and its python implementation pymultinest (Buchner et al. 2014).

4.3 An explanation for HAT-P-18 b’s transmission spectrum

We first provide an intuitive explanation for HAT-P-18 b’s transmission spectrum. Fig. 7 shows a spectral decomposition of the maximum likelihood (best-fitting) model transmission spectrum from the Poseidon one heterogeneity retrieval (model (ii) in subsection 4.1.2). We detect multiple H₂O absorption features near 0.95 μm, 1.15 μm, 1.4 μm, 1.9 μm, and 2.6 μm (with a combined significance of 12.5σ). We further detect a CO₂ absorption feature near 2.7 μm (7.3σ) and infer evidence of Na at lower significance (2.7σ). The small amplitude of the absorption features is due to an optically thick cloud deck which is uniform around the terminator (7.4σ). The inference of non-patchy clouds is driven by the wing shape of molecular bands (MacDonald & Madhusudhan 2017) – especially the H₂O band centred near 1.4 μm – and the flat continuum from 1.6–1.8 μm. Finally, the spectral slope shortwards of 1.65 μm is caused by the combination of unocculted star-spots (5.8σ) and the uniform cloud deck. The quoted detection significances are from Bayesian model comparisons between the reference Poseidon one heterogeneity model and another retrieval with one model component removed (e.g. Benneke & Seager 2013). We do not identify statistically significant evidence for K, CO, CH₄, HCN, NH₃, nor a scattering haze. We note that these conclusions are independently retrieved by the SCARLET and Aurora retrievals, which produce similar best-fitting spectra (see Section 4.4). Finally, we determine below that a single high data point near 1.1 μm, not fit by our retrieval models, is attributable to metastable helium. We verified that removing this data point does not alter the retrieval results.

4.3.1 Evidence of helium absorption

We assessed evidence of helium absorption in HAT-P-18 b’s atmosphere via an independent fit to the pixel-resolution transmission spectrum near 1.083 μm. Fig. 8 shows that a Gaussian fit, with free amplitude and width, favours an excess absorption of 0.13 ± 0.03 per cent (4.3σ) centred around the near-infrared helium triplet at 1.083 μm. Our result is consistent within 1σ to Fu et al. (2022), Vissapragada et al. (2022), and Paragas et al. (2021). Since the helium line profile is unresolved with NIRISS/SOSS, even at pixel resolution, we do not conduct more complex modelling of the helium line. Consequently, we cannot unambiguously attribute this signature to the planetary atmosphere nor disentangle the layers probed by the helium triplet. Nonetheless, we confirm that HAT-P-18 b exhibits a longer transit duration with more absorption post transit – as also observed by Fu et al. (2022) – which could indicate an extended atmosphere with the presence of a cometary-like tail as exhibited by WASP-107 b (Allart et al. 2019).

Future observations of HAT-P-18 b at higher spectral resolution can confirm the helium absorption. A confirmation would indicate that HAT-P-18 b is losing its atmosphere and is undergoing significant mass-loss. We estimated the intrinsic excess absorption arising from the resolved helium triplet to be 3.17 ± 0.67 per cent (Figure  8). This value was derived by fitting a Gaussian with a fixed width of 1 Å (the typical width of the helium triplet; e.g. Allart et al. 2018; Nortmann et al. 2018; Salz et al. 2018; Allart et al. 2019) and a free intrinsic amplitude, which was then convolved to the pixel-resolution NIRISS/SOSS data. Such a strong helium signature would be readily detectable with ground-based high-resolution spectrographs, which would provide complementary observations to study the unique processes shaping HAT-P-18 b’s upper atmosphere.

4.4 The atmosphere of HAT-P-18 b

Our retrieved atmospheric properties for HAT-P-18 b are summarized in Fig. 9. All three retrieval codes obtain retrieved transmission spectra consistent with the picture presented previously: HAT-P-18 b’s observed spectrum is shaped by absorption from H₂O and CO₂ (and likely Na), unocculted star-spots, and a cloud deck. The posterior distributions in Fig. 9, and corresponding confidence regions in Table 4, provide our quantitative constraints on HAT-P-18 b’s atmospheric properties.

We obtain precise constraints on several atmospheric properties, including the H₂O and CO₂ abundances, cloud deck pressure, and atmospheric temperature. The retrieved H₂O abundance (⁠|$\log X_{\mathrm{H_2 O}} = -4.47_{-0.28}^{+0.30}$| (Poseidon); |$-4.11_{-0.22}^{+0.27}$| (SCARLET); |$-4.65_{-0.22}^{+0.24}$| (Aurora)) is ∼ 10× lower than the expected abundance for a solar-composition atmosphere under chemical equilibrium at HAT-P-18 b’s equilibrium temperature (Hartman et al. 2011). Conversely, the retrieved CO₂ abundance (⁠|$\log X_{\mathrm{CO_2}} = -4.86_{-0.45}^{+0.45}$| (Poseidon); |$-4.39_{-0.26}^{+0.38}$| (SCARLET); |$-5.16_{-0.36}^{+0.40}$| (Aurora)) is > 100× higher than expected for a similar solar-composition atmosphere (cf. Woitke et al. 2018). The retrieved Na abundance (⁠|$\log X_{\mathrm{Na}} = -6.99_{-2.09}^{+0.59}$| (Poseidon); |$-5.45_{-1.59}^{+0.69}$| (SCARLET); |$-6.65_{-1.68}^{+0.42}$| (Aurora)) is more weakly constrained, being broadly consistent with a range of sub-solar to solar (Poseidon and Aurora) or even a super-solar composition (SCARLET). An opaque cloud deck, consistent with uniform terminator coverage (f_cloud > 0.83, see Table 4), is localized with a cloud top pressure of ≈ 30 mbar (⁠|$\log P_{\mathrm{cloud}} = -1.52_{-0.13}^{+0.10}$| (Poseidon); |$-1.53_{-0.13}^{+0.10}$| (SCARLET); |$-1.40_{-0.10}^{+0.07}$| (Aurora)). Finally, the retrieved P-T profile, over the pressure range probed by our NIRISS/SOSS observation, is isothermal at ≈ 600 K (⁠|$T_{\mathrm{ref}} = 616_{-78}^{+99}$| K (Poseidon); |$523_{-65}^{+81}$| K (SCARLET); |$673_{-85}^{+90}$| K (Aurora)). We discuss the sensitivity of these inferences to the assumed stellar contamination model in subsection 4.5.3. We note that, while SCARLET favours slightly higher abundances and lower temperatures than Poseidon and Aurora, all three codes produce consistent retrieval results to within 1σ. The small differences between our retrieval results may be due to the slightly different model configuration and priors used by SCARLET (see Table 3).

We also obtain robust upper limits on the abundances of several important chemical species. Our retrievals strongly disfavour the presence of K (log X_K ≲ −10 to 2σ), which could be due to condensation out of the gas phase. Similarly, we establish that CH₄ is at least 100 × less abundant than the expectation for a solar-composition atmosphere in chemical equilibrium (⁠|$\log X_{\mathrm{CH_4}} \lesssim -6$| to 2σ versus the expected |$\log X_{\mathrm{CH_4}} \sim -4$|⁠; Woitke et al. 2018), which could indicate photochemical dissociation (e.g. Moses et al. 2011; Venot et al. 2012). Although we see a tentative hint of CO in our posteriors (see Fig. 9), additional observations at longer wavelengths, covering the 4.6 μm CO band, would be required to robustly constrain the CO abundance (and hence the atmospheric C/O ratio). We also find upper limits on HCN and NH₃ (log X_HCN ≲ −5; |$\log X_{\mathrm{NH_3}} \lesssim -6$|⁠), which are consistent with both equilibrium and disequilibrium expectations for hot giant planets (e.g. Moses et al. 2011; Venot et al. 2012). These upper limits underline that SOSS data is of sufficiently high quality to provide scientifically meaningful constraints even on the abundances of non-detected chemical species.

4.5 Unocculted stellar heterogeneities

We next turn to constraints on unocculted stellar heterogeneities. We first present results under the assumption of a single unocculted heterogeneity population with fixed log g – a common approach in the literature – before relaxing this assumption by considering more complex treatments of the TLSE. Finally, we examine the impact of different unocculted stellar heterogeneity models on the retrieved atmospheric properties of HAT-P-18 b.

4.5.1 One heterogeneity stellar contamination model

Assuming a single population of unocculted stellar active regions, our retrievals require unocculted star-spots to explain HAT-P-18 b’s transmission spectrum (see Fig. 9 and Table 4). The spots are ≈ 300 K cooler than the photosphere (⁠|$T_{\mathrm{het}} = 4487_{-84}^{+79}$| K (Poseidon); |$4408_{-12}^{+16}$| K (SCARLET); |$4513_{-92}^{+85}$| K (Aurora), compared to T_phot ≈ 4800 K) and cover ≈ 10 per cent of the stellar surface (⁠|$f_{\mathrm{het}} = 0.10_{-0.02}^{+0.03}$| (Poseidon); |$0.10_{-0.01}^{+0.01}$| (SCARLET); |$0.12_{-0.03}^{+0.04}$| (Aurora)). We note that SCARLET has smaller uncertainties on the retrieved stellar parameters, which we attribute to the different parametrization of the stellar heterogeneity and hazes in the prior.

4.5.2 More complex stellar contamination models

We find that the best-fitting stellar contamination model for HAT-P-18 b is two distinct populations of unocculted heterogeneities – spots and faculae – with different surface gravities. We established this via a series of Bayesian model comparisons between the five stellar contamination models listed in subsection 4.1.2. Our results, summarized in Table 5, demonstrate a significant improvement in the Bayesian evidence for the two heterogeneities with free surface gravity model compared to the fiducial one heterogeneity model we have focused on thus far (a Bayes factor of 77.5, equivalent to 3.4σ). We also find that the preference for unocculted heterogeneities versus no stellar contamination is greater for the two heterogeneities with free surface gravity model compared to the one heterogeneity model (6.5σ versus 5.8σ). Nevertheless, we show below that adopting the preferred stellar contamination model does not significantly alter the atmospheric constraints for HAT-P-18 b presented in Fig. 9.

We show in Fig. 10 that the choice of stellar contamination model can significantly alter inferences about unocculted stellar active regions. The most significant change is that retrievals with free log g favour much cooler spots (ΔT ≈ −800 K versus ≈ −300 K) than those with log g_spot = log g_phot. Indeed, the retrieved spot surface gravity is lower than the photosphere by > 2σ. We similarly find that fitting for log g_fac allows a more accurate determination of the faculae temperature compared to fixing the surface gravities. Overall, our preferred solution for HAT-P-18’s unocculted active regions is as follows: spots and faculae cover |$f_{\mathrm{spot}} = 12^{+5}_{-3}$| per cent and f_fac > 4 per cent (2σ lower limit) of HAT-P-18’s surface, with temperatures of |$T_{\mathrm{spot}} = 4045^{+162}_{-116}$| K and |$T_{\mathrm{fac}} = 4952^{+62}_{-64}$| K and surface gravities of log g_spot < 4.17 (2σ upper limit) and log g_fac > 4.45 (2σ lower limit). These compare to a background photosphere with |$T_{\mathrm{phot}} = 4840^{+50}_{-51}$| K and |$\log g_{\mathrm{phot}} = 4.56^{+0.03}_{-0.03}$|⁠. These results demonstrate that more complete stellar contamination models can provide deeper insights into the nature of unocculted stellar active regions.

4.5.3 Sensitivity of atmospheric properties to stellar contamination models

Finally, we find that, provided the spectroscopic imprint of unocculted heterogeneities is considered in the atmospheric model, the adopted stellar contamination model only minimally impacts the retrieved atmospheric properties for this data. Fig. 11 shows that our preferred stellar contamination model (two heterogeneities with free surface gravity) produces retrieved temperatures, cloud pressures, and H₂O and CO₂ abundances consistent with the simpler spot-only model with fixed surface gravity (one heterogeneity). However, a retrieval not including stellar contamination would attribute the spectral slope (caused by spots and a cloud deck) to an atmospheric haze and erroneously retrieve a H₂O abundance 1–2σ higher than when spots are included. Given that we do not see evidence of a haze when including a stellar heterogeneity model, we conclude that our NIRISS/SOSS data have sufficient wavelength coverage and precision to lift the degeneracy between an atmospheric haze power law and unocculted star-spots. This mirrors a conclusion from Rackham et al. (2023), where they demonstrated that haze-spot degeneracies can be lifted with sufficiently high-quality transmission spectra (see their Fig. 22). Therefore, while the finer details of the stellar contamination model are less crucial if one only wishes to constrain the planetary atmosphere, a simple stellar contamination model (e.g. one heterogeneity, as considered in studies such as Pinhas et al. 2018 and Rathcke et al. 2021) is vital to obtain reliable atmospheric inferences for exoplanets transiting active stars.

4.6 A comparison between HST+Spitzer and JWST

We additionally ran retrievals on our newly reduced HST and Spitzer observations to contextualize our JWST results. We consider two Poseidon retrievals, a model without stellar contamination and with a single-unocculted stellar heterogeneity, adopting the same retrieval configuration and priors as used in the JWST analysis (see subsection 4.2 and Table 3). The only difference compared to our previous approach is an altered model wavelength grid (0.9–5.3 μm at R = 20 000) to cover the Spitzer data at longer infrared wavelengths.

Fig. 12 compares our retrieved spectrum and atmospheric constraints from HST and Spitzer to our previously presented JWST results. We focus on the abundances of H₂O, CO₂, and the cloud pressure for this comparison since these are the most well-constrained model parameters from our JWST/NIRISS data (see subsection 4.4). Both retrieval models infer an attenuated absorption feature in the HST/WFC3 data, suggesting a degree of cloud opacity (log P_cloud ≲ −2 to 1σ). The model without stellar contamination attributes the absorption near 1.1 and 1.4 μm to H₂O, with an abundance constraint at the order-of-magnitude level (⁠|$\log X_{\mathrm{H_2 O}} = -2.43_{-1.05}^{+0.87}$|⁠). However, the unocculted star-spot model, which obtains an improved fit to the Spitzer data, finds the HST/WFC3 data can be explained by either H₂O absorption or CH₄ in combination with star-spots. This three-way degeneracy between f_het, |$\log X_{\mathrm{H_2 O}}$|⁠, and |$\log X_{\mathrm{CH_4}}$| – caused by the lack of data shortwards of 1.0 μm – removes the H₂O detection one could claim under the assumption of no stellar contamination. Therefore, HAT-P-18 b’s H₂O abundance is essentially unconstrained from the HST and Spitzer observations. In contrast, our JWST/NIRISS observation has the wavelength coverage to break the degeneracies between unocculted star-spots and the atmospheric composition, affording detailed characterization of HAT-P-18 b’s atmosphere.

5 SUMMARY AND DISCUSSION

This work has led to important inferences on both the stellar activity of HAT-P-18 and the atmosphere of its hot Saturn. Our main results are as follows:

• The transit light curves show evidence of a spot-crossing event (> 5σ); we inferred that the most likely parameters for that occulted spot are a position on the projected stellar surface of (x, y) = (0.090 ± 0.005, 0.42 ± 0.05) R_*, with a radius of 0.116 ± 0.014 R_*, and a temperature colder than the star of ΔT = –93 ± 15 K.

• The main features in the transmission spectrum retrieved from NIRISS/SOSS observation are multiple absorption features produced by water vapour (12.5σ) with a retrieved abundance of log  H₂O ≈ −4.4 ± 0.3, a rise at redder wavelengths due to a CO₂ absorption feature (7.3σ) and evidence of Na (2.7σ). Also, there is a slope towards bluer wavelengths caused by unocculted star-spots (5.8σ) and a uniform, grey cloud deck (7.4σ). The latter also mutes some absorption features.

• Modelling stellar heterogeneities led to four slightly different solutions for the occulted spot, and we showed that different treatments of the unocculted active regions could be used. Fortunately, we found that the different solutions for the active regions had no significant impact on the retrieved transmission spectrum and atmospheric properties of HAT-P-18 b.

• Modelling spot spectra is best achieved with stellar models with lower surface gravities than the stellar photosphere. For the most likely solution of the occulted spot, that is a Δlog g = 1.16 ± 0.19 dex, which is in perfect agreement with our solutions for the unocculted spots with free log g.

We proceed to discuss the implications of our results and their potential impact on exoplanetary atmosphere studies in transmission spectroscopy.

5.1 The atmosphere of HAT-P-18 b in context

The best-fitting atmosphere model confirms the presence of H₂O and clouds in the terminator region of HAT-P-18 b, which is in agreement with previous work (Tsiaras et al. 2018). As in Fu et al. (2022), our transmission spectrum shows a slope towards bluer wavelengths; however, to explain it, instead of unocculted stellar heterogeneities, these authors inferred Rayleigh haze scattering, as also suggested by Kirk et al. (2017) from a ground-based observation. The TLSE, which was not taken into consideration in those studies, is likely responsible for that difference, as well as the 10× higher water abundance in Fu et al. (2022). We also find 10× less CO₂ than was reported by Fu et al. (2022), and we furthermore cannot confirm their detection of CH₄; both differences could come from the data reduction and modelling differences.

We show in Fig. 13 our retrieved chemical abundances in comparison to an equilibrium chemistry model from FastChem 2 (Stock, Kitzmann & Patzer 2022) at 600 K with a metallicity of 0.1 dex and a solar C/O ratio. We selected 10× sub-solar metallicity for this reference equilibrium model based on our retrieved sub-solar H₂O abundance (which is our most precisely constrained abundance), allowing us to compare the other retrieved molecular abundances to equilibrium expectations. The retrieved Na abundance is consistent with this equilibrium model, as are the non-detections of CO, HCN, and NH₃. On the other hand, our retrieved abundance of CO₂ is significantly more abundant, whereas CH₄ and K are depleted by at least two orders of magnitude. We also computed equilibrium chemistry models with different C/O ratios, but the abundances did not change significantly. The differences between our free-retrieval results and the expectations from equilibrium chemistry may be indicative of either 1) non-equilibrium effects, like photochemistry for CH₄ and clouds for K or 2) a demonstration of the sensitivity of free-retrievals to individual points which could be driving an anomalously high CO₂ abundance. Methane, which would be expected at the temperature of HAT-P-18 b, was also not detected in the atmosphere of the slightly warmer hot-Saturn WASP-39 b (Feinstein et al. 2023). The C/O ratio is still uncertain for HAT-P-18 b’s atmosphere and the need for longer wavelength observations with other instruments, such as NIRSpec G395H, are required to measure all the main carbon-bearing species and further investigate possible CH₄ signatures in the infrared.

Our detection of CO₂ may be spurious, and the abundance estimate unreliable. An edge effect in the NIRISS observation may be the culprit, as the CO₂ band is right at the red end of the SOSS spectrum, even exceeding it slightly. Longer wavelength data with the Near-Infrared Spectrograph (NIRSpec) G395H grating would allow the detection of CO₂ to be confirmed and its abundance refined.

Our atmospheric analyses consistently find a uniform cloud deck as the only aerosol required to explain HAT-P-18 b’s NIRISS/SOSS transmission spectrum. Atmospheric studies of hot Jupiters often predict a degree of patchy clouds in hot Jupiter atmospheres (e.g. Parmentier et al. 2016; Komacek et al. 2022). However, despite allowing for patchy clouds, our retrievals all favour terminator cloud fractions |$\gt 83~{{\ \rm per\ cent}}$| (to 2σ). Future general circulation model studies focused on HAT-P-18 b can further investigate whether uniform cloud coverage is more favoured for this Saturn-mass planet with a colder equilibrium temperature (∼ 850 K) than many hot Jupiters. Additional observations in the infrared will also be informative.

5.2 Legacy of Hubble and Spitzer in light of JWST results

For our HST and Spitzer transmission spectrum, the model without stellar contamination finds an abundance of H₂O (log H₂O = |$-2.43_{-1.05}^{+0.87}$|⁠) and a cloud pressure (log P_cloud ≲ −2 ≈ 10 mbar) in good agreement with the findings of Tsiaras et al. (2018) (log H₂O = −2.63 ± 1.18 and log P_cloud = −2.18 ± 0.91 ≈ 6.6 mbar). However, we showed that a model considering the stellar contamination with a single heterogeneity removes the H₂O detection because of a three-way degeneracy between the coverage fraction of heterogeneities and the abundances of H₂O and CH₄ due to the lack of the visible spectrum with HST/WFC3.

The impact of spots on optical transmission spectra has been known for some time to mimic a Rayleigh scattering slope (e.g. McCullough et al. 2014; Rackham, Apai & Giampapa 2018). Also, for exoplanets orbiting red dwarf stars, the 1.4 μm water absorption feature observed with HST/WFC3 was sometimes debated to come from potential water molecules in their cold star-spots (e.g. Barclay et al. 2021). Now, with NIRISS/SOSS, we can break the degeneracy between a hazy atmosphere and star-spots, and potentially for colder stars, we could confirm if water absorption features come from star-spots or the exoplanetary atmosphere.

HST and Spitzer inference studies have shown that water vapour, alkali metals, as well as clouds are common in hot giant planets (e.g. Madhusudhan 2019). Still, there often remained a degeneracy between water and clouds, preventing robustly retrieving the abundance of H₂O and the location of the cloud deck (e.g. Benneke & Seager 2012, 2013). In this work, we confirmed these detections of water vapour, clouds and potentially sodium in a hot giant atmosphere, and we demonstrate that the improved wavelength coverage provided by NIRISS/SOSS is sufficient to break the cloud-metallicity degeneracy.

5.3 Challenges of occulted spot analysis in the era of JWST

Our work adds to the growing body of literature demonstrating that the inference of the transit and occulted spots properties with JWST observations can and should be done to avoid possible biases in atmospheric estimates (e.g. Bixel et al. 2019; Rackham et al. 2023). A semi-analytic tool like spotrod (Béky, Kipping & Holman 2014), used here, is less computationally expensive compared to other tools that use a pixelation approach to model the stellar disc. Parallelization of those tools could be a simple solution to further improve the computation time.

Moreover, we found that modelling active features from light curves, with the current tools and knowledge, is a degenerate problem. The standard practice of fixing the orbital parameters using a white light-curve fit, which also fixes the occulted spot size and position, may prevent ruling out some families of solutions. Therefore, we are currently developing a package to simultaneously fit all spectral channels using both wavelength-dependent and wavelength-independent parameters. Still, theoretical advances are needed to understand the limits of inferences from occultations and break the degeneracy between star-spot size and temperature. Fortunately, those degeneracies do not impact the retrieved transmission spectrum significantly. Furthermore, spotrod assumes that active features have a circular shape and the same limb-darkening law as the star. This can also impact the retrieved properties and make it more difficult to lift some degeneracies. The occulted feature may, moreover, have a more complex structure consisting of a cool spot and a hot facula, which could explain the smaller temperature difference retrieved for the occulted spot compared to the unocculted spots. Nonetheless, the retrieved stellar fraction of the occulted spot |$f_{\mathrm{oc.\, spot}} = 1.35 \pm 0.02~{{\ \rm per\ cent}}$| is coherent with our preferred solution for the coverage fraction of unocculted spots |$f_{\mathrm{unoc.\, spot}} = 12^{+5}_{-3}$| per cent.

5.4 Accounting for unocculted heterogeneities in JWST transmission spectra

Our best-fitting model suggests that the slope towards bluer wavelengths in this NIRISS/SOSS observation comes mainly from unocculted stellar heterogeneities. Our work shows that we should constrain planetary atmosphere and stellar contamination properties simultaneously with JWST NIRISS/SOSS data; otherwise, there is a risk of biasing atmospheric inferences. However, we could not distinguish clearly which of the four stellar contamination treatments investigated (one and two heterogeneities, one and two heterogeneities with free surface gravity) is the most plausible, but we were able to show that a particular choice of treatment does not significantly impact the atmospheric inferences.

Furthermore, the spectra of active regions are currently modelled with 1D stellar models, which is another limitation, particularly for faculae (Norris et al. 2017; Johnson et al. 2021), because they fail to capture aspects of facular contrasts, such as the centre-to-limb variation (e.g. Yeo, Solanki & Krivova 2013) and the dependence on stellar metallicity (Witzke et al. 2018). More theoretical work, such as magnetohydrodynamic simulations of magnetic features, is needed to understand better and, thus, model the contrast of active regions more accurately as a function of stellar fundamental parameters and activity. In the meantime, we suggest that spots’ flux spectra be modelled with stellar models with lower surface gravities than the star.

6 CONCLUSION

We presented atmospheric retrieval analyses of two transmission spectra of HAT-P-18 b obtained with JWST NIRISS/SOSS, and HST/WFC3 with Spitzer/IRAC. We also modelled a spot-crossing event in the transit observed with NIRISS/SOSS. We confirmed that the wavelength coverage and spectral resolution of NIRISS/SOSS considerably improve the detection significance and the abundance constraints on chemical species compared to observations from HST and Spitzer combined. In including stellar heterogeneities in transit fits and atmospheric retrievals, we implemented new model considerations designed to fit the local surface gravity of stellar heterogeneities. This work is instructive for the upcoming JWST transit observations, particularly by informing the community on disentangling stellar and planetary atmosphere signals.

7 SOFTWARE

• astropy; Astropy Collaboration et al. (2013, 2018)

• batman; Kreidberg (2015)

• dynesty; Speagle (2020)

• emcee; Foreman-Mackey et al. (2013)

• ExoTiC-LD; Wakeford & Grant (2022)

• i python; Pérez & Granger (2007)

• juliet; Espinoza, Kossakowski & Brahm (2019b)

• matplotlib; Hunter (2007)

• nestle; Skilling (2004)

• numpy; Harris et al. (2020)

• Poseidon; MacDonald (2023)

• P yMSG; Townsend & Lopez (2023 )

• pymultinest; Buchner et al. (2014)

• scipy; Virtanen et al. (2020)

• spotrod; Béky, Kipping & Holman (2014)

• supreme-SPOON; Radica et al. (2023)

ACKNOWLEDGEMENTS

We thank Steven Rogowski for his initial work on HST and Spitzer data. This work is based on observations made with JWST, as well as the Hubble and Spitzer Space Telescopes. This project is undertaken with the financial support of the Canadian Space Agency. MFT acknowledges financial support from the Fonds de Recherche du Québec – Nature et technologies (FRQNT) and funding from the Trottier Family Foundation in their support of iREx. RJM acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant number HST-HF2-51513.001, awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. MR acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) as well as from FRQNT and iREx. CP acknowledges support from FRQNT and TEPS Ph.D. scholarships and the NSERC Vanier scholarship. KM acknowledges financial support from the FRQNT. OL acknowledges financial support from FRQNT and iREx. DJ is supported by NRC Canada and by an NSERC Discovery Grant. JDT was supported for this work by NASA through the NASA Hubble Fellowship grant number HST-HF2-51495.001-A awarded by the STScI, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract number NAS5-26555.

8 DATA AVAILABILITY

All data used in this study are publicly available from the Barbara A. MAST⁷ and the Spitzer Heritage Archive.⁸

Footnotes

References

APPENDIX A: ADDITIONAL REDUCTION

We carried out an independent reduction on the HAT-P-18 b SOSS TSO using NAMELESS (described in depth in Coulombe et al. 2023). We applied all stage 1 steps of the jwst pipeline except for the dark current subtraction. We subsequently applied the world coordinate system, source type, and flat field steps of stage 2. We proceeded with the background subtraction, scaling independently the two regions separated by the jump of the model background provided by STScI, as described in detail in Lim et al. (2023). Treatment of the 1/f noise was performed by scaling each individual column of the trace independently and finding the value that results in the lowest Chi-square for a given column and integration (Coulombe et al. 2023). Finally, we extracted the 2D spectra using the getSimpleSpectrum function of the transitspectroscopy pipeline⁹ with an aperture width of 30 pixels.

We perform the light-curve fitting on the extracted spectrophotometric observation in an identical manner to the one described in Section 3. We fix the orbital parameters (T₀, b, a/R_*) to the best-fitting values from the order 1 white light-curve fit of the supreme-SPOON pipeline to ensure consistency. A comparison of the best-fitting white light-curve parameters for NIRISS/SOSS with the two different pipelines is available in Table A1. We note that the parameter values of R_p/R_*, b and a/R_* differ by 1.2–1.3σ. These parameters are correlated (e.g. see the corner plot for one of the broad-band light-curve fits; Fig. B4), and each reduction may simply prefer a slightly different place in the same family of solutions with a comparable agreement. The retrieved transmission spectrum, along with the one from the reference supreme-SPOON reduction, are shown in Fig. A1. These transmission spectra obtained through two different pipelines are in overall good agreement; showing consistent transit depths and features. However, the amplitude of the CO₂ and Na features are different between the two reductions (similarly to, e.g. Radica et al. 2023) and would lead to some differences in the retrieved abundances.

APPENDIX B: ADDITIONAL MATERIALS

Author notes