The Hubble PanCET program: Transit and Eclipse Spectroscopy of the Strongly Irradiated Giant Exoplanet WASP-76b

WASP-76A has a companion star. WASP-76B was first discovered by (Wollert & Brandner, 2018) through lucky imaging with a separation of 0.425” $\pm$ 0.012” and position angle of 216.9^∘ $\pm$ 2.93^∘. Due to the small separation, light from WASP-76B is well mixed with WASP-76A in our HST spatial scan spectrum which causes a dilution effect on the extracted planet spectrum (Crossfield et al., 2012). To correct for this dilution effect, the companion stellar spectral type needs to be determined, and the extra flux contribution removed. The temperature of WASP-76A is 6250 $\pm$ 100K (West et al., 2016) and the updated distance from GAIA (Gaia Collaboration, 2018b) is 195.31 $\pm$ 6.03 parsecs. There are a total of three spatially resolved images of the WASP-76 system in the archive, taken with different filters (Table 1) using the Space Telescope Imaging Spectrograph (P.I.s: David Sing & Mercedes López-Morales), and Keck-AO with NIRC2 (P.I.: Brad Hansen), all shown in Figure 1. To determinate the spectral type of WASP-76B, we performed a two-component SED fit (Rodriguez et al., 2019) for WASP-76 and we determine the radius and temperature of WASP-76B to be $R_{\star}=0.795\pm{0.055}R_{\sun}$ and $T_{eff}=4850\pm{150K}$ which we then used as the prior in a EXOFASTv2 (Eastman, 2019) global analysis (Table 7). Within the EXOFASTv2 fit, the host star parameters are constrained using the we use the MESA Isochrones and Stellar Tracks (MIST) stellar evolution models (Paxton et al., 2015; Dotter et al., 2016; Choi et al., 2016). For the EXOFASTv2 fit, we included six new light curves (Fig. 16) from EulerCAM (Ehrenreich et al., 2020), Hazelwood and MVRC observations in addition to the transit and RV data used in the discovery paper (West et al., 2016) to refine and update the system ephemeris.

With the best-fit radius and effective temperature for both WASP76A and WASP76B, we can then use PHEONIX stellar models to calculate the flux contribution from both stars and the dilution effect of transit and eclipse depth can be corrected as follows:

where $F_{B}$ and $F_{A}$ are the flux contribution from the companion and the primary star at a given wavelength range. Since the companion star is spatially resolved at different levels with being mostly resolved in STIS spectra while completely blended in at Spitzer bands, we purposefully choose larger aperture sizes at all wavelength when extracting stellar spectra to ensure all companion flux contributions are included. Finally, the dilution is applied across the entire transit and eclipse spectra for consistent correction.

To propagate the uncertainties on the effective temperature of both stars into dilution factors and the final planet spectrum, we adopted the bootstrapping method used in (Stevenson et al., 2014) by generating 10000 PHOENIX stellar models for each star with $T_{eff}$ randomly sampled from a Gaussian distribution based on the $T_{eff}$ uncertainty. By calculating the corresponding dilution factors for each PHOENIX model pair, we obtain a 10000 sample size distribution of dilution factors at each wavelength bin. The final dilution factors are the median values of each distribution and the uncertainties will be the corresponding one sigma values which are then propagated into the reduced planet transit and eclipse spectra.

We observed WASP-76b in transit with 2 visits using HST/STIS G430L and 1 visit using the G750L grating (Table 1). Both gratings were observed using the ACCUM mode with the 50X2 aperture to minimize any slit losses. CCD subarray of 128 x 1024 pixels was used to reduce readout time and maximize observing efficiency. Each frame has an exposure time of $\sim$148 seconds and each orbit has $\sim$16 exposures. The combination of these two gratings provided a complete wavelength coverage from 2900 to 10300 Å. One prominent source of systematics in STIS light curves comes from the orbital motion of the telescope during the observations (Nikolov et al., 2014). As the telescope orbits between the day side and night side of the Earth, it experiences thermal expansion and contraction. This effect manifests as a varying observed flux as a function of telescope orbital phase.

Our data analysis process follows the standard methodology detailed in Sing et al. (2011) and Nikolov et al. (2015). We fit the STIS transit light curves using a combination of transit and instrument systematics models. The transit model is based on the analytic formula developed by Mandel & Agol (2002), and the systematics model is a fourth-order polynomial of the telescope orbital phase, a linear time term and wavelength shift ($\omega$) for each frame. Orbital inclination and $a/R_{star}$ are both fixed at the best-fit values derived in this paper during the fit. For limb darkening we calculate the relevant coefficients with ATLAS stellar models in the same way as detailed in Nikolov et al. (2015). The raw, corrected light curves and corresponding residuals for all three visits are shown in Figures 6, 7 and 8.

We observed both transits and eclipses using HST/WFC3 G141 in spatial scan mode to maximize photon-collecting efficiency (Deming et al., 2013). All frames used SPARS10 and NSAMP=16, with an exposure time of $\sim$104 seconds, and a forward and backward scan to maximize observing efficiency. Due to occultation of the telescope by the Earth, a $\sim$45 min gap exists between every HST orbit. In total, there are 5 orbits per visit and $\sim$19 spectra per orbit. Two orbits are pre-transit, two are in-transit, and one is post-transit.

The automatic CalWF3 pipeline does not include spatial scan mode, therefore additional processing is required before extracting the 1D spectra. We followed the standard procedures of background subtraction and energetic particle removal by flagging outliers relative to the median value along the vertical scan direction (Wakeford et al., 2013). Next, we corrected for the wavelength shift of each spectrum in the horizontal direction. To calculate the sub-pixel level shifts between each frame, we first summed each frame in the vertical direction to obtain a 1D spectrum and normalized it by its own median flux. Then we used $scipy.interpolate.interp1d$ function to interpolate normalized flux of each 1D spectrum in the wavelength direction relative to its pixel positions. Next we applied sub-pixel shifts to each 1D spectrum relative to a reference spectrum and calculated the shifts by minimizing the normalized flux differences between them. Finally we applied the calculated shifts on every 1D spectrum to obtain the wavelength shifts corrected 1D spectra. The hydrogen Paschen-beta line at 1.28 $\mu$m in the star is used to establish the zero-point of the wavelength calibration.

HST/WFC3 time series spectra often exhibit a ramp-like systematic shape when observing bright stars in high-cadence (Wilkins et al., 2014). This effect is attributed to charge trapping in the WFC3 HgCdTe infrared detector (Kreidberg et al., 2014b; Zhou et al., 2017). As initial photons arrive at the beginning of each orbit, some charge carriers can be trapped by impurities in the detector and cause lower readout signals. When all available traps are filled during the orbit, the measured signals asymptotically approach a constant level (Fig 4). The double-ramp shape per orbit is due to differences in exposure timing and telescope pointing between forward and backward scan. The timings for when each pixel receives light are different in forward and backward scan, and that can affect the ramp shape. Moreover, the illumination pattern on the pixel grid is slightly different from forward to backward scan. Since each pixel has a different number of charge traps, a constant offset in measured flux can occur when different portions of the detector are illuminated by forward and backward scan.

We observed transits and eclipses of WASP-76b with Spitzer at 3.6 and 4.5 $\mu$m (Table 1). Unlike HST, Spitzer is able to continuously observe targets for the entire transit and eclipse duration. One eclipse at each of 3.6 and 4.5 $\mu$m are reported by Garhart et al. (2020), and we do not re-analyze those data here. We here analyze the transit data, and two additional eclipses at 3.6 $\mu$m from program 13038. Our analysis of two additional eclipses at 3.6 $\mu$m followed the exact same procedures used by Garhart et al. (2020), in fact the same codes (implemented by D.D.). The new eclipse depths are included in Table 3.

Spitzer’s primary systematic effect comes from intra-pixel sensitivity variations coupled to pointing jitter, overlaid by temporal ramps. We correct for this combination of systematic effects using the pixel-level decorrelation (PLD) technique developed by Deming et al. (2015), with the implementation of the fit being the same as described by Garhart et al. (2020). PLD takes advantage of the total flux conservation within the aperture containing the star, and utilizes the relative flux contribution of individual pixels as basis vectors in the fit. This technique eliminates the need for finding the centroid position of the star while being capable of effectively removing red noise and flat-fielding inaccuracies.

Our solutions for the Spitzer transit depths incorporate quadratic limb darkening coefficients calculated for the Spitzer bands by Claret et al. (2013). These produce excellent agreement with the observed transit curves. Given that limb darkening is a minimal effect at Spitzer’s wavelengths, we adopt the Claret coefficients without further perturbation. Our initial procedure was to also freeze the orbital parameters at previously-determined values, since our experience with other data shows that this simple method usually produces excellent agreement with the shape of Spitzer’s observed transit curves. However, atmospheric characterization can be sensitive to alternate treatments of the orbital parameters (Alexoudi et al., 2018). Given also that we find some differences between the transit depths observed at 3.6- versus 4.5 $\mu$m, and between two transits at 3.6 $\mu$m, we explored other treatments of the orbital parameters. We used independent Gaussian priors for the two parameters that most affect the transit shape (orbital inclination and $a/R_{s}$), based on the discovery results from West et al. (2016). Those fits produced transit depths that differed minimally from fits that froze the orbital parameters at the West et al. (2016) values. Those differences (orbital priors minus orbital freeze), were 157 and 34 ppm in $R_{p}^{2}/R_{s}^{2}$ for the two transits at 3.6 $\mu$m, and -93 ppm at 4.5 $\mu$m. Our best-fit values of inclination and $a/R_{s}$ differed from West et al. (2016) by less than $1\sigma$. We also explored freezing the orbital parameters at the values derived in this paper, noting that our values for inclination and $a/R_{s}$ are within $1\sigma$ of West et al. (2016). Those transit depths differed from our initial values by 84 and 3 ppm at 3.6 $\mu$m, and 128 ppm at 4.5 $\mu$m.

In the various solutions for Spitzer transit depths described above, differences persist between 3.6 and 4.5 $\mu$m, and between the two transits observed at 3.6 $\mu$m. Those differences are minimized by our default solutions, i.e. freezing the orbital parameters at the values given by West et al. (2016) and solving for $R_{p}^{2}/R_{s}^{2}$. Given that the orbital parameters we derive in this paper are closely consistent with West et al. (2016), we adopt our default solutions for transit depths. Those values are listed in Table 3, and the best-fit transit times are included. Figure 5 illustrates the fits, after removal of the systematic effects, and binning the data for visual clarity.

The WASP-76b spectrum shows a steep slope in the G430L spectrum. In other cooler hot Jupiters, the STIS blue part of the spectrum has been used to probe the Rayleigh scattering in the atmosphere as it usually exhibits larger transit depths and slopes down into longer wavelengths. However, 0.3 to 0.4 $\mu m$ of WASP-76b spectrum shows a much steeper slope compared to the rest of the spectrum which means one continuous Rayleigh scattering slope can not sufficiently explain the observed spectrum. To understand the origin of unexpected excess transit depth, we performed retrieval analysis with ATMO (Amundsen et al., 2014; Drummond et al., 2016; Goyal et al., 2017; Tremblin et al., 2015, 2016), which has been widely used before for retrieval analyses of transmission (Wakeford et al., 2017) and emission (Evans et al., 2017) spectra. We performed a cloud-free free-element equilibrium-chemistry retrieval with free abundance of specific species (C, O, Na, Ti, V, Fe) and a fitted TP profile (5). All other elements were varied with a single metallicity parameter. ATMO is able to fit the STIS blue part of the spectrum with a solar Fe abundance and the best-fit model has a $\chi^{2}_{\nu}$ of 1.64. However, the first observed point extending from 0.29 to 0.37 $\mu m$ is still $\sim$2$\sigma$ higher than the ATMO model.

The next modeling tool we applied is PHOENIX (Lothringer et al., 2018) atmosphere forward model, It self-consistently solves layer by layer radiative transfer assuming chemical and radiative-convective equilibrium based on the irradiation received at the top of the atmosphere from the host star (Lothringer & Barman, 2019). PHOENIX is equipped with a comprehensive EUV-to-FIR opacity database of atomic opacity due to their importance in modeling stellar spectra, which makes it particularly suitable on predicting ultra-hot Jupiter atmospheres in the bluer wavelengths (Lothringer et al., 2020). We generated a grid of PHOENIX models with various metallicity, heat redistribution and internal temperature. The best-fit model (Fig. 12) is at solar metallicity with a terminator temperature of 2000K which has a $\chi^{2}_{\nu}$ of 2. With the additional opacity from metals and molecules (Fe I, Fe II, Ti I, Ni I, Ca I, Ca II, and SiO) included in the PHOENIX model, it is able to fully fit the short-wavelength slope. However, it predicts larger absorption depth between the 0.4 to 0.5 $\mu m$ region which is likely due to the assumption of solar metallicity and elemental abundances. The lower than expected abundances of NUV absorbers such as TiO, V I and Fe I could be due to condensation and/or rain-out on the day-to-night terminator detected by Ehrenreich et al. (2020).

This similar feature of steep slope in the NUV has also been observed in WASP-121b (Evans et al., 2017) with Sodium Hypochlorite (SH) proposed as the missing opacity source. With more recent observations (Sing et al., 2019) with STIS E230M from 228 and 307 nm, multiple atomic lines including Mg II and Fe II have been detected and resolved in WASP-121b. This indicates neutral and ionized atomic metal lines are more likely to be the cause of the strong NUV absorption signatures in the STIS G430L spectrum. With both WASP-76b and WASP-121b showing strong NUV absorption features, neutral and ionized metals may exist in many more ultra-hot Jupiter atmospheres (Lothringer et al., 2020).

We detected TiO and H₂O in the transmission spectrum of WASP-76b. The 0.4 to 1 $\mu m$ part of the spectrum where TiO opacity dominates shows significantly deeper transit depth ($\sim$500 ppm) compared to the WFC3/G141 spectrum. This feature is well explained by all two models with TiO absorption features. At this temperature range, TiO is expected to be in gaseous and abundant in the atmospheres as shown in the top-right panel of Fig. 15. Water vapor absorption feature at 1.4 $\mu m$ has also been observed in the spectrum, which is expected as thermal dissociation of water starts at temperatures $>$2500K (Bottom-left Fig. 15).

WASP-76b shows blackbody-like WFC3/G141 emission spectrum with muted water features but strong CO emission feature at Spitzer 4.5 $\mu m$ band. The best-fit PHOENIX model shows dayside heat-redistribution and solar metallicity assuming equilibrium chemistry. We also ran a comparison PHOENIX model with the same setup but excluding TiO/VO to demonstrate the data strongly favors the presence of gaseous TiO/VO, as the $\chi^{2}_{\nu}$ is larger by 4.32, which is consistent with our finding in the transmission spectrum. In addition, we performed ATMO free-element equilibrium chemistry retrieval similarly to the transmission spectrum, though isotopic scattering was also included along with the thermal emission. The resulting ATMO best-fit model is highly consistent compared to the PHOENIX model with both models showing similar emission spectra and TP profiles (See Fig. 14). ATMO also favors solar metallicity in the retrieval posterior distribution (Fig. 18) but with less certainty at the C/O ratio since the muted water feature limits the constrains on the oxygen abundance. Both models favor a dayside temperature range of 2500 to 2600K around 1 bar and an inverted TP profile with temperature increasing to around 3000K at 0.1 mbars. Water starts to dissociate at such high temperature and low pressure region of the atmosphere as shown in Figure 15, therefore we do not see prominent water emission features. At deeper levels ($\sim$1 bar) of the atmosphere, water vapor should still survive, but any absorption features will be obscured by the hotter continuum emission in the upper atmosphere layers. On the other hand, CO is able to survive in much higher altitude and temperature due to the strong triple bond structure. Indeed, we see clear CO emission features in the Spitzer 4.5 $\mu m$ band.

We found clear temperature inversion (Figure 14) confirmed by ATMO and PHOENIX models. The Spitzer 4.5 $\mu m$ CO emission feature strongly favors an inverted TP profile with higher temperature CO gas presence in the upper atmospheres. The transmission spectrum also favors an inverted TP profile as the retrievals need the higher temperature at the low pressures to boost the scale heights and the size of spectral features to better match the data. Theories have indicated inversion is caused by a combination of optical absorbers such as TiO (Hubeny et al., 2017; Fortney et al., 2008) and atomic metal absorption heating the upper layers with the lack of cooling from molecules like water (Lothringer et al., 2018; Gandhi & Madhusudhan, 2019). Our observed spectrum supports this paradigm with detection of TiO and atomic metal opacity in the transmission and muted water emission feature due to thermal dissociation at the highest altitudes. The detection of TiO and temperature inversion in the emission spectrum is also consistent with the independent analysis from Edwards et al. (2020) which reported similar findings.

ATMO best-fit model has the lower $\chi^{2}_{\nu}$ but PHOENIX generates a remarkable good fit in both transit and eclipse especially as only two parameters were varied in our grid of forward models. Retrieval frameworks find the best-fit spectrum through minimizing the likelihood which allows it to fine tune model parameters and better respond to smaller features in the data. Therefore, despite using an incomplete NUV opacity database, ATMO is able to produce better overall $\chi^{2}_{\nu}$ best-fit spectrum than the PHOENIX forward model. However, it is more important that ATMO and PHOENIX show good agreement on the general physical parameters including temperature structure, C/O ratio and chemical abundance. This gives us increased confidence in our conclusion, as both the retrieval and forward modeling methods agree.