Solar-to-supersolar sodium and oxygen absolute abundances for a ‘hot Saturn’ orbiting a metal-rich star

ABSTRACT

We present new analysis of infrared transmission spectroscopy of the cloud-free hot-Saturn WASP-96b performed with the Hubble and Spitzer Space Telescopes (HST and Spitzer). The WASP-96b spectrum exhibits the absorption feature from water in excellent agreement with synthetic spectra computed assuming a cloud-free atmosphere. The HST-Spitzer spectrum is coupled with Very Large Telescope (VLT) optical transmission spectroscopy which reveals the full pressure-broadened profile of the sodium absorption feature and enables the derivation of absolute abundances. We confirm and correct for a spectral offset of |$\Delta R_{{\rm p}}/R_{\ast }=(-4.29^{+0.31}_{-0.37})\, \times 10^{-3}$| of the VLT data relative to the HST-Spitzer spectrum. This offset can be explained by the assumed radius for the common-mode correction of the VLT spectra, which is a well-known feature of ground-based transmission spectroscopy. We find evidence for a lack of chromospheric and photometric activity of the host star which therefore make a negligible contribution to the offset. We measure abundances for Na and O that are consistent with solar to supersolar, with abundances relative to solar values of |$21^{+27}_{-14}$| and |$7^{+11}_{-4}$|⁠, respectively. We complement the transmission spectrum with new thermal emission constraints from Spitzer observations at 3.6 and 4.5 |$\mu$|m, which are best explained by the spectrum of an atmosphere with a temperature decreasing with altitude. A fit to the spectrum assuming an isothermal blackbody atmosphere constrains the dayside temperature to be T_p = 1545 ± 90 K.

techniques: spectroscopic, planets and satellites: composition, planets and satellites: individual: WASP-96b, planets and satellites: atmospheres, stars: abundances, stars: individual: WASP-96

1 INTRODUCTION

From massive gas-giants down to Earth-mass worlds, from ultra-hot planets to much cooler and potentially habitable Earth-like analogues, exoplanets have now been found orbiting many of the stars across the solar neighbourhood. Following nearly three decades of exoplanet detections, the exoplanet census now shows that planet formation is ubiquitous with Super-Earths being common and hot Jupiters, planets similar to Jupiter in the Solar system but with orbital periods shorter than about 10 d, uncommon (Howard et al. 2010, 2012; Mayor et al. 2011; Burke et al. 2015; Wang et al. 2015; Brewer et al. 2016). Due to their large radii and hot temperatures, hot Jupiters host the only class of exoplanet atmosphere that can currently be observationally characterized in detail with spectroscopic observations over an entire orbital phase (Fortney, Dawson & Komacek 2021). Even the first steps to characterize exoplanets via transits have started to reveal astoundingly diverse worlds very different from the planets of the Solar system with: (i) evaporating and escaping planetary atmospheres resembling comet-like tails (Ehrenreich et al. 2015); (ii) atmospheric hazes and clouds composed of minerals and metal oxides, which tend to be ubiquitous (Lecavelier Des Etangs et al. 2008; Huitson et al. 2013; Sing et al. 2013; Nikolov et al. 2015; Alam et al. 2020; Carter et al. 2020; Wilson et al. 2020; Spyratos et al. 2021); and (iii) a huge diversity and a continuum from clear to cloudy atmospheres (Sing et al. 2016; Gibson et al. 2017; Evans et al. 2018; Wakeford et al. 2018; Alam et al. 2021; Fu et al. 2021; Sheppard et al. 2021; Spake et al. 2021). From all types, exoplanet atmospheres free of clouds offer a unique opportunity to constrain elemental abundances (e.g. Na, K, water, and carbon-bearing molecules and metallicity objectively), as clouds and hazes obscure the amplitudes of absorption features (Fortney 2005; Fortney et al. 2010). Such rare exoplanets enable unbiased absolute abundance constraints, which may inform theory of irradiated gas-giant exoplanet atmospheres and giant planet formation (Öberg, Murray-Clay & Bergin 2011; Madhusudhan 2012; Madhusudhan et al. 2017). Because the composition of protoplanetary discs vary significantly with distance from the star, an exoplanet’s atmospheric composition may contain a fossil record of the conditions of the disc at the location where the planet formed. Madhusudhan (2012) showed that the C/O ratio carries important information on the chemistry of volatile species, which has started to provide definitve constraints (Line et al. 2021). Lothringer et al. (2021) explored the insight that the measurement of refractory abundances can provide into a planet’s origins. Through refractory-to-volatile elemental abundance ratios, they demonstrated that a planet’s atmospheric rock-to-ice fraction can be estimated, which can constrain planet formation and migration scenarios.

Motivated by the recent observational results from space and the ground, theoretical predictions for strong alkali features, stratosphere inducing TiO and VO, and diversity of clouds hazes and dust, we initiated a large transmission spectral survey of 20 exoplanets with the FORS spectrograph on the Very Large Telescope (VLT). Optical transmission spectra from this program have been reported for WASP-96b, WASP-103b, WASP-88b, and WASP-110b (Nikolov et al. 2018b, 2021; Wilson et al. 2020; Spyratos et al. 2021). In this study, we performed follow-up reconnaissance observations of WASP-96b using Hubble Space Telescope (HST)and Spitzer, aiming at detecting water in the transmission spectrum of the planet.

WASP-96b is a ‘hot Saturn’ with a planetary mass M_p = 0.48 ± 0.03M_J, where M_J is the mass of Jupiter, planet radius R_p = 1.20 ± 0.06R_J, where R_J is the radius of Jupiter and an equilibrium temperature T_eq = 1285 ± 40 K on a 3.4-d circular orbit around a V = 12.2 chromospherically quiet and photometrically stable (Section 3.4) G8 star at a distance of 356 ± 5 pc in the southern constellation Phoenix (Hellier et al. 2014). Unlike the optical transmission spectrum for a majority of hot Jupiters known to date, the spectrum of WASP-96b shows very clear pressure-broadened wings of the sodium line, evidence of potassium, and a near-UV Rayleigh scattering slope, with the latter defining the hydrogen continuum level and proving the planet has a clear atmosphere at the limb. The VLT spectrum provides an absolute sodium abundance, 2–18 × the solar value (1σ) and an atmospheric metallicity in agreement with the mass–metallicity trend observed for Solar-system planets and exoplanets (Thorngren et al. 2016; Nikolov et al. 2018b). This suggests that the absorption signature of molecules, such as water with absorption bands at ∼0.95, 1.15, and 1.4 |$\mu$|m, carbon monoxide, and methane at 3.6 and 4.5 |$\mu$|m are expected to also be free from the suppressing effect of clouds and hazes in the infrared (Sing et al. 2016). We carried out (i) transit observations with the HSTin spatial scanning mode, and (ii) Spitzer Space Telescope (Spitzer) transit and eclipse observations, to extend the planet’s optical transmission spectrum at these wavelengths, constrain the atmospheric composition and day side brightness temperature of the planet (GO-15469 and 14255; Nikolov et al. 2018c, 2019). An analysis of the transit observations has been reported by Yip et al. (2021), finding absorption from water and an offset between the optical and infrared.

In this paper, we present new results for the infrared transmission spectrum, which overall agree with the results from Yip et al. (2021), finding solar to supersolar abundances for sodium and oxygen. We additionally find the planet’s potassium and carbon abundances to be low with broad uncertainties owing to lower precision of the spectrum at the relevant wavelengths. We present new results for the elemental abundances of the host star, its chromospheric activity index, constrain the dayside temperature of the planet, and refine the system parameters and planet orbital ephemeris. We demonstrate in this work, as has been detailed in Nikolov et al. (2018b), that the observed offset of the optical spectrum is explained by the assumed absolute planet-to-star radius ratio for the calculation of common mode (wavelength independent) systematics model for the ground-based spectra. As the hot Jupiter with the clearest known atmosphere, WASP-96b is poised to become an important prototype exoplanet target for the James Webb Space Telescope and provide the first accurate abundance measurements, unbiased by clouds and hazes.

We describe the observations in Section 2, discuss the data reduction and light-curve analysis in Section 3, present results from a retrieval analysis in Section 4, discuss our findings in the context of literature result in Section 5, and summarize and conclude in Section 6.

2 OBSERVATIONS AND REDUCTIONS

2.1 HST WFC3 spatial scan mode transit spectroscopy

We observed two primary transits of WASP-96b with the Wide Field Camera 3 (WFC3) aboard the HSTas part of Program 15469 (Nikolov et al. 2018c) in spatial scan mode. Both visits included five consecutive HST orbits, with the first, second, third, and fifth orbits covering out-of-transit phases and the fourth orbit during the transit. The HST orbital sequence was planned such that the in-transit orbit covers the planet phase between second and third contacts (end of ingress and start of egress, respectively) to maximize the precision of the measured planet radius (Figs A1 and A2). Each observation was performed in spatial scan mode using the 256 × 256 subarray centred on the target to reduce read-out overheads. The spatial scan mode spreads the stellar flux across more pixels compared to a fixed pointing (stared) observation, allowing for a longer exposure times with higher duty cycle. Prior to the spectral observations, we acquired direct images of the target to facilitate a wavelength solution and to inspect the vicinity of the target for field stars that can contaminate the spectra. We exposed for 1.111 and 1.389s using the F127M and F139M filters for the G102 and G141 grisms, respectively. An inspection of the direct images showed the field around WASP-96 to be free of stars.

The first observation covers the transit on UT 2018 December 18. We used the dispersive element G102, which covers the spectral range from 0.8 to 1.15 |$\mu$|m at a resolving power of R ∼ 210 (at 1 |$\mu$|m) and a dispersion of 2.45 nm px⁻¹. We used forward scanning with a rate of 0.013 arcsec s⁻¹ and the SPARS25 sampling sequence with 11 iterations per exposure (NSAMP9) resulting in total integration times of ∼179 s and scans across 2.3 arcsec (∼18 px). Typical count levels reached a maximum of ∼2.9 × 10⁴ analogue-to-digital units (ADU), i.e. well within the linear regime of the detector.

The second observation covered the transit on UT 2018 December 28. We exploited the dispersive element G141, which covers the spectral range from 1.075 to 1.7 |$\mu$|m at a resolving power of R ∼ 130 (at 1.4 |$\mu$|m) and a dispersion of 4.65 nm px⁻¹. We used forward scanning with a rate of 0.022 arcsec s⁻¹ and the SPARS25 sampling sequence with 12 iterations per exposure (NSAMP8) resulting in total integration times of ∼157 s and scans across 3.4 arcsec (∼268 px). The spectra reached typical maximum count levels of ∼2.8 × 10⁴ ADU, i.e. within the linear regime of the detector.

We performed data reduction and analysis following the procedures outlined in Nikolov et al. (2018a). Our analysis started with the 2D ima spectra, produced by the calwf3 pipeline (v3.1), which are corrected for bias, dark current, flat-field, and detector non-linearity effects. We extracted flux for WASP-96b from each exposure by taking each successive non-destructive read difference. We removed the background on each individual read difference by taking the median flux in a box away from the stellar spectrum. We then determined the flux weighted centre of the WASP-96 scan and set to zero all pixel values located more than 20 pixels above and below along the cross-dispersion axis. Application of this top hat filter had the effect of eliminating most of the pixels affected by cosmic rays. The final reconstructed images were produced by adding together the read differences for each exposure. Any remaining cosmic rays were identified and corrected, following the procedure described in Nikolov et al. (2014, 2018a). We find a total of 2 or 3 pixels that sample the spectra to be affected by cosmic rays for each reconstructed image.

We performed a spectral extraction using a fixed-size box by summing the flux of all pixels. The box had dimensions of 171 × 28 and 169 × 39 pixels for the G102 and G141 data, respectively, and centred for each individual exposure. To identify the box positions along the dispersion and cross-dispersion axis we took the flux-weighted mean of each 2D spectrum. The target drift along the dispersion and cross-dispersion axis are shown in Fig. A1.

We established wavelength solutions for each grism time series by cross-correlating each spectrum with an ATLAS synthetic stellar spectrum (Kurucz 1979) that matches the astrophysical properties of the WASP-96 host star i.e. |${\mathrm{T_{\rm {eff}}=5500\pm 150}}$| K, log g = 4.42 ± 0.02, and [Fe/H] = 0.14 ± 0.19, reported in Hellier et al. (2014). We obtained synthetic spectra from the HST exposure time calculator, which also accounts for the instrument sensitivity.

2.2 Spitzer IRAC transit and eclipse photometry

Photometric data were collected during two primary transits and two secondary eclipses using the Spitzer Space Telescope (Spitzer, Werner et al. 2004) Infrared Array Camera (IRAC, Fazio et al. 2004), as part of Program 14255, Nikolov et al. (2019). Each IRAC observation covers equal duration of in-transit, pre-ingress, and post-egress time series (Figs A3 and A4). Each exposure consisted of 64 subarray frames of 32 × 32 pixels. We employed an exposure time of 1.92 s and obtained a total of 9600 images for each observation. Details of the IRAC observations are summarized in Table 1.

Table 1.

Spitzer IRAC observations.

Channel	Start UTC	End UTC	Event
ch1/3.6	2019-10-22 10:29:59	2019-10-22 15:51:20	transit
ch2/4.5	2019-10-29 06:59:07	2019-10-29 12:20:28	transit
ch1/3.6	2019-10-30 23:54:54	2019-10-31 05:16:15	eclipse
ch2/4.5	2019-11-06 20:15:01	2019-11-07 01:36:22	eclipse

Our data reduction and analysis procedures follow the methodology detailed in Nikolov et al. (2015). The data were calibrated by the Spitzer pipeline version 21 and are available in the form of Basic Calibrated Data (.bcd) files. After organizing the data, we converted the images from flux in mega-Jansky per steradian (MJy sr⁻¹) to photon counts i.e. electrons using the information provided in the fits headers. We then performed an outlier filtering for hot (energetic) or lower pixels in the data by following each pixel through time. This task was performed in two steps, first flagging all pixels with intensity more than 8σ compared to the median value computed from the five preceding and five following images. The values of these flagged pixels were replaced with the local median value. In the second pass, we flagged and replaced outliers above the 4σ level, following the same procedure. The total fraction of corrected pixels was ∼0.4 per cent for the 3.6 |$\mu$|m and ∼0.1 per cent for the 4.5 |$\mu$|m channel.

We then subtracted the background flux from the time series by performing an iterative 3σ outlier clipping. For each image we removed the pixels within the stellar point spread function (PSF), background stars, or hot pixels. We then created a histogram from the remaining pixels and fitted a Gaussian to determine the sky background. Prior to aperture photometry, we measured the target’s PSF position on each image using the flux-weighted centroiding method with a circular region with a radius of 3 pixels. The variation of the x and y positions of the PSF on the detectors were measured to be ∼0.2 for the 3.6 |$\mu$|m and ∼0.1 pixels for the 4.5 |$\mu$|m transits. We measure variation of the x and y positions of the PSF of ∼0.5 and ∼0.1 for the 3.6 |$\mu$|m and ∼0.3 and ∼0.4 for the 4.5 |$\mu$|m eclipses.

We performed photometry with both fixed and time variable apertures and filtered the resulting light curves for 5σ outliers with a width of 20 data points. In the first method we covered the range of fixed (in time) radii from 1.5 to 3.5 pixels in increments of 0.1 pixels. In the second approach, we scaled the size of the extraction aperture by the value of a quantity known as the noise pixel parameter (Mighell 2005; Knutson et al. 2012). The best results from both methods were identified by examining both the residual root mean square (rms) after fitting the light curves from each channel, as well as the white and red noise components measured with the wavelet technique detailed in Carter & Winn (2009). While we find no evidence for a particular aperture size that produces significantly lower light curve scatter compared to all apertures, the fitted transit depths are consistent within their 1σ uncertainties. The time variable approach produced the lowest light curve scatter for the 3.6 |$\mu$|m observations with an aperture radii of 1.97 and 2.05 pixels with sky annuli from 2.10 to 5.72 and from 2.05 to 5.54 for the transit and eclipse, respectively. The lowest light curve scatter for the 4.5 |$\mu$|m is found from the fixed and time variable aperture methods for the transit and eclipse, respectively. We find aperture radii of 2.25 and 2.43 pixels and sky annuli defined between radii 2.25 and 6.08 pixels and from 2.43 to 6.55 pixels for the transit and eclipse, respectively (Figs A3 and A4).

2.3 TESS transit photometry

WASP-96 was observed by the Transiting Exoplanet Survey Satellite (TESS Ricker et al. 2014) on Camera 2 CCD 1 during the Prime and Extended Mission between 2018 August 23 (Sector 2) and 2018 September 20, and between 2020 August 26 and 2020 September 21 (Sector 29) with a 2-min cadence in the full-frame images (FFIs). We obtained light curves for both observing campaigns from the MAST archive. The light curves were produced by the TESS Science Processing Operations Center (SPOC) pipeline, which was used to calibrate FFIs and to assign world-coordinate system information to the FFI data delivered to the MAST (Caldwell et al. 2020). First, we inspected the simple aperture photometry (SAP) fluxes and identified seven transits in each campaign. We chose to perform light-curve fits with the goal to identify and remove residual systematic effects. We obtained portions of the light curves centred on the primary transits with 6 h out-of-transit baseline data prior ingress and additional 6 h after egress. One transit from the second campaign has been discarded in our analysis due to a significantly higher flux scatter as compared to the rest of the data.

2.4 ASSASN phase curve photometry

WASP-96 was monitored for photometric variability by the Ohio State University’s All-sky Automated Survey for Supernovae (ASAS-SN) Photometry Database (Shappee et al. 2014; Jayasinghe et al. 2019). A total of 233 photometric measurements were collected and processed during five consecutive observing campaigns starting on 2014 May 12 and ending on 2018 September 24. Each campaign, except the last one, covers the full visibility window of WASP-96 from May to January.

2.5 ESO/MPG 2.2 FEROS high-resolution spectroscopy

The host star of WASP-96b was observed with the MPG/ESO 2.2-metre telescope equipped with the FEROS high-resolution spectrograph (ESO program: 098.A-9007, PI P. Sarkis). Two observations with an exposure time of 900 s were obtained on UT 2016 December 12 and 19, using a fibre with a diameter of 2 arcsec and a spectral resolution of R = λ/δλ = 48 000. The spectra were reduced with ESO’s FEROS pipeline to phase three products i.e. extracted one-dimensional spectra with wavelength solution. Both spectra have a signal-to-noise ratio of 27 over a wavelength bin of 0.03Å. A magnified view centred on the Ca ii H&K lines from the combined spectra is shown in Fig. B1.

3 LIGHT-CURVE ANALYSIS

Our light-curve analysis of the HST, Spitzer, and TESS data follows the methods detailed in Nikolov et al. (2015, 2018a, 2021).

3.1 HST WFC3 spectrophotometry

3.1.1 White-light curves

We produced white-light curves by summing the flux of each spectrum from 0.78 to 1.13 and from 1.12 to 1.65 |$\mu$|m for the G102 and G141 transit time series, respectively (Fig. A2). Both data sets exhibit the well-known ‘hook’ systematic, which correlates with the HST orbital phase and is considered to originate from charge-trapping in the WFC3 detector (Deming et al. 2013; Huitson et al. 2013; Zhou et al. 2017). The out-of-transit flux also exhibits a longer term drift, which is approximately linear in time.

We fit the transit and systematic effects of the white-light curves by treating the data as a Gaussian process¹ (Gibson et al. 2012). The transit parameters orbital period, eccentricity, inclination i, and normalized semimajor axis a/R_* (where R_* is the radius of the star) were held fixed to the previously determined values of Hellier et al. (2014) and Nikolov et al. (2018b), which are detailed in Table 2. We accounted for the stellar limb darkening by adopting the four-parameter non-linear limb-darkening law (Claret 2000) and computed the values of the coefficients (c₁, c₂, c₃, c₄) using a three-dimensional stellar atmosphere grid (Magic, Weiss & Asplund 2015). In these calculations, we adopted the closest match to the effective temperature, surface gravity, and metallicity of the exoplanet host star found in Hellier et al. (2014). The mid-time T_mid and planet-to-star radius ratio R_p/R_* were allowed to vary in the fit to each of the two white-light curves.

Table 2.

System parameters from HST WFC3.

Parameter, symbol, unit	Value
Orbital period, P, (d)	3.4252602 (fixed)
Eccentricity, e	0 (fixed)
Inclination, i, (°)	85.14 (fixed)
a/R*	8.84 (fixed)
|${\it HST}\, G102$|	0.8–1.15 |$\mu$|m
Mid-transit time, Tmid, JD	|$2458470.7795^{+0.0037}_{-0.0015}$|
Planet-to-star radius ratio, |$\mathrm{R_{\rm {p}}/R_{\ast }}$|	|$0.1186^{+0.0015}_{-0.0013}$|
First limb-darkening coefficient, c1	0.631621 (fixed)
Second limb-darkening coefficient, c2	−0.251583 (fixed)
Third limb-darkening coefficient, c3	0.417851 (fixed)
Fourth limb-darkening coefficient, c4	−0.190682 (fixed)
Characteristic amplitude, A (ppm)	|$2226^{+5369}_{-991}$|
Characteristic length-scale, ηϕ	4 ± 3
Characteristic length-scale, ηx	|$7^{+4}_{-5}$|
Characteristic length-scale, ηy	|$2^{+5}_{-2}$|
Intercept, a0	|$0.9995^{+0.0016}_{-0.0018}$|
Slope, a1	|$0.00049^{+0.00048}_{-0.00034}$|
Re-scaled uncertainty, σw, (ppm)	|$158^{+78}_{-64}$|
|${\it HST}\, G141$|	1.075–1.700 |$\mu$|m
Mid-transit time,Tmid JD	|$2458481.05626^{+ 0.00041}_{-0.00046}$|
Planet-to-star radius ratio, Rp/R*	|$0.11962^{+0.00029}_{-0.00037}$|
First limb-darkening coefficient, c1	0.657415 (fixed)
Second limb-darkening coefficient, c2	−0.105797 (fixed)
Third limb-darkening coefficient, c3	0.031906 (fixed)
Fourth limb-darkening coefficient, c4	−0.019303 (fixed)
Characteristic amplitude, A (ppm), A (ppm)	|$2038^{+3167}_{-1156}$|
Characteristic length-scale, ηϕ	4 ± 2
Characteristic length-scale, ηx	8 ± 4
Characteristic length-scale, ηy	4 ± 2
Intercept, a0	0.9994 ± 0.002
Slope, a1	|$0.00012^{+0.00028}_{-0.00025}$|
Re-scaled uncertainty, σw, (ppm)	128 ± 24

Under the Gaussian process assumption, the data likelihood is a multivariate normal distribution with a mean function μ describing the deterministic transit signal and a covariance matrix K that accounts for stochastic correlations (i.e. poorly constrained systematics) in the light curves:

$$\begin{eqnarray} p(\boldsymbol {f} \vert \boldsymbol {\theta }, \gamma) = \mathcal {N}(\mu ,K), \end{eqnarray}$$

(1)

where p is the probability density function, |$\boldsymbol {f}$| and |$\boldsymbol {\theta }$| are vectors containing the flux measurements and mean function parameters, respectively, γ is a function containing the covariance parameters, and |$\mathcal {N}$| is a multivariate normal distribution. The mean function μ is defined as

$$\begin{eqnarray} \mu (\boldsymbol {t},\boldsymbol {\hat{t}};a_0,a_1,\boldsymbol {\theta })=[a_0+a_1\boldsymbol {\hat{t}}]\, T(\boldsymbol {\it {t}},\boldsymbol {\theta }), \end{eqnarray}$$

(2)

where |$\boldsymbol {t}$| is a vector of all central exposure time stamps in Julian Date, |$\boldsymbol {\hat{t}}$| is a vector containing all standardized times, that is, with subtracted mean exposure time and divided by the standard deviation, a₀ and a₁ describe a linear baseline trend, |$T(\boldsymbol {\theta })$| is an analytical expression describing the transit, and |$\boldsymbol {\theta }=(i, a/R_{\ast }, T_{\rm {mid}}, R_{p}/R_{\ast }, c_1, c_2, c_3, c_4)$|⁠. We used the analytical formulae of Mandel & Agol (2002) and Kreidberg (2015).

Similar to our earlier HST and VLT studies, we defined the covariance matrix as |$K = \sigma _i^2\delta _{ij} + k_{ij}$|⁠, where σ_i contains the photon noise uncertainties, δ_ij is the Kronecker delta function, and k_ij is a covariance function. The white noise term σ_w was assumed to have the same value for all data points and was allowed to freely vary. We chose to use the Matérn ν = 3/2 kernel with the HST orbital phase (ϕ) and spectral dispersion (x) and cross-dispersion (y) drifts, respectively, as input variables. Our kernel choice is motivated by the study of Gibson et al. (2013), where the Matérn ν = 3/2 kernel is empirically motivated using simulated data, and is the first to use this kernel for light-curve analysis. The GP free parameters were the characteristic correlation amplitude (A) and correlation length-scales for each input parameter (τ_ϕ, τ_x, τ_y). Similar to our previous studies, we also accounted for the long-term out-of-transit baseline trend with a linear time term. As with the linear time term, we also standardized the input parameters before the light-curve fitting. The covariance function was defined as

$$\begin{eqnarray} k_{ij} = A^2(1+\sqrt{3}D_{ij})\exp {-\sqrt{3}D_{ij}}, \end{eqnarray}$$

(3)

where A is the characteristic correlation amplitude and

$$\begin{eqnarray} D_{ij} = \sqrt{\frac{(\boldsymbol {\hat{\phi }}_i-\boldsymbol {\hat{\phi }}_j)^2}{\tau _\phi ^2}+\frac{(\boldsymbol {\hat{x}}_i-\boldsymbol {\hat{x}}_j)^2}{\tau _x^2}+\frac{(\boldsymbol {\hat{y}}_i-\boldsymbol {\hat{y}}_j)^2}{\tau _y^2}}, \end{eqnarray}$$

(4)

where τ_ϕ, τ_x, and τ_y are the correlation length-scale and the hatted variables are standardized. The parameters |$\boldsymbol {X} = (a_0, a_1, T_{\rm {mid}}, R_{p}/R_{\ast }, \sigma _w)$| and |$\boldsymbol {Y} = (A, \tau _\phi , \tau _x, \tau _y)$| were allowed to vary and fixed the orbital period P and eccentricity e to the values reported in Hellier et al. (2014) and the system parameters (i) and (a/R_*) to the values found in Nikolov et al. (2018b). Our choice to fix the two system parameters was determined by the fact that both HST observations lack data points (except one for the last orbit of the G141 visit, Fig. A2) obtained during ingress and egress. We adopted uniform priors for x and log-uniform priors for y. We marginalize the posterior distribution using the Markov Chain Monte Carlo (MCMC) software package emcee (Foreman-Mackey et al. 2013). We identified the maximum-likelihood solution using the Levenberg–Marquardt least-squares algorithm (Markwardt 2009) and initialized three groups of 250 walkers close to that maximum. The first two groups were run for 450 samples and the third one for 3500 samples. To ensure faster convergence, we resampled the positions of the walkers in a narrow space around the position of the best walker from the first run before running for the second group. This helps prevent some of the walkers starting in a low-likelihood area of parameter space, which can require more computational time to converge. Fig. A2 and Table 2 show transit models for each of the two observations computed using the marginalized posterior distributions and the fitted parameters, respectively. We find residual dispersion of 83 and 94 parts-per-million for the blue (G102) and red (G141) light curves, respectively.

3.1.2 Spectroscopic light-curve analyses

We produced spectroscopic light curves by summing the flux of WASP-96 in bands with variable widths in the range from 0.7903 to 1.6396 |$\mu$|m. Our choice for the bin widths was governed by the goal to obtain the highest spectral resolution with spectrophotometric precision that is comparable to the precision of the VLT data. We produced 9 and 56 bins for the G102 and G141 grisms, respectively. The G102 and G141 grisms overlap in wavelength over a narrow region from 1.0195 to 1.1333 |$\mu$|m, which allows a direct comparison of the measured transit depths, Figs 3, 4 and 7.

We established wavelength-independent i.e. common-mode, systematic correction factors for each data set, similar to our VLT and HST studies detailed in Nikolov et al. (2014, 2016). We obtained the correction factors by dividing the white-light curve from each visit by the transit model computed with the parameters detailed in Table 2. The common-mode factors for each visit are shown in Fig. A2.

We fit the spectroscopic light curves using a two-component function that takes into account the transit signal and systematics simultaneously, similar to our studies in Nikolov et al. (2016, 2018b). Prior to fitting for the systematics, we removed the common mode factors from each set of spectroscopic light curves by dividing the raw flux of the spectroscopic light curves by the corresponding common-mode factors for the relevant visit. We computed transit models using the analytical models of Mandel & Agol (2002). In the fits, we allowed only the planet-to-star radius ratio R_p/R_* to vary. We held the four coefficients of the non-linear limb-darkening law to their theoretical values. We obtained these values following the same method as for the white light curves.

We accounted for the systematics using a fourth-order polynomial of HST orbital phase (ϕ) and up to a second-order polynomial of dispersion and cross-dispersion drift, Fig. A1. We produced all possible combinations of systematics models and performed separate fits with each of them included in the two-component function. This approach has been preferred as opposed to GP regression, as the CM-corrected HST spectroscopic light curves exhibit a lower level of systematic effects (Sing et al. 2015, 2016; Nikolov et al. 2018a). We computed the Akaike Information Criterion (AIC) for each fit and estimated the statistical weight of the model depending on the number of degrees of freedom (Akaike 1974). We marginalized the resulting relative radii following Gibson (2014). The highest evidence models included a fourth-order polynomial of HST orbital phase.

3.2 Spitzer IRAC photometry

Our light-curve analysis of the 3.6 and 4.5 |$\mu$|m IRAC transits and eclipses is similar to the approach for the WFC 3 spectroscopic light curves and our previous Spitzer studies (Nikolov et al. 2015; Sing et al. 2016). We analysed the primary/secondary eclipses and instrumental systematics simultaneously. We used the analytic formulae of Mandel & Agol (2002) to fit the transit light curves and assumed the non-linear limb-darkening law with coefficients computed from the three-dimensional stellar atmosphere model grid of Magic et al. (2015). To obtain an analytical eclipse model, we relied on the same formulae, but computed model with limb darkening set to zero. We fixed the remaining system parameters to literature values listed in Table 2 and fitted only for the eclipse depth and central time. To correct for the intrapixel sensitivity-induced flux variations, we fit a polynomial function of the stellar centroid position of the form:

$$\begin{eqnarray} f(t)=a_0 +a_1x+a_2x^2 +a_3y+a_4y^2 +a_5xy+a_6t, \end{eqnarray}$$

(5)

where f(t) is the stellar flux as a function of time, t; x and y are the positions of the stellar centroid on the detector; and a₀ to a₆ are the free parameters of the fit. We produced all possible combinations of terms from equation (5) and performed separate fits. Using the AIC statistic, we marginalized the models to obtain parameters with uncertainties following Gibson (2014).

3.2.1 Primary transits

First, we fit the Spitzer transit light curve with the goal to measure the planet system parameters. We allowed the planet orbital inclination, i, normalized semimajor axis, a/R_*, planet-to-star radius ratio, R_p/R_*, and central transit time, T₀, to vary free in the fit. The planet orbital period and eccentricity were fixed to the values reported in Table 2 and the four limb-darkening coefficients held fixed for the 3.6 and 4.5 |$\mu$|m were c₁ = 0.5493, c₂ = −0.4770, c₃ = 0.4054, c₄ = −0.1443 and c₁ = 0.5480, c₂ = −0.6455, c₃ = 0.6196, c₄ = −0.2232, respectively. We found inclination values of i = 84.75 ± 0.75° and 85.35 ± 0.41° and normalized semimajor axis of a/R_* = 8.48 ± 0.71 and a/R_* = 8.99 ± 0.42 for the 3.6 and 4.5 |$\mu$|m transits, respectively. We combine these measurements with other results to refine the system parameters and planet orbital ephemeris in Section 4.7.

To obtain measurements for the transmission spectrum, we first allowed the mid-transit time (T₀) and planet-to-star radius (R_p/R_*) to vary and held fixed the orbital period (P), eccentricity(e), inclination, and scaled semimajor axis (a/R_*) to the values reported in Table 2. We found central transit times T_mid = 2458779.05194 ± 0.00061 and 2458785.90231 ± 0.00045 (BJD_TDB) for the 3.6 and 4.5 |$\mu$|m channel, respectively. Finally, we repeated the fit allowing only the planet radius to vary and fixed the mid-transit times to their best-fitting values from the first fit. Figs 3, 4 and 7 summarize our results. Applying the wavelet method on the light-curve residuals, as detailed in Carter & Winn (2009), we found white and red noise components σ_w = 0.010 and σ_r = 0.0010 for the 3.6 |$\mu$|m transit data, and σ_w = 0.012 and σ_r = 0.0010 for the 4.5 |$\mu$|m transit, respectively.

3.2.2 Secondary eclipses

First, we fit the secondary eclipses by allowing the mid-eclipse time (T_ecl) and eclipse (or occultation) depth, δ_occ = (B_λ(T_p)/B_λ(T_*))(R_p/R_*)² to vary freely. We find T_ecl(3.6) = 2458787.6139 ± 0.0033 and T_ecl(4.5) = 2458794.4777 ± 0.0075. An analysis of the departures of these measurements from the predicted ephemeris is presented in Section 4.7. We measure white and red noise components of σ_w = 0.010 and σ_r = 0.003 for the 3.6 |$\mu$|m transit data, and σ_w = 0.013 and σ_r = 0.0013 for the 4.5 |$\mu$|m transit, respectively.

In a second pass, we fixed the central eclipse times to these values and allowed only the occultation depths to vary in the light-curve fits. We find δ_occ(3.6) = 1600 ± 250 and δ_occ(4.5) = 926 ± 334 parts-per-million (Fig. 8).

3.3 TESS photometry

We performed fits to the TESS light curves using the same method as for the Spitzer light curves. We assumed the following limb-darkening coefficients: c₁ = 0.6228, c₂ = −0.1870, c₃ = 0.4742, c₄ = −0.2224 and kept them fixed throughout the analysis. In addition, the planet orbital eccentricity and period were held fix to the values reported in Table 2. We first allowed all transit parameters to vary freely in the fit aiming to measure the system parameters. A weighted mean to the orbital inclination and normalized semimajor axis measured from the 2018 and 2020 light curves gave i = 85.11 ± 0.27° and a/R_* = 8.66 ± 0.26. These values are in excellent agreement with the VLT results, Table 2. Results for the central times are detailed in Table 3. In a second pass, we fixed the transit central times to the best-fitting values and allowed only the planet-to-star radius (R_p/R_*) to vary. Combining all measured TESS radii, we find a weighted mean value of R_p/R_* = 0.11986 ± 0.00061.

Table 3.

Central transit times.

Epoch	Mid-transit time, BJDTDB	O−C, d	Data	Note
−498	2456258.06290 ± 0.00020	−0.000027	WASP	1
0	|$2457963.84120^{+0.00034}_{-0.00033}$|	0.000534	VLT	2
7	2457987.81699 ± 0.00029	−0.000472	VLT	2
114	2458354.32081 ± 0.00094	0.000902	TESS	3
115	2458357.74406 ± 0.00086	−0.001104	TESS	3
116	2458361.17133 ± 0.00085	0.000909	TESS	3
117	2458364.59511 ± 0.00076	−0.000567	TESS	3
119	2458371.44639 ± 0.00081	0.000200	TESS	3
120	2458374.87130 ± 0.00082	−0.000147	TESS	3
121	2458378.29731 ± 0.00102	0.000607	TESS	3
148	|$2458470.77900^{+0.00370}_{-0.00150}$|	0.000371	HST	3
151	|$2458481.05503^{+0.00041}_{-0.00046}$|	0.000632	HST	3
238	2458779.05194 ± 0.00061	0.000226	Spitzer	3
240	2458785.90231 ± 0.00045	0.000083	Spitzer	3
329	2459090.7484xs2 ± 0.00084	−0.001636	TESS	3
330	2459094.17473 ± 0.00093	−0.000582	TESS	3
331	2459097.60105 ± 0.00085	0.000481	TESS	3
333	2459104.45013 ± 0.00098	−0.000952	TESS	3
334	2459107.87539 ± 0.00098	−0.000948	TESS	3
335	2459111.30260 ± 0.00098	0.001005	TESS	3

Note. 1: Hellier et al. (2014), 2: Nikolov et al. (2018b), 3: This work

3.4 Stellar elemental abundances, chromospheric activity, photometric variability

We obtained abundances of the main atmospheric constituents of the host star using the FEROS spectra, using methods similar to those given in Doyle et al. (2013). Results from our analysis are shown in Table 4. Abundances are given with respect to the Solar abundances from Asplund et al. (2009), except for lithium which is given as logarithm of the number ratio with respect to hydrogen plus 12 (e.g. Bergemann & Serenelli 2014, p.247).

Table 4.

Stellar elemental abundances of WASP-96.

Abundance	Value
[C/H]	+0.25 ± 0.19
[N/H]	+0.32 ± 0.16
[O/H]	+0.37 ± 0.11
[Na/H]	+0.31 ± 0.13
[K/H]	+0.27 ± 0.19
logA(Li)	1.69 ± 0.15
[Mg/H]	+0.47 ± 0.10
[Ca/H]	+0.18 ± 0.16
[Sc/H]	+0.14 ± 0.21
[Ti/H]	+0.18 ± 0.12
[V/H]	+0.22 ± 0.15
[Cr/H]	+0.25 ± 0.11
[Mn/H]	+0.49 ± 0.28
[Fe/H]	+0.21 ± 0.15
[Co/H]	+0.20 ± 0.15
[Ni/H]	+0.24 ± 0.15
[Y/H]	+0.17 ± 0.23
[Si/H]	+0.24 ± 0.16
[S/H]	+0.25 ± 0.11

We compared the tabulated and observed wavelengths of the Ca ii H&K line cores finding evidence for an offset of the K line wavelength as indicated in Fig. B1. This has been observed towards other exoplanet host stars observed with FEROS, including WASP-110 and can be attributed to the relatively low signal-to-noise ratio (Nikolov et al. 2021). We measure the activity index, S_HK, as estimated from the spectrum using the method described by Vaughan, Preston & Wilson (1978) and transformed to the |$\log R^{\prime }_{\rm HK}$| system using the relations from Noyes et al. (1984). We find the star to be chromospherically quiet with (⁠|$\log {\rm {R^{^{\prime }}_{H\&K}}}=-5.3\pm 0.1$|⁠). This places WASP-96 among some of the most chromospherically quiet stars hosting a transiting exoplanet.

To assess the level of photometric variability due to spots on the host star, we performed Lomb–Scargle periodogram analysis on the ASAS-SN light curve with trial frequencies ranging between 0.005 and |$1\, \rm {d^{-1}}$|⁠, corresponding to a period range between 1 and 200 d. We find a rotational modulation with a period of 17.8 d with a false alarm probability of ∼0.1. We phase-folded the light curve and fitted four-parameter sine wave to measure the amplitude of the flux variation. We found an amplitude of |$0.49\pm 0.11{{\ \rm per\ cent}}$| and residual scatter of 11 parts-per-thousand (ppt), Fig. B2. This amplitude exceeds the <1mmag (⁠|$95{{\ \rm per\ cent}}$|⁠), or |$\lt 0.1{{\ \rm per\ cent}}$|⁠, rotational modulation reported by Hellier et al. (2014) by a factor of ∼5. However, given the low value for chromospheric activity index in combination with the low photometric variability from the WASP photometry we conclude that the ASAS-SN result should be taken with a low weight. The TESS SAP PDC corrected light curves from sectors 2 and 29 have an average scatter of 3.4 ppt. They were examined using period04 (Lenz & Breger 2005). No variations above 0.3 ppt were found at periods longer than the orbital period of the system. Therefore, these data confirm the lack of variation found in the WASP photometer and rule out the variation suggested by the ASAS-NS data.

4 RESULTS AND INTERPRETATION

4.1 HST-Spitzer transmission spectrum

The measured HST and Spitzer wavelength-dependent relative planet radii, which comprise the planet’s transmission spectrum, are shown in Fig. 1. The HST spectrum exhibits the full absorption signature of the rotational–vibrational absorption band of water with peaks at ∼0.95, 1.15, and 1.4 |$\mu$|m covering ∼7 atmospheric pressure scale heights (with 1 scale height corresponding to ∼610 km, assuming an equilibrium temperature of T_eq = 1285 K, surface gravity of |$g=7.6\, \mathrm{m\,s^{-2}}$|⁠, and solar mean molecular weight of μ = 2.3 a.m.u.), in a wavelength region from ∼0.8 to 1.7 |$\mu$|m.

Figure 1.

Infrared transmission spectrum of the cloud-free WASP-96b. The error bars indicate the 1σ uncertainties. Top panel: Results from the two HST WFC3 spectroscopic observations. Lower panel: Combined HST WFC3 transit spectroscopy and TESS and Spitzer IRAC 3.6 and 4.5 |$\mu$|m transit photometry.

4.2 Combining VLT with HST and Spitzer transmission spectroscopy

Following the ultimate goal of this study i.e. extending the optical ground-based transmission spectrum of WASP-96b to infrared wavelengths, we complemented the VLT spectrum with the HST-Spitzer transit observations analysed here. Fig. 2 shows the combined ground-based and space-based spectra, presented in Nikolov et al. (2018b) and in this work, respectively. The infrared spectrum is found to be at a substantially larger radius of ΔR_p/R_* ∼ 0.004, or a change of transit depth of |$\sim 0.1{{\ \rm per\ cent}}$| (1000 ppm). Factors with the potential to contribute to such spectral offsets include: (1) stellar activity and photometric variability of the host star; (2) discrepancy in the system parameters assumed in the light-curve fits of the separate data sets; (3) accurate accounting of the limb-darkening; (4) unaccounted systematic errors in the light-curve analysis. The |$\lt 0.1{{\ \rm per\ cent}}$| photometric variability of the host star reported by Hellier et al. (2014), or even the higher variability of |$0.49\pm 0.11{{\ \rm per\ cent}}$| tentatively seen in the ASAS-SN monitoring, can only account for a negligible radius offset. In this study, we assumed identical system parameters for each of the data set and the limb darkening was properly treated in both studies. We note that the limb darkening of the host star is predicted to substantially vary from optical to infrared wavelengths. The VLT data have been fit with a quadratic limb-darkening law with the first coefficient allowed to vary to account for the variable transmission of the Earth atmosphere. In this study, we fixed the limb-darkening coefficients to their theoretical values from the same synthetic spectra assumed for the limb darkening of the ground-based data (Magic et al. 2015). A visual inspection of the light-curve fits shown in Extended Data figs 2 and 3 of Nikolov et al. (2018b) shows no evidence for significant unaccounted systematics. However, in that study (and here too), the authors performed a common-mode correction, which strongly depends on the assumption for system parameters and particularly on the planet-to-star radius measurement. Nikolov et al. (2018b) measured planet-to-star radii of R_p/R_* = 0.1147 ± 0.0014 and |$0.1168^{+0.0015}_{-0.0014}$| for the GRIS600B and GRIS600RI white light curves, respectively. While these values differ by only ΔR_p/R_* = 0.0021 ± 0.0021, the choice for either of them displaces the mean radius of the transmission spectrum by a substantial margin. The authors chose to rely on the blue (GRIS600B) transmission spectrum, which is less affected by telluric lines, and used the overlapping region of the two grisms to offset the red (GRIS600RI) transmission spectrum to the inherent for the blue (GRIS600RI) spectrum smaller planet radii. Choosing the red spectrum as a reference would account for |$50{{\ \rm per\ cent}}$| of the difference between the HST and VLT spectra. While Yip et al. (2021) concluded that combining data sets is risky and could not ascertain if the VLT and HST data were compatible, we point out that techniques such as common-mode corrections are well known to compromise the absolute transit depths, such that an overall transmission spectrum can have a large vertical uncertainty which is not fully captured by the point-to-point uncertainties of the relative transit depths (typically the quantity of interest). This is largely because the common-mode correction effectively assumes a broad-band R_p/R_* with no uncertainty, which ignores the actual uncertainty on the broad-band R_p/R_* from the white light-curve fit.

Figure 2.

Combined transmission spectrum of WASP-96b. Plotted are the VLT and HST spectra with orange and black symbols, respectively. A radius measurement from a combined analysis of TESS photometric light curves is indicated with the green symbol. Error-bars mark the 1σ uncertainties.

With a TESS measurement from a broad-band covering the wavelength region from 0.6 to 1 |$\mu$|m, one could employ this radius to offset the VLT spectrum. However, the TESS band is centred in the small overlapping region of the red (GRIS600RI) VLT and blue HST (G102) data, which hampers quantifying the size of the offset. We also constructed a band in the HST G102 data corresponding to the narrow overlapping region with the VLT data. However, we found that the uncertainty of such a narrow band of only 3 nm is insufficient to derive an accurate radius. Instead, we derive an offset correction for the VLT data by a simultaneous retrieval of the VLT and HST-Spitzer spectra as detailed in the next section.

4.3 Retrieval of the combined VLT, HST, and Spitzer transmission spectrum

We considered inverse methods to constrain atmospheric properties at the day–night terminator of WASP-96b and the radius difference between the VLT optical and HST-Spitzer infrared transmission spectra. We employed the 1D–2D radiative–convective equilibrium atmo model (Tremblin et al. 2015, 2016). Previous retrieval results using atmo can be found in Wakeford et al. (2017a), Evans et al. (2017), Mikal-Evans et al. (2019), Mikal-Evans et al. (2020), Spake et al. (2021), and Evans et al. (2018). atmo is capable of solving for the pressure–temperature (P-T) profile and chemical abundances that satisfy hydrostatic equilibrium and conservation of energy given a set of opacities. The code computes isotropic multigas Rayleigh scattering, H₂−H₂ and H₂−He collision-induced absorption, as well as opacities for all major chemical species taken from the most up-to-date high-temperature line list sources, including Na, K, Li, Rb, Cs, H₂, He, H₂O, CO₂, CO, CH₄ (Sauval & Tatum 1984; Barber et al. 2006; Heiter et al. 2008; Rothman et al. 2010; Tashkun & Perevalov 2011; Amundsen et al. 2014, 2017; Yurchenko & Tennyson 2014; Tremblin et al. 2015, 2016; Drummond et al. 2016; Goyal et al. 2017). Rainout of condensate species is also included (see discussion in Goyal et al. 2019). atmo can also fit a parametrized P-T profile and chemical abundances to transmission spectra. For the retrieval analysis in our study, we used the parametrized P-T profile detailed in Guillot (2010), which gives three free parameters: the Planck mean thermal infrared opacity (κ_IR); the ratio of optical to infrared opacities (γ_O/IR); and an irradiation efficiency factor (β). We set the internal temperature to 400 K based on Thorngren, Gao & Fortney (2019).

atmo treats hazes as parametrized enhanced Rayleigh scattering and clouds are assumed to add grey opacity. The code uses the correlated-k approximation with the random overlap method to compute the total gaseous mixture opacity, which has been shown to agree well with a full line-by-line treatment (Amundsen et al. 2017).

We first performed a retrieval assuming chemical equilibrium and allowed the radius, opacity from clouds and hazes, and elemental abundances of Na, K, C, and O allowed to freely vary. In the synthetic spectra, alkali line-wing shapes were computed assuming the predictions from Allard, Homeier & Freytag (2012). We allowed four elemental abundances to vary independently, as they are major species that are also likely to be sensitive to spectral features in the data, while the rest were varied by a trace metallicity parameter log (Z_trace/Z_⊙). This quantity is not the overall bulk metallicity but contains the abundances for trace species not otherwise individually fit (i.e. all except H, He, Na, K, C, and O). By varying the elemental abundances of Na, K, C, and O separately, we allow for non-solar compositions but with chemical equilibrium imposed such that each model fit has a chemically plausible mix of molecules given the retrieved temperatures, pressures, and underlying elemental abundances. Importantly, by varying both O and C separately (rather than a single C/O value) we alleviate an important modelling assumption that can affect the retrieved C/O value (see Drummond et al. 2019). To quantify the radius difference between the VLT and HST-Spitzer spectra, we included a parameter that controls the offset of the VLT spectrum (δ_VLT) and was allowed to freely vary in the retrieval. We chose to discard the TESS radius measurement owing to its large uncertainty and broad wavelength region.

The best-fitting transmission spectra are shown in Figs 3 and 4. The P-T profile is shown in Fig. 5, and the posteriors of all retrieved quantities are shown in Fig. 6. Our retrieval analysis finds a negligible contribution from cloud or haze opacity, which indicates that the atmosphere of WASP-96b is cloud and haze-free at the pressures being probed at the limb. This result is in agreement with the retrieval results from our VLT optical spectrum alone detailed in Nikolov et al. (2018b). We find an offset in the planet-to-star radius of |$\delta _{\mathrm{VLT}}=0.0043^{+0.00031}_{-0.00037}$| (or a transit depth change of 1004 ppm), which is in agreement with the results of Yip et al. (2021).

Figure 3.

Comparison of the VLT, HST, and Spitzer measured transmission spectra of WASP-96b (dots with 1σ vertical error bars; the bin widths are indicated with horizontal bars) with the best-fitting model from the retrieval analysis (red line) with resolution identical to the data, with the 1σ, 2σ, and 3σ confidence levels (dark to pale blue regions).

Figure 4.

Same as Fig. 3 but showing a zoom around the combined optical-infrared spectrum from VLT and HST. We binned with the weighted mean the VLT, G102, and G141 spectra in 2, 2, and 4 points, respectively.

Figure 5.

Pressure-temperature profile results from atmo assuming chemical equilibrium. Shown are the median profile of the posterior distribution (purple line) and the 1σ, 2σ, and 3σ model distributions with dark, medium, and light blue shaded regions, respectively. Cloud condensation curves, assuming solar abundance, are indicated with the dashed lines from Visscher, Lodders & Fegley (2006, 2010) and Wakeford et al. (2017b).

Figure 6.

Atmospheric properties for WASP-96b obtained from a retrieval with atmo, assuming chemical equilibrium. Marginalized posterior distributions of trace metallicity, planet radius, cloud and haze opacity, cloud top pressure, index of haze power, infrared opacity, optical-to-infrared offset, β and elemental abundances of carbon, oxygen, sodium, and potassium in solar units. Lack of clouds and hazes at the pressures being probed by our observations is evidenced by the high value of cloud top pressure and negligible haze opacity of the marginalized posterior distributions. The abundances of oxygen and sodium are the two constrained quantities. The mean abundances and 1σ uncertainties are indicated with the vertical continuous and dashed lines, respectively. The posterior distribution of the C/O ratio is reconstructed from the |$\mathrm{\log (C/C_{\odot })}$| and |$\mathrm{\log (O/O_{\odot })}$| abundances and is shown, with the solar value indicated with the blue dotted line in the top panel.

We obtain constraints of |$\mathrm{log(Na/Na_{\odot })}=1.32^{+0.36}_{-0.45}$| for the sodium abundance, |$\mathrm{log(K/K_{\odot })}=-0.05^{+0.49}_{-0.39}$| for the potassium abundance, |$\mathrm{log(C/C_{\odot })}=0.03^{+0.53}_{-0.45}$| for the carbon abundance, and |$\mathrm{log(O/O_{\odot })}=0.88^{+0.38}_{-0.43}$| for the oxygen abundance (Table 5 and Fig. 6). Sodium and oxygen are the sole two elements for which we place definitive constraints, which correspond to solar to supersolar values with abundances relative to the solar values of |$21^{+27}_{-14}$| (2.4σ) and |$7^{+11}_{-4}$| (2σ), respectively. Our results are in agreement with the measured supersolar abundances of the host star Table 4. The best-fitting model gives χ² = 120.96 for 103 degrees of freedom. We note that the posterior distribution of the C/O ratio shown in Fig. 6 has been reconstructed from the best-fitting retrieved log(C/C_⊙) and log(O/O_⊙) abundances. This implies that the prior ranges for the C/O ratio are constructed for the log(O/O_⊙) and log(C/C_⊙) parameters which are fit separately. One limitation of the data and retrieval modelling is that the carbon-bearing species are not resolved by the Spitzer data, including CO₂, CO, and CH₄. While our modelling assumes equilibrium chemistry to model the mixture of these gasses, it has previously been found that the C/O ratio can vary dramatically depending on whether free chemistry or equilibrium chemistry is assumed in the model (e.g. Spake et al. 2021).

Table 5.

Retrieval results.

Parameter, unit	Value	Prior range
|$\chi ^2_{\mathrm{min}}$|	120.96
Nfree	14
Ndata	117
log (Ztrace/Z⊙)	|$-0.40^{+0.48}_{-0.25}$|	−1 to 2
Rpl, (RJup)	|$1.2493^{+ 0.0039}_{-0.0036}$|	1.13 to 1.39
ln δcloud	|$-1.30^{+4.54}_{-3.36}$|	−10 to 10
ln δhaze	|$-2.85^{+2.79}_{-2.90}$|	−10 to 10
Cloud Top, log10(bar)	|$-1.84^{+0.83}_{-1.53}$|	−6 to 0.7
αhaze	|$2.45^{+1.13}_{-1.15}$|	0 to 5
δVLT × 103	|$4.29^{+0.31}_{-0.37}$|	0.29 to 8.3
log κIR	|$-2.73^{+0.96}_{-0.93}$|	−5 to 0.5
log γO/IR	|$-1.37^{+0.96}_{-1.15}$|	−4 to 1.5
β	|$0.78^{+0.13}_{-0.09}$|	0 to 2
|$\mathrm{\log (C/C_{\odot })}$|	|$0.03^{+0.53}_{-0.45}$|	−1 to 2
|$\mathrm{\log (O/O_{\odot })}$|	|$0.88^{+0.38}_{-0.43}$|	−1 to 2
|$\mathrm{\log (Na/Na_{\odot })}$|	|$1.32^{+0.36}_{-0.45}$|	−1 to 2
|$\mathrm{\log (K/K_{\odot })}$|	|$-0.05^{+0.49}_{-0.39}$|	−1 to 2

In addition to the retrieval with alkali metals profile from Allard et al. (2012), we performed a second retrieval assuming Lorentzian profile from Burrows, Marley & Sharp (2000). Our results indicate an excellent agreement between the retrieved atmospheric properties of the two assumed profiles without a preference of one versus the other.

We also performed retrievals excluding the Spitzer data. Because little-to-no carbon species are needed to fit the transmission spectrum (with low log(C/C_⊙) and resulting C/O ratio found), and the Spitzer data points have fairly large error bars, excluding the Spitzer data had little effect on the retrieved parameters. The log(C/C_⊙) parameter was found to increase by about 0.1 dex when excluding the Spitzer data, but this is well within the retrieval uncertainties (0.5 dex).

The combined VLT-HST-Spitzer transmission spectrum after a correction of the VLT spectrum is detailed in Table 6.

Table 6.

WASP-96b VLT, HST, and Spitzer transmission spectrum after a correction of the VLT spectrum.

λc, |$\mu$|m	Δλ, |$\mu$|m	Rp/Rs	|$\mathrm{\sigma _{R_p/R_s}}$|		λc, |$\mu$|m	Δλ, |$\mu$|m	Rp/Rs	|$\mathrm{\sigma _{R_p/R_s}}$|
0.37565	0.02565	0.11908	0.00152		1.12345	0.00465	0.11741	0.00120
0.40530	0.00400	0.11876	0.00182		1.13275	0.00465	0.11813	0.00124
0.41330	0.00400	0.11811	0.00167		1.14205	0.00465	0.11972	0.00111
0.42130	0.00400	0.11983	0.00167		1.15135	0.00465	0.12159	0.00116
0.42930	0.00400	0.11888	0.00205		1.16065	0.00465	0.12027	0.00114
0.43730	0.00400	0.11750	0.00163		1.16995	0.00465	0.11739	0.00123
0.44530	0.00400	0.11934	0.00129		1.17925	0.00465	0.11790	0.00114
0.45330	0.00400	0.11895	0.00123		1.18855	0.00465	0.11724	0.00127
0.46130	0.00400	0.11891	0.00125		1.19785	0.00465	0.11794	0.00100
0.46930	0.00400	0.11790	0.00116		1.20715	0.00465	0.11980	0.00105
0.47730	0.00400	0.11711	0.00100		1.21645	0.00465	0.11834	0.00106
0.48530	0.00400	0.11753	0.00130		1.22575	0.00465	0.11900	0.00129
0.49330	0.00400	0.11870	0.00104		1.23505	0.00465	0.11600	0.00122
0.50130	0.00400	0.11880	0.00138		1.24435	0.00465	0.11965	0.00088
0.50930	0.00400	0.11755	0.00112		1.25365	0.00465	0.11953	0.00105
0.51730	0.00400	0.11747	0.00101		1.26295	0.00465	0.11730	0.00092
0.52530	0.00400	0.11746	0.00096		1.27225	0.00465	0.11929	0.00095
0.53330	0.00400	0.11873	0.00089		1.28155	0.00465	0.12127	0.00126
0.54130	0.00400	0.11794	0.00095		1.29085	0.00465	0.11891	0.00117
0.54930	0.00400	0.11790	0.00089		1.30015	0.00465	0.11765	0.00128
0.55730	0.00400	0.11751	0.00102		1.30945	0.00465	0.11952	0.00129
0.56530	0.00400	0.11970	0.00095		1.31875	0.00465	0.11893	0.00124
0.57330	0.00400	0.11928	0.00095		1.32805	0.00465	0.12073	0.00109
0.58130	0.00400	0.12082	0.00084		1.33735	0.00465	0.11894	0.00094
0.58930	0.00400	0.12114	0.00092		1.34665	0.00465	0.11730	0.00136
0.59730	0.00400	0.12077	0.00086		1.35595	0.00465	0.11908	0.00139
0.60530	0.00400	0.12036	0.00086		1.36525	0.00465	0.12145	0.00117
0.61330	0.00400	0.11975	0.00081		1.37455	0.00465	0.12052	0.00106
0.62130	0.00400	0.11916	0.00096		1.38385	0.00465	0.12164	0.00136
0.62930	0.00400	0.12077	0.00099		1.39315	0.00465	0.12144	0.00118
0.63730	0.00400	0.11886	0.00096		1.40245	0.00465	0.12082	0.00097
0.64530	0.00400	0.11928	0.00091		1.41175	0.00465	0.11989	0.00109
0.65330	0.00400	0.11982	0.00100		1.42105	0.00465	0.12124	0.00125
0.66130	0.00400	0.11713	0.00089		1.43035	0.00465	0.11941	0.00106
0.66930	0.00400	0.11961	0.00081		1.43965	0.00465	0.12087	0.00135
0.67730	0.00400	0.11871	0.00117		1.44895	0.00465	0.12023	0.00097
0.68930	0.00800	0.11618	0.00078		1.45825	0.00465	0.11984	0.00122
0.70130	0.00400	0.11740	0.00090		1.46755	0.00465	0.11852	0.00099
0.70930	0.00400	0.11912	0.00096		1.47685	0.00465	0.12035	0.00099
0.71730	0.00400	0.11857	0.00085		1.48615	0.00465	0.11982	0.00113
0.72530	0.00400	0.11892	0.00097		1.49545	0.00465	0.12122	0.00094
0.73330	0.00400	0.12001	0.00094		1.50475	0.00465	0.11976	0.00136
0.74130	0.00400	0.11786	0.00081		1.51405	0.00465	0.11897	0.00112
0.74930	0.00400	0.11834	0.00096		1.52335	0.00465	0.11942	0.00095
0.76130	0.00800	0.11902	0.00075		1.53265	0.00465	0.11795	0.00128
0.77330	0.00400	0.11741	0.00109		1.54195	0.00465	0.11757	0.00113
0.78130	0.00400	0.11771	0.00096		1.55125	0.00465	0.11615	0.00105
0.78930	0.00400	0.11765	0.00123		1.56055	0.00465	0.11766	0.00127
0.79730	0.00400	0.11716	0.00117		1.56985	0.00465	0.11889	0.00112
0.82460	0.03430	0.11857	0.00099		1.57915	0.00465	0.11759	0.00124
0.87600	0.01710	0.11682	0.00079		1.58845	0.00465	0.11823	0.00115
0.91025	0.01715	0.11868	0.00074		1.59775	0.00465	0.11736	0.00091
0.94455	0.01715	0.11833	0.00076		1.60705	0.00465	0.11872	0.00139
0.97885	0.01715	0.11860	0.00069		1.61635	0.00465	0.11842	0.00129
1.01320	0.01720	0.11864	0.00072		1.62565	0.00465	0.11829	0.00140
1.04755	0.01715	0.11822	0.00080		1.63495	0.00465	0.11817	0.00130
1.08185	0.01715	0.11693	0.00104		3.60000	0.35460	0.12140	0.00130
1.11415	0.00465	0.11829	0.00101		4.50000	0.45020	0.11890	0.00120
1.11615	0.01715	0.11984	0.00063

4.4 Forward model analysis of the combined VLT, HST, and Spitzer transmission spectrum

To interpret the combined WASP-96b observed transmission spectrum, we also used a forward model grid of transmission spectra generated using self-consistent P-T profiles in atmo for a range of metallicity, C/O ratio, and energy redistribution, described in detail in Goyal et al. (2020). Here, self-consistent imply radiative–convective equilibrium P-T profiles consistent with equilibrium chemical abundances. Prior to fitting the data, we offset the VLT spectrum with by δ = 4.29, found in the retrieval analysis. The best-fitting forward model transmission spectrum is shown in Fig. 7. This fit corresponds to χ² = 129.3 and reduced χ² = 1.11 with 116 degrees of freedom (117 observed points minus 1 for the vertical offset in the transmission spectrum). This fit could be considered as a relatively good fit given that the model transmission spectrum is purely forward with self-consistent P-T profiles. The best-fitting transmission spectrum has a recirculation factor of 0.25 (high-energy redistribution), indicating a relatively cooler P-T profile as expected for transmission spectrum, since we probe the day–night terminator region of the exoplanet using the transmission spectrum. We also determine that the best-fitting forward model transmission spectrum has solar metallicity and sub-solar C/O ratio (0.35). Only one other model spectrum from the grid lie within the 1σχ² value of the best-fitting transmission spectrum and this model also has solar metallicity and sub-solar C/O ratio (0.35).

Figure 7.

Comparison of the VLT, HST, and Spitzer measured transmission spectrum (black dots with 1σ error bars) with the best-fitting forward model transmission spectrum (blue) from the self-consistent planet-specific grid of WASP-96b across a range of recirculation factor (RCF), metallicity (Z), and C/O ratio. χ² = 129.3 for the best-fitting model spectrum. The best-fitting model parameters are RCF = 0.25, solar metallicity, and C/O ratio of 0.35.

4.5 Thermal emission spectrum

The two measured Spitzer IRAC eclipse depths are plotted in Fig. 8 and comprise the planet’s dayside thermal spectrum. We find the 4.5 |$\mu$|m eclipse depth to be lower compared to the 3.6 |$\mu$|m data, which suggests a decreasing temperature with altitude.

Figure 8.

Emission spectrum of WASP-96b, synthetic emission spectra, and assumed pressure-temperature profiles. Left-hand panel: WASP-96b thermal spectrum from Spitzer IRAC 3.6 and 4.5 |$\mu$|m observations (black symbols with 1σ error bars) compared to forward synthetic spectra assuming an isothermal (dashed line), atmo and petitRADTRANS decreasing (continuous and dotted lines, respectively), and inverted (dash–dotted line) pressure-temperature profiles. Model predictions, averaged over the wavelengths for each channel, are indicated with the open symbols and colour coded to the relevant model. The lowest chi-square spectrum corresponds to the atmo synthetic spectrum, assuming a non-inverted pressure-temperature profile (continuous line). Given the relatively low level of irradiation for WASP-96b, an inverted temperature profile is rather unlikely. Yet, we include a spectrum for the inverted P-T profile merely for reference purposes, not expecting it to fit the data. Right-hand panel: Pressure-temperature profiles assumed in the calculation of the thermal emission spectra.

We first fit a blackbody synthetic spectrum to determine the dayside temperature of the planet. The blackbody temperature was the sole parameter allowed to freely vary in this fit. We assumed a blackbody spectrum for the host star with the effective temperature of WASP-96. We found that the best-fitting spectrum for the planet, with χ² = 6.7 for 1 degree of freedom, and BIC = 7.4, corresponds to a temperature of T_p = 1545 ± 90 K. This spectrum corresponds to the dashed line in Fig. 8.

We then employed the publicly available radiative transfer code petitRADTRANS (Mollière et al. 2019) to compute synthetic spectra for the planet’s dayside. petitRADTRANS computes planet emission spectra for an assumed pressure-temperature (P-T) profile and gravity. Because exoplanet emission spectra are known to depend on the underlying P-T profile, we chose to compute spectra, assuming P-T profiles of a decreasing (non-inverted), and an increasing (inverted) temperature with altitude. WASP-96b is irradiated rather modestly, and we would not expect to encounter an inverted temperature profile. We note that the spectrum for the inverted case is included merely for reference purposes, not expecting it to fit the data. We used the analytic P-T profile of Guillot (2010). We assumed solar abundances for the main atmospheric constituents, including H₂O, CO, CO₂, and CO₄ and a solar mean molecular weight of 2.33. Both spectra were computed assuming unity for the recirculation factor i.e. no redistribution of the input stellar energy in the planet’s atmosphere. Prior to comparing synthetic with observed data, we averaged the synthetic spectra within the wavelength bins of the observed spectrum and computed the corresponding χ². We find the observed 3.6 and 4.5 |$\mu$|m eclipse depths to be best explained by a spectrum assuming the non-inverted profile with a chi-square of χ² = 4.6. The spectrum assuming an inverted profile resulted in χ² = 8.8. We note that complementary observations at a higher signal-to-noise ratio and spectral resolution would be needed to confidently confirm this result.

We also interpreted the Spitzer IRAC eclipse depths with a self-consistent planet-specific forward model grid of emission spectra as shown in Goyal et al. (2021). These emission spectra were generated using radiative–convective equilibrium P-T profiles consistent with equilibrium chemical abundances in atmo for WASP-96b, for a range of metallicity, C/O ratio, and energy redistribution, described in detail in Goyal et al. (2020). The best-fitting model emission spectrum and the corresponding P-T profile is shown in Fig. 8. We find the observed 3.6 and |$4.5\, \mu$|m eclipse depths to be best explained by a spectrum with non-inverted P-T profile with χ² = 3.8. The best-fitting model emission spectrum corresponds to a recirculation factor of 1, implying none of the energy is redistributed to the night side and therefore the P-T profile corresponding to this best-fitting model emission spectrum is one of the hottest in the grid of models for WASP-96b. The best-fitting model emission spectrum corresponds to a metallicity of 200 times solar and C/O ratio of 0.7. However, by looking at all the model spectra that lie within the 1σχ² value of the best-fitting emission spectrum we determine that a large range of metallicities (From 1x to 200x solar) and C/O ratio (0.35 to 0.7) can explain observed eclipse depths due to large uncertainty in the |$4.5\, \mu$|m eclipse depth. Interestingly, all these model spectra similar to the best-fitting model spectrum have non-inverted P-T profiles and correspond to a recirculation factor of 1 (no-redistribution), driven by the higher |$3.6\, \mu$|m eclipse depths requiring higher temperatures. The strong absorption in the |$4.5\, \mu$|m band is a result of CO, leading to lower eclipse depth in this band as seen in left side of Fig. 8 for the best-fitting model. Therefore, the best-fitting model has a non-inverted P-T profile which results in this absorption feature (instead of an inverted P-T profile that would lead to an emission feature due to CO, in the |$4.5\, \mu$|m band). All the model spectra that lie within the 1σχ² value of the best-fitting emission spectrum show this CO absorption feature but with varying degree of strength (i.e model |$4.5\, \mu$|m eclipse depth).

4.6 Refined system parameters

To improve the planet orbital inclination and semimajor axis of WASP-96b, we combined our measurements with literature results. In particular, our Spitzer and TESS measurements were complemented with the VLT results reported in Nikolov et al. (2018b). Weighted mean results for the orbital inclination i and normalized semimajor axis a/R are reported in Table 7. Both results agree within 2σ with the measurements of Hellier et al. (2014).

Table 7.

System parameters and orbital ephemeris.

Parameter, symbol, unit	Value
Inclination, i, (°)	85.135 ± 0.081
a/R*	8.820 ± 0.090
Orbital period, P, (d)	3.42525650 ± 0.00000037
Transit central time, T0, (BJD)	2457963.84067 ± 0.00012
ΔtO-C(3.6), (min)	−1.4 ± 4.8
ΔtO-C(4.5), (min)	18 ± 11
ecos (ω)	−0.00790 ± 0.00010

4.7 Improved planet orbital ephemeris

We combined the central transit times from our HST, Spitzer, and TESS light-curve analysis with the transit times reported by Hellier et al. (2014) and Nikolov et al. (2018b) to derive an updated transit ephemeris of WASP-96b. These two studies, respectively, include a single transit mid-time from a combined analysis of WASP, EulerCam, and TRAPPIST light curves, and central times from two VLT FORS2 white light curves. All central times were converted to BJDs (using the utilities detailed in Eastman, Siverd & Gaudi 2010) and fitted with a linear function of the planet orbital period (P), zero-epoch central transit time (T₀), and transit epoch (E): T_C(E) = T₀ + EP. Our results for the fitted parameters are reported in Table 7 and resulted in a reduced chi-square of |$\chi ^2_r=1.03$| for 18 degrees of freedom. The updated ephemeris is in an excellent agreement with the ephemeris reported by Hellier et al. (2014) and Yip et al. (2021) and does not indicate the presence of transit timing variations. Results for the measured mid-transit times, observed minus computed transit times (O−C), are tabulated in Table 3 and plotted in Fig. A7.

We used the refined orbital ephemeris to predict the central times of the observed secondary eclipses. Significant departures of mid-eclipse times, measured as the difference (Δt_O-C) between the observed and computed times, would indicate a non-circular planet orbit. Table 7 summarizes the values for Δt_O-C(3.6) and Δt_O-C(4.5), which agree with the predictions at 1σ and ∼1.6σ and indicate no departures from a circular orbit. Finally, we used equation (6) from Wallenquist (1950) to calculate ecos ω for the two eclipses. Table 7 details the weighted mean value from the two observations.

5 DISCUSSION

The sodium abundance obtained from our retrieval analysis of the combined optical-infrared transmission spectrum, |${\rm {\log (Na/Na_{\odot })}}$|  = |$1.32^{+0.36}_{-0.45}$|⁠, is about one order of magnitude larger, but is still compatible with the previous determination using only the VLT spectrum, |$\log {\epsilon _{\mathrm{\, Na}}}=6.9^{+0.6}_{-0.4}$| i.e. |${\rm {\log (Na/Na_{\odot })}}$|  = |$0.66^{+0.6}_{-0.4}$| at the ∼1.2σ level.

5.1 A comparison with an earlier result

An analysis of the HST G102 and G141 transit spectroscopy and the Spitzer IRAC 3.6 and 4.5 |$\mu$|m photometric transits has been presented by Yip et al. (2021). We find the two spectra to be slightly offset by a mean difference between the spectrum from this work and the literature spectrum of ΔR_p/R_s = (47.36 ± 0.24) × 10⁻⁵, measured using the VLT wavelengths. This offset is consistent the uncertainty of R_p/R_s found in our analysis of the white light curve (see Table 2). After applying the mean difference, the two spectra are in a respectable agreement with regards to the spectral shape (Fig. 9). The Spitzer IRAC radii are found to be in a good agreement at the ∼1σ confidence level. We note significantly smaller uncertainties, a factor of ∼2.2 of the Spitzer radii measured in our work compared to the results of Yip et al. (2021). Factors that can contribute to such significant difference include unaccounted systematics effects or fits performed to binned light curves. Yet, we are unsure of the exact reasons. We additionally note a marginal difference in the radii of the VLT spectrum and their uncertainties reported by Yip et al. (2021). Finally, the G102 bin widths of the literature result reported in their table 6 are a factor of 10 smaller than the plotted bin widths throughout that study. We are unsure for the reason that led to this misquotation of the published results. Despite the differences between the literature transmission spectrum and the spectrum reported in this work, both studies find abundance results in excellent agreement using different analysis and retrieval methods.

Figure 9.

Comparison between the result in this work with those of Yip et al. (2021). An offset of ΔR_p/R_s = (47.36 ± 0.24) × 10⁻⁵ has been added to the literature result measured using the VLT spectrum. The horizontal bars represent 1σ uncertainties and the horizontal bars are the wavelength bins used to measure the transmission spectra.

5.2 Exoplanet abundances in the Solar system context

Elemental abundances of planetary atmospheres are a consequence of the joint influence of atmospheric chemical processes, formation conditions, and migration (Madhusudhan 2012; Fortney et al. 2013; Welbanks et al. 2019; Line et al. 2021; Lothringer et al. 2021). In the standard core-accretion paradigm, as they form, giant planets accrete H/He-dominated gas, as well as planetesimals that enrich their H/He envelopes in metals. Higher metal enrichment is expected for low-mass H/He envelopes, as these envelopes have smaller amount of gas for the in-falling metals to be mixed into. This is the case of the gas giants Jupiter and Saturn and the ice giants Uranus and Neptune in our own Solar system, which have elemental abundances that increase with a decreasing mass and also coincides with an increasing distance from the Sun (Fig. 10, Anders & Grevesse 1989; Asplund et al. 2009; Atreya et al. 2016, 2020).

Figure 10.

Elemental abundance ratios in the atmospheres of our own Solar system planets Jupiter, Saturn, Uranus, and Neptune, meteorites and the Halley’s comet, relative to the protosolar values (adapted from Atreya et al. 2020). The oxygen abundance of Jupiter has been updated to the latest measurement reported by  Li et al. (2020). Elemental abundance ratios for the atmosphere of WASP-96b relative to the protstellar abundances have been obtained with the same assumption as for the protosun i.e. by increasing the protostellar value by +0.04  (Asplund et al. 2009). Oxygen and sodium are the only two abundances obtained for WASP-96b (red symbols) and the grey triangles indicate unbounded abundances for potassium and carbon.

Our abundance and retrieval analysis show that both the host star and planet in the WASP-96 system have solar-to-supersolar elemental abundances of their atmospheres (Tables 4 and 5). While metal enrichment of planets does not appear to correlate with the metallicity of the parent star (Teske et al. 2019), the case of WASP-96 is a prime example for a significant enrichment of both, which likely reflects the initial formation conditions of this system. Comparing the oxygen abundance ratio for the atmosphere of WASP-96b relative to its host star, we find an oxygen ratio of 3.2 ± 0.4, which is similar to the oxygen ratio of Jupiter’s atmosphere and the Halley’s comet in our Solar system relative to the oxygen abundance of the Sun (Anders & Grevesse 1989). We find a higher ratio of |$10^{+0.4}_{-0.5}$| for the abundance of sodium in the atmosphere of WASP-96b relative to its host star, which is in agreement with the corresponding value for the Halley’s comet relative to the Sun, which does not hold for Jupiter. A recent study by Bhattacharya et al. (2021) shows evidence for alkali elements in the deep atmosphere of Jupiter from the Juno microwave data pointing towards severely subsolar values. This implies a much lower abundance ratio of sodium for the atmosphere of Jupiter compared to the solar value.

6 SUMMARY AND CONCLUSIONS

We presented HST WFC3 transit and Spitzer IRAC transit and eclipse observations of the atmosphere of the hot Saturn WASP-96b. This atmosphere has been identified by our team as cloud-free from a reconnaissance optical transmission spectroscopy made with the VLT. The ground-based spectrum exhibits the full pressure-broadened sodium profile, absorption from potassium and lithium, and scattering from H₂ in the near-UV. Our team has motivated the HST WFC3 observations with the goal to extend the planet’s transmission spectrum at infrared wavelengths to search for the theoretically predicted absorption from H₂O, CO, and CO₂, and to constrain atmospheric composition. The secondary eclipse observations were aimed at constraining the planet’s dayside temperature. The HST-Spitzer observations reveal the full absorption signature of water, which is an independent new evidence that the planet’s atmosphere is free of clouds at the terminator and at pressures being probed. Similar to earlier analysis of the HST-Spitzer transits, we find the HST-Spitzer spectrum to be at a significantly larger planet-to-star radius relative to the VLT level. We confirm and correct for a spectral offset of |$\Delta R_{{\rm p}}/R_{\ast }=(-4.29^{+0.31}_{-0.37})\, \times 10^{-3}$| of the VLT data relative to the HST-Spitzer spectrum. This offset can be explained by the assumed radius and corresponding uncertainty for the common-mode correction of the VLT spectra, which is a well-known feature of ground-based transmission spectroscopy. We find evidence for a lack of chromospheric and photometric activity of the host star and therefore have an insufficient contribution to the offset. We measure abundances for Na and O that are consistent with solar-to-supersolar, with abundances relative to solar values of |$21^{+27}_{-14}$| (2.4σ) and |$7^{+11}_{-4}$| (2σ), respectively. While combining transmission spectra from space and the ground is not trivial, finding the first evidence for a cloudless atmosphere from the ground and the following independent confirmation from space-based observations at new, highly complementary wavelengths, our results demonstrate the key role of ground-based low-resolution transmission spectroscopy, particularly the VLT with FORS2, in exploring the atmospheres of transiting exoplanets. A fit of a spectrum of a blackbody isothermal atmosphere to the thermal spectrum constrains a dayside brightness temperature of T_p = 1545 ± 90 K. While metal enrichment of exoplanets has been shown not to correlate with the metallicity of their parent stars, the case of WASP-96 is a prime example for a significant enrichment of both the star and planet, which likely reflects the initial formation conditions of this system.

ACKNOWLEDGEMENTS

Some/all of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This work is based on observations with the NASA/ESA Hubble Space Telescope (GO-15469, Nikolov et al. 2018c), obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. This work is based on observations made with the Spitzer Space Telescope (Nikolov et al. 2019), which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Based on observations collected at the European Southern Observatory under ESO programme 199.C-0467(H). This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA’s Science Mission directorate. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The authors are grateful to the anonymous referee for their thoughtful comments and suggestions, which helped to improve the manuscript. NN acknowledges support for this work by NASA through grants under the HST-GO-15469 program from the STScI. The research leading to these results received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 336792. JMG and NJM acknowledge support from a Leverhulme Trust Research Project Grant. This work used the Python Gaussian process library George. NM acknowledges funding from the UKRI Future Leaders Scheme (MR/T040866/1), Science and Technology Facilities Council Consolidated Grant (ST/R000395/1), and Leverhulme Trust research project grant (RPG-2020-82). ALC is supported by a grant from STScI (JWST-ERS-01386) under NASA contract NAS5-03127.

DATA AVAILABILITY

Raw and calibrated Hubble Space Telescope spectral transit time series and Spitzer Space Telescope transit and eclipse time series photometry are publicly available at the Mikulski Archive for Space Telescopes (MAST; https://archive.stsci.edu) and the NASA/IPAC Infrared Science Archive (IRSA; https://sha.ipac.caltech.edu/applications/Spitzer/SHA/), respectively. TESS light curves are publicly available at the MAST archive. Calibrated and extracted high-resolution FEROS spectra are publicly available via the European Southern Observatory’s Spectral Data Products Query Form (http://archive.eso.org/wdb/wdb/adp/phase3_spectral/form). Broad-band light curves are publicly available at the webpage of the All-Sky Automated Survey for Supernovae (ASAS-SN; https://asas-sn.osu.edu).

Footnotes

1

We used the publicly available george Gaussian Process Python suite (Foreman-Mackey 2015).

APPENDIX A: LIGHT CURVES AND AUXILIARY PARAMETERS

Plots of the separate time series observations are shown in Figs A1–A6, complementary detrending data is presented in Fig. A1 and observed minus computes (O-C) diagram is exhibited Figure A7.

Figure A1.

HST WFC3 white light curves and spectroscopic properties. The left- and right-hand panels refer to the G102 (blue) and G141 (red) data, respectively. The filled and open symbols, respectively, indicate retained and discarded data in our light-curve analysis. a-b, Normalized raw flux. c-d, Measured background for each frame in units of electrons per pixel, using the last individual read of each full scan. e-f, Median-subtracted drift of the spectra along the dispersion axis in units of pixels, estimated by cross-correlating the spectra with respect to the first spectrum. g-h, Same as the previous panel, but for the cross-dispersion axis measured with flux-weighted mean.

Figure A2.

HST  WFC3 white light curves. The left- and right-hand panels show the G102 (blue) and G141 (red) observations, respectively. The top row shows normalized raw light curves along with the marginalized Gaussian process model, assuming Matern 3/2 kernel (grey). The latter has been used as a common mode (wavelength invariant) correction for the spectroscopic light curves. The second row shows the detrended light curves along with the best-fitting transit model (red). The third row shows the best-fitting light-curve residuals and 1σ error bars, obtained by taking the difference of the raw flux and marginalized transit and systematics models. The dashed lines indicate the light curve dispersion, ∼75 parts per million.

Figure A3.

Spitzer  IRAC 3.6 and 4.5 |$\mu$|m photometry (left and right, respectively).Top panels: raw flux and the best-fitting transit and systematics model; middle panels: detrended light curves and the best-fitting transit models and binned by 8 min; lower panels: light-curve residuals.

Figure A4.

Spitzer  IRAC 3.6 and 4.5 |$\mu$|m eclipse photometry (left and right, respectively). Top panels: raw flux and the best-fitting transit and systematics model; middle panels: detrended light curves and the best-fitting eclipse models and binned by 10 min; lower panels: light-curve residuals.

Figure A5.

Spectrophotometric light curves from grism G102 offset by a constant for clarity. The first panel displays raw flux. We used dots and square symbols for the odd and even light curves, respectively, to better distinguish the individual light curves. The second panel shows the common-mode (CM)-corrected light curves with transit and systematics models (lines) with the highest statistical weight. The third panel is similar to the second one, but showing the detrended light curves and transit models. The fourth panel shows residuals with 1σ uncertainties. The dashed lines indicate the median residual values, with the dotted lines indicating the residual dispersion level and the percentage of the theoretical photon noise limit reached (blue) for each channel.

Figure A6.

Same as Fig. A5 but for grism G141. Shown is a subset of 20 light curves.

Figure A7.

Observed minus computed (O-C) transit times from the best-fitting orbital period and central transit times reported in the literature (orange: combined analysis of WASP, EulerCam and TRAPPIST light curves; cyan: VLT transits) and derived from our analysis of TESS (magenta), HST (green), and Spitzer (yellow) light curves. The error bars and dotted lines indicate the 1σ uncertainties of the transit times and derived ephemeris, respectively.

APPENDIX B: STELLAR ACTIVITY AND PHOTOMETRIC VARIABILITY

This section contains figures from our analysis of the FEROS high-resolution spectroscopy, shown in Fig. B1 and an ASAS-SN light curve is presented in Fig. B2.

Figure B1.

FEROS spectroscopy around the Ca ii H&K lines for WASP-96. The green line presents the original data and the black line indicates a Gaussian smoothed version for clarity. The tabulated laboratory wavelengths of the cores of the H|$\&$|K lines are plotted with the purple vertical lines.

Figure B2.

Phase-folded photometry from the ASAS-SN observing campaigns with the period of highest significance, as determined in our periodogram analysis. Shown are the unbinned photometry (red), binned light curve (blue), and the best-fitting sine function (green). The vertical bars indicate 1σ uncertainties.