HST PanCET Program: A Complete Near-UV to Infrared Transmission Spectrum for the Hot Jupiter WASP-79b

WASP-79b was observed during three primary transit events with HST STIS, two with the G430L grating and one with the G750L grating. The specific observing dates and instrument settings are summarized in Table 1. Combined, the two gratings cover the wavelength regime from 2900 to 10270 Å, with an overlapping region from ∼5260 to 5700 Å. The two gratings have a resolving power of ∼2.7 and ∼4.9 Å per pixel for the G430L and G750L gratings, respectively. Thus, they offer a resolution of R = λ/Δλ = 500–1000.

Each transit event consists of 57–77 spectra, spanning five HST orbits, where each HST orbit takes about ∼95 minutes. Because HST is in a low-Earth orbit, the data collection are truncated for ∼45 minutes in each orbit when HST is occulted by the Earth. The observations were scheduled such that the transit event occurs in the third and fourth HST orbits, while the remaining orbits provide an out-of-transit baseline before and after each transit. All observations were made with the 52 × 2 arcsec² slit to minimize slit losses. Readout times were reduced by only reading out a 1024 × 128 pixel subarray of the CCD. This strategy has previously been found to deliver high signal-to-noise ratios (S/Ns) near the Poisson limit during the transit events (e.g., Huitson et al. 2012; Sing et al. 2013). We show an example G430L and G750L spectra of WASP-79 in Figure 1.

Zoom In Zoom Out Reset image size

The relatively long exposure times (149–207 s) meant that our images were contaminated by multiple cosmic rays. Similar to previous studies, we found that correcting for cosmic rays using the CALSTIS ¹² pipeline did not yield satisfactory results. We, therefore, performed a custom cosmic-ray correction procedure based largely on the method described by Nikolov et al. (2014), which we explain here. For all the .flt images to be corrected, we created four difference images between the image itself and its two neighboring images on both sides in time. This effectively canceled out the stellar flux and left only the cosmic rays, which were seen as positive values for the image we were correcting and negative values for the neighboring image. Next, we created a median combined image from the four difference images, leaving only the cosmic-ray events that we sought to identify and replace. For all pixels in the median image, we then computed the median of the 20 closest pixels in that row and flagged the pixel in question if it exceeded a 4σ threshold in that window. When all pixels in an image were analyzed, we replaced all the flagged pixels with a corresponding value obtained from the four nearest images in time. All pixels flagged as "bad" by CALSTIS in the corresponding data quality frames were replaced in the same manner. Additionally, we inspected the extracted (see Section 2.3) 1D spectra for any potential cosmic-ray hits missed by the procedure applied on the .flt images. This was done by comparing every pixel with the corresponding value in all of the other spectra and flagging values more than 5σ above the median of that pixel. We found that a few cosmic-ray hits still persisted, which further investigations revealed to be located primarily close to the peak of the stellar point-spread function in the .flt images. These were corrected by replacing them in the same manner as before, by using the four nearest images in time.

We performed a uniform data reduction for all the STIS data. The data was bias-, flat-, and dark-corrected using the latest version of CALSTIS v3.4 and the associated relevant calibration frames.

One-dimensional spectra were extracted from the calibrated and corrected .flt science frames using the APALL procedure in IRAF. To determine what aperture size to use when running the APALL procedure, a number of different widths were tested, ranging from 9 to 17 pixels with a step size of 2. The aperture that provided the smallest out-of-transit baseline flux photometric scatter was then chosen. For all data sets, we found that this was achieved with an aperture of width 13. Like previous studies (e.g., Sing et al. 2013), no background subtraction was used as the background contribution is known to have a negligible effect. Ignoring the background can even help minimize the out-of-transit residual scatter (Sing et al. 2011; Nikolov et al. 2015). The extracted spectra were then mapped to a wavelength solution obtained from the .x1d files. The discrepancy between exposure times for the G750L visit during the last HST orbit (exposure times of 149 s) compared to the preceding orbits (exposure times of 150 s) was corrected by extrapolating for the missing second under the assumption that the detector was still operating within its linearity regime.

To allow for analyses to be done in a fully Bayesian framework, all fits were carried out by (1) treating each light curve as a Gaussian process (GP), which we implemented through the use of the GP Python package george (Ambikasaran et al. 2015); and (2) using the Nested Sampling (NS; Skilling et al. 2004) algorithm MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019), implemented by the Python package PyMultiNest (Buchner et al. 2014), which we combined with the GP likelihood function to conduct parameter inference.

GPs have been widely applied by the exoplanet community in recent years. Common applications include modeling stellar activity signals in radial velocity data (e.g., Rajpaul et al. 2015; Jones et al. 2017) and the correction of instrumentally induced systematics in transit data (e.g., Gibson et al. 2012; Sedaghati et al. 2017; Evans et al. 2018). GPs offer a nonparametric approach that finds a distribution over all possible functions that are consistent with the observed data. The main assumption behind GPs is that the input comes from infinite-dimensional data where we have observed some finite-dimensional subset of that data and this subset then follows a multivariate normal distribution. This yields the key result that input data that lie close together in input space will also produce outputs that are close together. Formally, a GP is fully defined by a mean function and a kernel (covariance) function, and it is the kernel that determines the similarity of the inputs and how correlated the corresponding outputs are. Given the rapidly growing and already extensive use of GPs applied in the literature, we refer readers unfamiliar with GPs and their applications to this type of analysis to Gibson et al. (2012), which gives a good introduction to their uses in transmission spectroscopy data.

We adopted the analytic transit model of Mandel & Agol (2002) for our GP mean function. This model is a function of midtransit times (t₀), the orbital period (P), the planet–star radius ratio (R_p /R_⋆), the semimajor axis in units of stellar radii (a/R_⋆), the orbital inclination (i), and limb-darkening coefficients. We used the BATMAN package to implement the transit model (Kreidberg 2015). Uncertainties for each data point were initially derived based solely on Poisson statistics.

Nested sampling is a numerical method for Bayesian computation targeted at the efficient calculation of the Bayesian evidence, with posterior samples produced as a by-product. Compared to traditional Markov Chain Monte Carlo (MCMC) techniques, NS is able to sample from multimodal and degenerate posteriors efficiently. It is a Monte Carlo algorithm that explores the posterior distribution by initially selecting a set of samples from the prior, called live points. The live points are then iteratively updated by calculating their individual likelihoods and replacing the live point with the lowest likelihood. This procedure ensures an increasing likelihood as the prior volume shrinks through each iteration and runs until a specified tolerance level is achieved.

3.1.1. Model Comparison

Our overall approach offers several advantages, some of which we briefly highlight here. Rather than enforcing a parameterized function to model the systematic effects, the GP allows for a nonparametric approach that simultaneously fits for both the transit and systematics. Picking an optimal systematics model or, in our case, an optimal kernel function requires conducting a model comparison. It is a general problem that such optimization routines can lead to overfitting. This problem has well-established solutions, such as the Bayesian information criterion (BIC; Schwarz 1978) and the Akaike information criterion (AIC; Akaike 1974). These two corrective terms to maximum likelihoods both enforce a penalty, based on the number of model parameters, but are slightly different in the way they introduce the penalty term. One potential critical flaw of the BIC and AIC approaches is that the model selection is based on a single maximum likelihood estimate, which does not consider the uncertainties of the model parameters, θ . Rather than relying on these methods, our application of NS allowed us to not only carry out posterior inference but also model comparison in a Bayesian framework. The model comparison was done as follows: With θ being the parameter vector, D the data, and ${ \mathcal M }$ the model, Bayes' theorem is then given by

$$\begin{eqnarray}&&p({\boldsymbol{\theta }}| {\boldsymbol{D}},{ \mathcal M })=\displaystyle \frac{p({\boldsymbol{D}}| {\boldsymbol{\theta }},{ \mathcal M })p({\boldsymbol{\theta }}| { \mathcal M })}{p({\boldsymbol{D}}| { \mathcal M })},\end{eqnarray} \tag{ 1 }$$

where p( θ ∣ D , M) is the posterior probability distribution for θ , p( D ∣ θ , M) (hereafter ${ \mathcal L }({\boldsymbol{\theta }})$) the likelihood, p( θ ∣M) (hereafter π( θ )) the prior probability, and p( D ∣M) (hereafter ${ \mathcal Z }$) is called the evidence or marginal likelihood.

${ \mathcal Z }$ is a normalization constant for the posterior and is computed from samples produced from the posterior probability distribution of θ as

$$\begin{eqnarray}&&{ \mathcal Z }=\int { \mathcal L }({\boldsymbol{\theta }})\pi ({\boldsymbol{\theta }})d{\boldsymbol{\theta }}.\end{eqnarray} \tag{ 2 }$$

It follows that the posterior probability of model ${ \mathcal M }$ is

$$\begin{eqnarray}&&p({ \mathcal M }| {\boldsymbol{D}})=\displaystyle \frac{p({\boldsymbol{D}}| { \mathcal M })p({ \mathcal M })}{p({\boldsymbol{D}})}.\end{eqnarray} \tag{ 3 }$$

To perform a relative comparison between two models we then took the ratio of the model posterior probabilities and canceling the term p( D ), yielding

$$\begin{eqnarray}&&\displaystyle \frac{p({{ \mathcal M }}_{i}| {\boldsymbol{D}})}{p({{ \mathcal M }}_{j}| {\boldsymbol{D}})}=\displaystyle \frac{p({\boldsymbol{D}}| {{ \mathcal M }}_{i})\,p({{ \mathcal M }}_{i})}{p({\boldsymbol{D}}| {{ \mathcal M }}_{j})\,p({{ \mathcal M }}_{j})}=\displaystyle \frac{{{ \mathcal Z }}_{i}\,\pi ({{ \mathcal M }}_{i})}{{{ \mathcal Z }}_{j}\,\pi ({{ \mathcal M }}_{j})}.\end{eqnarray} \tag{ 4 }$$

With no a priori model preferences the $\pi ({{ \mathcal M }}_{i})/\pi ({{ \mathcal M }}_{j})$ term cancels out, leaving us with only the evidence ratio ${{ \mathcal Z }}_{i}/{{ \mathcal Z }}_{j}$. This ratio is commonly referred to as the Bayes factor (Kass & Raftery 1995) and is what we used to directly compare two models. This model comparison comes with the benefit of incorporating Occam's razor, automatically penalizing unreasonable model complexity that in turn would lead to overfitting (and worse predictive power). Hence, our approach eliminates the need for methods such as BIC or AIC to help perform kernel choices. However, we note that the evidence calculation is based not only on the choice of kernel but also on the optimization of the hyperparameters, which can potentially get caught in bad local optima. This is usually accounted for by running the optimization routine multiple times with different starting conditions for the hyperparameters (see, e.g., the Mauna Loa atmospheric CO₂ example in chapter 5 of Rasmussen & Williams 2006), but the additional benefit of utilizing NS is its ability to handle irregular likelihood surfaces while still being efficient compared to MCMC methods. While this does not guarantee finding the global optima, we chose to rely on its ability to handle such a likelihood surface, as this provided us with a significant computational speedup.

3.1.2. Kernel Selection

With a way of evaluating the comparative performance between kernels, we set out to determine what kernel to use in the light-curve fits. We used time as the input variable and included the following five different "standard" kernels in this investigation:

• 1.  
  Squared Exponential:

  $$\begin{eqnarray}&&k({x}_{n},{x}_{m})={\sigma }^{2}\exp \left(-\displaystyle \frac{{\left({x}_{n}-{x}_{m}\right)}^{2}}{2{\ell }^{2}}\right).\end{eqnarray} \tag{ 5 }$$
• 2.  
  Rational Quadratic:

  $$\begin{eqnarray}&&k({x}_{n},{x}_{m})={\sigma }^{2}{\left(1+\displaystyle \frac{{\left({x}_{n}-{x}_{m}\right)}^{2}}{2\alpha {\ell }^{2}}\right)}^{-\alpha }.\end{eqnarray} \tag{ 6 }$$
• 3.  
  Mátern 3/2:

  $$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}k({x}_{n},{x}_{m}) & = & {\sigma }^{2}\left(1+\displaystyle \frac{\sqrt{3| {x}_{n}-{x}_{m}{| }^{2}}}{\ell }\right)\\ & & \times \exp \left(-\displaystyle \frac{\sqrt{3| {x}_{n}-{x}_{m}{| }^{2}}}{\ell }\right).\end{array}\end{array}\end{eqnarray} \tag{ 7 }$$
• 4.  
  Periodic:

  $$\begin{eqnarray}&&k({x}_{n},{x}_{m})={\sigma }^{2}\exp \left(-\displaystyle \frac{2{\sin }^{2}(\pi | {x}_{n}-{x}_{m}| /p)}{{\ell }^{2}}\right).\end{eqnarray} \tag{ 8 }$$
• 5.  
  Linear:

  $$\begin{eqnarray}&&k({x}_{n},{x}_{m})={\sigma }^{2}\left({x}_{n}-c\right)\left({x}_{m}-c\right),\end{eqnarray} \tag{ 9 }$$

where x_n , x_m refer to the elements in the covariance matrix, σ is the maximum variance allowed, ℓ is the characteristic length scales, α is the Gamma distribution parameter, p is the period between repetitions, and c is a constant term. In addition, we also incorporated a white-noise kernel in all fits, which has the form:

$$\begin{eqnarray}&&\mathrm{White}\ \mathrm{Noise}:\ k({x}_{n},{x}_{m})={\sigma }_{{wn}}^{2}\delta ({x}_{n},{x}_{m}),\end{eqnarray} \tag{ 10 }$$

where σ_wn is the amplitude of the white noise (i.e., photon noise) for each data point, and δ(x_n , x_m ) is the Kronecker delta function.

To decide which kernel to use, we tried all these kernels separately in the white-light-curve fits (see Section 3.3) and compared them by their evidence. We then expanded upon this by utilizing the fact that any additive or multiplicative combination of these five kernels are still valid kernels. This allowed us to construct more complex kernels built from these "standard" kernels, which offered the possibility of modeling many different properties that some of the "standard" kernels would struggle with. Rather than enforcing the structural form of our kernel of choice, we performed a comprehensive test of different kernels in the white-light-curve fits. This search was carried out by setting up a grid consisting of the above-mentioned five kernels and then trying all two-component additive and multiplicative combinations. Following this step, we allowed once again for an additional "standard" kernel to be added or multiplied onto the existing two-component kernels. At this stage, we kept only the best-performing kernel and attempted to expand the remaining kernel even further. We found the evidence did not improve (slightly worsened, in fact), indicating that no (or very little) structure was left. This was further reinforced by investigating the residuals after adding the third "standard" kernel, which was well described by a normal distribution with a standard deviation similar to the photon noise.

The treatment of stellar limb-darkening effects can have a significant impact on derived transmission spectra. Optimally, these effects could be accounted for by fitting for the coefficients defining a given limb-darkening model. However, the incomplete phase coverage and the relatively low temporal sampling rate of the observations make it difficult to derive the coefficients directly from the data. Instead, we used the Limb Darkening Toolkit (LDTk) (Parviainen & Aigrain 2015) Python package, which utilizes the PHOENIX stellar models of Husser et al. (2013) to calculate limb-darkening coefficients. This procedure allowed us to fix (rather than fit) the limb-darkening coefficients in our light-curve models to theoretical values, which gave us the advantage of freely picking a limb-darkening law parameterized by a higher number of coefficients. Therefore, we made use of the nonlinear limb-darkening law described by four parameters (Claret 2003) and estimated the parameters based on the PHOENIX stellar model grid point closest to that of WASP-79. The resulting coefficients used in the light-curve fits are shown in Table 6.

To refine system parameters for the planet, we initially performed fits for the light curves produced by a summation of the entire dispersion axis, commonly referred to as a white light curve. As several of the physical parameters of the system are wavelength independent, we unsurprisingly obtain the most precise system parameters when including the entire wavelength range as this ensured the highest possible S/N.

In accordance with common practice, we discarded all exposures from the first HST orbit and the first exposure of each subsequent orbit as they are known to suffer from unique and complex systematics arising from the telescope thermally relaxing into its new pointing position (Brown et al. 2001). We conducted the white-light-curve fit jointly for the two G430L grating visits, but separately from the G750L grating visit. Furthermore, we assumed a circular orbit (zero eccentricity), in correspondence with the results of Smalley et al. (2012), for all light-curve fits. This enabled us to perform our fits by allowing t₀, R_p /R_⋆, a/R_⋆, i, and the hyperparameters related to the kernel function of the GP to vary as free parameters.

As described above (see Section 3.1) we set out to find a suitable kernel for the GP, and this search resulted in a composite kernel consisting of the white-noise kernel, a periodic kernel times a squared exponential kernel (hereafter referred to as a locally periodic kernel), and a Matérn-3/2 kernel. The resulting multi-component kernel takes the form

$$\begin{eqnarray}&&\begin{array}{l}k({x}_{n},{x}_{m})={\sigma }_{{wn}}^{2}\delta ({x}_{n},{x}_{m})\\ +\ {\sigma }_{a}^{2}\exp \left(-\displaystyle \frac{2{\sin }^{2}(\pi | {x}_{n}-{x}_{m}| /p)}{{{\ell }}_{a}^{2}}\right)\exp \left(-\displaystyle \frac{| {x}_{n}-{x}_{m}{| }^{2}}{2{{\ell }}_{b}^{2}}\right)\\ +\ {\sigma }_{b}^{2}\left(1+\displaystyle \frac{\sqrt{3| {x}_{n}-{x}_{m}{| }^{2}}}{{{\ell }}_{c}}\right)\exp \left(-\displaystyle \frac{\sqrt{3| {x}_{n}-{x}_{m}{| }^{2}}}{{{\ell }}_{c}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where σ_a , σ_b and ℓ_a , ℓ_b , ℓ_c are the allowed variance and correlation length scales for the corresponding part of the composite kernel.

We note that while the search for the structural form of the kernel was determined by the data itself, the individual components of the kernel will reflect fits to physically introduced systematic effects (which could be of instrumental or astrophysical origin). Here, the locally periodic kernel component is physically motivated by the well-known breathing effect, which introduces substantial correlated systematics in the data (Brown et al. 2001). This effect is the product of the spacecraft suffering from significant thermal variations in its low-Earth ∼95 minute orbits. The Matérn-3/2 kernel component is a flexible kernel, which performs well in modeling correlations on shorter length scales, and thus is implemented to deal with residual correlated systematic effects of unknown origin. The fits assume uniform priors on all parameters. For the free transit parameters (i.e., the parameters of the GP mean function) we applied a bound on our uniform priors at ±20σ from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) inferred values of Sotzen et al. (2020). For the GP kernel parameters, we set the length-scale prior lower limit at zero and the upper limit at the value corresponding to the time between the first and the last observation, and the amplitude parameters were only restricted to not be larger than the difference between the minimum and maximum flux measurements.

The results of these fits are summarized in Table 2 and visualized in Figure 2. We found that our inferred wavelength-independent system parameters resulted in better fits as well as lower standard deviations than those quoted in the discovery paper of Smalley et al. (2012). We, therefore, chose to use these and the TESS photometry inferred values to calculate weighted average values (also shown in Table 2), which we use in the wavelength-binned fits.

Zoom In Zoom Out Reset image size

In order to assemble the transmission spectrum, we extracted light curves from wavelength bins for both gratings. Specifically, we produced wavelength-binned light curves by dividing the spectra into bins varying in size from 85 to 1000 Å. We chose to customize bin sizes based on the criterion that the S/N in each bandpass was sufficiently high not to be dominated by photon noise, yet still small enough to preserve valuable information from the underlying transmission spectra. This criterion was achieved in all bins with an average S/N of ∼1500. Additionally, we also made sure that the borders of the channels did not coincide with prominent stellar lines. Fits were then carried out in each spectrophotometric channel similar to the white-light-curve fits (see Section 3.3), but with the exception that we froze each wavelength-independent parameter to the weighted average values quoted in Table 2. Effectively, this meant that the fits carried out in the spectrophotometric channels only allowed for the parameter of interest, R_p/R_⋆, and the GP kernel parameters to vary as free parameters. Identical to the white-light-curve fits, we jointly fit the wavelength-binned light curves from our two data sets obtained with the G430L grating.

We accounted for potential systematic discrepancies in the absolute transit depths between the G430L and G750L gratings and put the combined STIS transmission spectrum on an absolute scale. First, we measured the offset between the two gratings using their overlapping region between 0.53 and 0.57 μm. We did this by fitting the light-curve jointly for the two G430L data sets and separately for the G750L data set in the overlapping region. Second, TESS observed 12 transits of WASP-79b in January and February of 2019, yielding a tight constraint of R_p/R_⋆ = 0.10675 ± 0.00014 (Sotzen et al. 2020), that we utilized to calibrate the transmission spectrum to an absolute scale. Therefore, we performed a similar fit for the G750L grating in the 0.59–1.02 μm range corresponding to the TESS bandpass. Finally, we stitched the combined transmission spectra together by uniformly offsetting the G430L transmission spectra by the difference between the inferred values for the two gratings in the same bandpass, and then offsetting the entire transmission spectrum in the same way, anchoring it to the TESS value. The inferred values are noted in Table 3. The detrended binned light curves for all three visits are shown in Figure 3 for visits 67 and 68, and Figure 4 for visit 69. To check for the sensitivity of our limb-darkening treatment, we repeated the analysis but applied the quadratic limb-darkening law instead and fit for the two coefficients in the light-curve models. We found the two treatments showed excellent consistency in the relative transit depths, and measurements in all channels agreed within 1σ (see the Appendix, Figures 11 and 12). Our final stitching-corrected transmission spectrum is presented in Figure 5 (alongside the observations from Sotzen et al. 2020) and summarized in Table 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A visual inspection of the transmission spectrum shows no obvious signs of absorption from sodium or potassium, but it does show decreasing transit depths toward blue wavelengths over the optical spectral range. This morphology also appears in ground-based LDSS3 data from Sotzen et al. (2020), though those data appear systematically vertically offset from the STIS data. It is not clear what is causing this vertical offset, but some explanations include instrumental systematics, stellar variability, different orbital parameters, and/or different limb-darkening coefficients. Nevertheless, we verified that excluding the LDSS3 data does not alter our atmospheric inferences in later sections.

We conducted a series of atmospheric retrievals using three different codes. Each code was free to choose its own set of molecular, atomic, and ionic opacities, along with a pressure–temperature (P-T) profile and cloud/haze parameterization. Our approach, considering multiple independent retrieval codes, ensures robust atmospheric inferences. The configurations used by each retrieval code are summarized in Table 4.

Our retrievals include chemical species with prominent spectral features over the observed wavelength range (Sharp & Burrows 2007; Tennyson & Yurchenko 2018) anticipated to be present in hot Jupiter atmospheres (Madhusudhan et al. 2016). For the near-infrared, we assess contributions from H₂O, CH₄, CO, CO₂, HCN, and NH₃. For optical wavelengths, we consider Na, K, H⁻, TiO, VO, FeH, and Fe. Each retrieval considered a subset of these potential species. Common opacity across all three codes are H₂O, H⁻, collision-induced absorption due to H₂–H₂ and H₂–He (Richard et al. 2012), and H₂ Rayleigh scattering.

Multiple studies have recently considered the inclusion of H⁻ opacity in atmospheric retrievals (e.g., Gandhi et al. 2020; Lothringer & Barman 2020; Sotzen et al. 2020). However, there remains no consensus on how to parameterize H⁻ opacity in a retrieval context. H⁻ is expected to become an important opacity source at high temperatures, when H₂ thermally dissociates to form atomic H (e.g., Bell & Cowan 2018; Parmentier et al. 2018). Atomic H absorbing a photon in the vicinity of a free electron produces free–free H⁻ absorption (Bell & Berrington 1987):

$$\begin{eqnarray}&&h\nu +{{\rm{e}}}^{-}+{\rm{H}}\to {\rm{H}}+{{\rm{e}}}^{-}.\end{eqnarray} \tag{ 12 }$$

Alternately, the photodissociation of a bound H⁻ ion results in bound–free H⁻ absorption (John 1988):

$$\begin{eqnarray}&&h\nu +{{\rm{H}}}^{-}\to {\rm{H}}+{{\rm{e}}}^{-}.\end{eqnarray} \tag{ 13 }$$

Although both processes are commonly referred to as "H⁻ opacity," only the bound–free contribution involves a H⁻ ion. Due to the distinct nature of these processes, we use a general treatment to parameterize H^– opacity. Considering their combined opacity,

$$\begin{eqnarray}&&{\kappa }_{{{\rm{H}}}^{-}}={n}_{{{\rm{H}}}^{-}}\,{\sigma }_{\mathrm{bf},{{\rm{H}}}^{-}}+{n}_{{\rm{H}}}\,{n}_{{{\rm{e}}}^{-}}\,{\alpha }_{\mathrm{ff},{{\rm{H}}}^{-}},\end{eqnarray} \tag{ 14 }$$

where ${\kappa }_{{{\rm{H}}}^{-}}$ is the H⁻ extinction coefficient (in cm⁻¹), ${\sigma }_{\mathrm{bf},{{\rm{H}}}^{-}}$ is the bound–free H⁻ cross section (given in cm² by Equations (4) and (5) from John 1988), ${\alpha }_{\mathrm{ff},{{\rm{H}}}^{-}}$ is the free–free H⁻ binary cross section (given in cm⁵ by Equation (6) from John 1988 multiplied by k_B T in cgs units), and n_i are the number densities (in cm⁻³) of H⁻, H, and free electrons. We propose that a parameterization suitable for free retrievals is to treat the mixing ratios of H⁻, H, and e⁻ as independent free parameters.

We also included the effects of unocculted spot/faculae in all our retrievals (e.g., Pinhas et al. 2018; Rackham et al. 2018). This was motivated by the negative slope toward blue wavelengths in our transmission spectrum (Figure 5)—atmospheric scattering would instead cause a positive slope—and indicators of stellar activity for WASP-79 (Section 4). We adopt a consistent prescription for stellar heterogeneity across all three codes. This invokes a three-parameter model, based on the approach of Pinhas et al. (2018):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda ,\mathrm{obs}}={{\rm{\Delta }}}_{\lambda ,\mathrm{atm}}\,{\epsilon }_{\lambda ,\mathrm{het}},\end{eqnarray} \tag{ 15 }$$

where Δ_λ , _obs is the observed transmission spectrum, Δ_λ , _atm is the transmission spectrum from the planetary atmosphere alone, and _λ , _het is the wavelength-dependent "contamination factor" from a heterogeneous stellar surface. For a two-component stellar disk with a photosphere and an excess heterogeneity (spots or faculae), the contamination factor can be written as (Rackham et al. 2018)

$$\begin{eqnarray}&&{\epsilon }_{\lambda ,\mathrm{het}}={\left(1-{f}_{\mathrm{het}}\left(1-\displaystyle \frac{{I}_{\lambda ,\mathrm{het}}({T}_{* ,\mathrm{het}})}{{I}_{\lambda ,\mathrm{phot}}({T}_{* ,\mathrm{phot}})}\right)\right)}^{-1},\end{eqnarray} \tag{ 16 }$$

where f_het is the fractional stellar disk coverage of the heterogeneous regions, and I_λ , _het and I_λ , _phot are the specific intensities of the heterogeneity and photosphere, respectively, with T_*,het and T_*,phot their corresponding temperatures. In our default retrieval prescription, we treat f_het, T_*,het, and T_*,phot as free parameters. We also investigated replacing the T_*,het parameter with the average temperature difference between heterogeneous regions and the photosphere: ΔT_* = T_*,het − T_*,phot. Because the stellar photosphere temperature is known a priori, we place an informative Gaussian prior on T_*,phot. We compute stellar spectra by interpolating models from the Castelli-Kurucz 2004 atlas (Castelli & Kurucz 2003) using the pysynphot package (STScI Development Team 2013).

All three codes conducted an atmospheric retrieval, as summarized in Table 4, for parameter estimation. The full posterior distributions from these retrievals are provided as supplementary material. Detection significances for key model components were also computed, via Bayesian model comparisons. We now provide a brief summary of each retrieval code.

5.1.1. NEMESIS

5.1.2. POSEIDON

5.1.3. ATMO

Here we present the results of our comparative retrievals. We first explore the best-fitting atmospheric and stellar interpretation matching the transmission spectrum of WASP-79b. Constraints on the atmospheric properties of WASP-79b are then presented, followed by inferences of stellar heterogeneity.

5.2.1. Explaining the Transmission Spectrum of WASP-79b

Our retrievals arrived at a consistent explanation for the transmission spectrum of WASP-79b. Our best-fitting model spectra, shown in Figure 7, are characterized by three components: (i) H₂O opacity (explaining the absorption feature around 1.4 μm), (ii) spectral contamination from unocculted faculae (producing the negative slope over optical wavelengths), and (iii) H⁻ bound–free absorption (resulting in a relatively smooth continuum from 0.4 to 1.3 μm). Our new STIS observations play a crucial role in the identification of faculae, extending the coverage of WASP-79b's transmission spectrum to wavelengths <0.6 μm where the effects of stellar contamination are more pronounced. We verified that retrievals excluding the LDSS3 data (i.e., STIS + WFC3 + Spitzer only) arrive at the same conclusion.

Zoom In Zoom Out Reset image size

The broad agreement between our retrievals, despite their quite different configurations (Table 4), motivated an exercise to identify the minimal model capable of explaining the present observations. Although the fit qualities shown in Figure 7 are comparable, the differing numbers of free parameters (30 for POSEIDON, 21 for NEMESIS, and 9 for ATMO) resulted in a range of best-fitting reduced chi-square values suggestive of model overcomplexity for the present data sets (${\chi }_{\nu ,\min }^{2}$ = 1.84, 1.68, and 1.25 for POSEIDON, NEMESIS, and ATMO, respectively). Taking the nine parameters defining the ATMO model as a starting point, we ran additional retrievals with progressively fewer free parameters to identify the simplest model capable of explaining WASP-79b's transmission spectrum—corresponding to the model with maximal Bayesian evidence (analogous to ${\chi }_{\nu ,\min }^{2}$ minimization). This process arrived at a seven-parameter "minimal" model: ¹³ H₂O and H⁻ in a clear, isothermal, H₂-dominated atmosphere transiting a stellar surface with unocculted faculae (${\chi }_{r,\min }^{2}$ = 1.12). We show the best-fitting spectrum from this retrieval in Figure 8, demonstrating consistency with the more complex models explored previously. With respect to this minimal model, Bayesian model comparisons yielded strong detections of faculae (4.7σ) and H₂O (4.0σ), along with moderate evidence of H⁻ (3.3σ).

Zoom In Zoom Out Reset image size

The contributions of these opacity/contamination sources to the best-fitting minimal model spectrum are shown in Figure 8. The features of the observed spectrum are reproduced by a combination of H₂O, H⁻ and H₂–H₂ collision-induced absorption within the planet's atmosphere, alongside contributions from unocculted faculae on the stellar surface. Faculae occupying regions of the star outside the transit chord result in the planet occulting a region of the star that is cooler and darker than the disk average, since faculae are relatively hot and bright. This results in an underestimation of the true planet-to-star radius ratio, as illustrated by the "atmosphere-only" model in Figure 8, which has greater transit depths than the composite spectrum. The magnitude of this effect varies with wavelength, such that the faculae/disk contrast is most pronounced at shorter wavelengths, resulting in shallower transit depths in the near-UV and optical relative to the infrared, thus giving rise to the negative slope at short wavelengths. This "transit light source effect" is discussed more generally in Rackham et al. (2018).

5.2.2. The Atmosphere of WASP-79b

Our retrieved atmospheric properties are displayed in Figure 9 (top row) and summarized in Table 5. All three codes reach excellent agreement on the atmospheric parameters of WASP-79b, which all agree within their respective 1σ confidence intervals.

Zoom In Zoom Out Reset image size

Table 5. Retrieved Planetary and Stellar Parameters

Retrieval	POSEIDON	NEMESIS	ATMO	Minimal
Planetary Atmosphere
T1 mbar (K)	${958}_{-168}^{+217}$	${936}_{-233}^{+272}$	${1028}_{-196}^{+303}$	${836}_{-139}^{+187}$
Rp,ref (RJ )	${1.68}_{-0.02}^{+0.02}$	${1.62}_{-0.03}^{+0.02}$	${1.75}_{-0.02}^{+0.01}$	${1.70}_{-0.02}^{+0.01}$
log(${X}_{{{\rm{H}}}_{2}\ {\rm{O}}}$)	$-{1.92}_{-0.53}^{+0.39}$	$-{2.37}_{-1.19}^{+0.64}$	$-{2.43}_{-1.03}^{+0.66}$	$-{2.43}_{-0.96}^{+0.66}$
log(${X}_{{{\rm{H}}}^{-}}$)	$-{7.80}_{-0.63}^{+0.66}$	$-{7.72}_{-0.79}^{+0.78}$	$-{8.37}_{-0.78}^{+0.98}$	$-{8.81}_{-0.65}^{+0.78}$
Derived Properties
O/H (×stellar)	${13.9}_{-9.7}^{+19.9}$	${4.7}_{-4.4}^{+16.0}$	${4.1}_{-3.8}^{+14.8}$	${4.1}_{-3.7}^{+14.8}$
Stellar Properties
T*,phot (K)	${6619}_{-80}^{+84}$	${6992}_{-271}^{+281}$	${6616}_{-71}^{+67}$	${6623}_{-92}^{+87}$
ΔT* (K)	${486}_{-173}^{+270}$	${480}_{-128}^{+153}$	${495}_{-192}^{+390}$	${452}_{-237}^{+605}$
fhet	${0.15}_{-0.07}^{+0.09}$	${0.17}_{-0.05}^{+0.07}$	${0.14}_{-0.07}^{+0.11}$	${0.15}_{-0.10}^{+0.20}$
Statistics
${\chi }_{\nu ,\min }^{2}$	1.84	1.68	1.25	1.12
Nparam	30	21	9	7
Degrees of freedom	31	40	52	54

Note. Only parameters with bounded constraints (i.e., both lower and upper bounds) are included in this summary table—see the supplementary material for the full posterior distributions. The "minimal" retrieval is the simplest model that can fit our transmission spectra of WASP-79b: a clear, isothermal, atmosphere containing H₂O and H⁻ alongside stellar contamination from unocculted faculae (seven free parameters)—also computed with POSEIDON. R_p,ref is defined at P = 10 bar for NEMESIS and POSEIDON, and 1 mbar for ATMO. ΔT_* = T_*,het − T_*,phot. The stellar O/H is assumed equal to WASP-79's stellar [Fe/H] (0.03, Stassun et al. 2017).

Download table as:  ASCIITypeset image

Our robust detection of H₂O allows the atmospheric metallicity of WASP-79b to be constrained. We derive statistical constraints on the atmospheric O/H ratio from the full set of posterior samples as in MacDonald & Madhusudhan (2019). Taking the stellar [Fe/H] ( = 0.03, Stassun et al. 2017) to be representative of the stellar [O/H], we compute the metallicity of WASP-79b (relative to its host star) via M = (O/H)_atm/(O/H)_*. Our retrieved metallicities are then as follows: ${13.9}_{-9.7}^{+19.9}\times$ stellar (POSEIDON), ${4.7}_{-4.4}^{+16.0}\times$ stellar (NEMESIS), and ${4.1}_{-3.8}^{+14.8}\times$ stellar (ATMO). The median retrieved values are suggestive of a superstellar metallicity for WASP-79b. However, a stellar metallicity remains consistent with the present observations to 1σ for NEMESIS and ATMO (due to the long tails in their H₂O abundance posteriors) and to 2σ for POSEIDON.

The bound–free absorption of H⁻ inferred from our retrievals produces corresponding constraints on its abundance. All three retrievals concur on a H⁻ abundance of log(${X}_{{{\rm{H}}}^{-}}$) ≈ −8.0 ± 0.7 (see Table 5). This precise H⁻ constraint arises from two principal features of our observations: (i) the high-precision WFC3 G141 data (∼50 ppm) closely follow the shape of the H⁻ bound–free opacity near the photodissociation limit (see Figure 8); and (ii) the long spectral baseline provided by our STIS observations, alleviating normalization degeneracies (see, e.g., Heng & Kitzmann 2017; Welbanks & Madhusudhan 2019). We discuss the plausibility of our inferred H⁻ opacity in Section 6.2.

Our retrievals additionally constrain the terminator temperature of WASP-79b. Despite the three different P-T profile prescriptions (see Table 4), all three codes arrived at the same conclusion: a near-isothermal terminator with T ∼ 1000 ± 300 K. We summarize the retrieved temperatures for each code in Table 5, where for the nonisothermal profiles we quote T_{1 mbar} as a photosphere proxy. This retrieved temperature is markedly colder than the equilibrium temperature of WASP-79b (1900 ± 50 K, Smalley et al. 2012). Recently, MacDonald et al. (2020) noted that most retrieved temperatures from transmission spectra are significantly colder than T_eq. They attributed this trend to a bias arising from 1D atmosphere assumptions. We discuss the impact of this bias on our retrieved metallicity in Section 6.3.

The atmospheric region probed by our transmission spectrum is consistent with a lack of detectable cloud opacity. However, the limits on cloud properties derived by our retrievals (e.g., P_cloud > 10⁻⁵ bar to 2σ from POSEIDON) still allow for the existence of deeper cloud decks. Nevertheless, for the present observations, our results are invariant to the chosen cloud prescription.

5.2.3. Stellar Properties

Our retrieved stellar properties are shown in Figure 9 (bottom row) and also summarized in Table 5. All three retrievals produce a consistent interpretation: a stellar surface with ∼15% faculae coverage at a temperature contrast of ∼500 K. The only disagreement found is the retrieved photospheric temperature, for which NEMESIS retrieves a value ∼400 K higher than ATMO and POSEIDON. This difference arises from the Gaussian prior used for T_*,phot by POSEIDON and ATMO having a standard deviation one-fifth that used by NEMESIS (100 K versus 500 K). However, this discrepancy does not influence the agreement for any of the other retrieved parameters, as it is the temperature contrast, ΔT_*, that governs the manifestation of spectral contamination from heterogeneous regions.

The degeneracies between our retrieved stellar parameters are shown in Figure 10. The fractional faculae coverage, f_het, and faculae temperature, T_*,het, are partially degenerate. The origin of this degeneracy is intuitive: a large coverage fraction with relatively cool faculae produces a similar contamination signal to a lower coverage fraction with hot faculae. The uncertainties in other parameters introduced by this degeneracy are already accounted for in the marginalized posteriors shown in Figure 9. The f_het − T_*,het degeneracy is not, however, an exact degeneracy. For sufficiently large ΔT_*, wavelength-dependent spectral signatures arise (from the intensity ratio in Equation (16)) that cannot be compensated by varying f_het (see, e.g., Pinhas et al. 2018). Such signatures are especially prominent at the short wavelengths sampled by our STIS G430L observations (0.3–0.6 μm; see Figure 7). This second-order effect, crucially probed by STIS, allows bounded constraints on the faculae coverage fraction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (5 panels)

Tarball of figure set images

An optical to infrared transmission spectrum of WASP-79b, combining results from Hubble/WFC3 with optical data from LDSS3 and photometry from TESS and Spitzer, has been published by Sotzen et al. (2020). A retrieval analysis on this data set was performed using the ATMO retrieval code. Our present study extends the blue wavelength coverage of the transmission spectrum of WASP-79b from 0.6 μm down to 0.3 μm with Hubble/STIS G430L data while adding STIS G750L data to complement the prior LDSS3 observations.

The ATMO free retrieval presented by Sotzen et al. (2020) considered many similar atmospheric parameters to those used here, including H₂O, CO, Na, K, VO, FeH, and H⁻. Our retrievals investigated a wider range of chemical species, cloud, and temperature profile parameterizations, and a stellar heterogeneity treatment (see Table 4). Despite our different retrieval prescriptions and expanded data set, we obtain consistent H₂O abundances and temperatures, and agree on the overall lack of significant cloud opacity. However, differences emerge when considering other gaseous species. Sotzen et al. (2020) infer the presence of Na, although they stress that this is driven by the TESS photometry point having a deeper transit than the LDSS3 data, rather than a resolved Na profile. We do not find evidence for Na. The other key difference is that Sotzen et al. (2020) do not detect H⁻ but instead find evidence for FeH.

These differences can be attributed to the additional information provided by our STIS observations. The combination of the WFC3 spectrum with the broad, relatively flat, STIS spectra leads our retrievals to favor H⁻ as the optical absorber (modulated by unocculted faculae) rather than FeH. While both FeH and H⁻ are capable of fitting the short-wavelength end of the WFC3 spectrum, FeH absorption does not extend shortwards of ∼0.7 μm (see Figure 14 in Sotzen et al. 2020 and Tennyson & Yurchenko 2018). The importance of H⁻ is seen clearly in Figure 8 and is identified by all three retrieval codes, as shown in Figure 7. We performed an additional POSEIDON retrieval without H⁻ opacity to investigate the differences between our results and Sotzen et al. (2020). In this case, FeH is recovered with a large abundance mode (log(FeH) = −2.8 ± 0.5) consistent with that found by Sotzen et al. (2020; see the supplementary material). However, our retrievals considering both H⁻ and FeH rule out such high FeH abundances (log(FeH) < −3.78 to 2σ) and have a higher Bayesian evidence. We also note that FeH abundances exceeding ∼10⁻⁷ are unexpected in thermochemical equilibrium for a giant planet at WASP-79b's equilibrium temperature (Visscher et al. 2010), and FeH abundances >10⁻⁶ were recently ruled out by high-resolution transmission spectra of 12 giant exoplanets (Kesseli et al. 2020). In summary, we find that the combined effect of H⁻ and unocculted faculae provides a more robust explanation of WASP-79b's transmission spectrum.

Recently, Skaf et al. (2020) retrieved an alternative reduction of WASP-79b's WFC3 transmission spectrum with the TauREx code. They also do not include H⁻ opacity and hence recover a similar FeH abundance to Sotzen et al. (2020). The retrieved temperature and H₂O abundance are consistent with our findings.

H⁻ is thought to become an important opacity source for high-temperature planets when H₂ thermally dissociates to form atomic hydrogen. H₂ dissociation is generally expected to occur on the day sides of "ultra-hot" Jupiters (e.g., Bell & Cowan 2018), with prominent H⁻ opacity expected around 2500 K (Arcangeli et al. 2018). Under the assumption of chemical equilibrium, H⁻ abundances decrease for lower temperatures, with observable signatures in transmission spectra not expected below ∼2100 K (Goyal et al. 2020).

Our retrieved terminator temperature for WASP-79b (∼1000 K) lies in a regime where H⁻ opacity would not be expected under equilibrium considerations. In comparison, the self-consistent models of Goyal et al. (2020) predict 1 mbar temperatures of ∼2000 K for WASP-79b (assuming C/O = 0.5, M/H = 10× solar, and recirculation factor of unity). However, our three retrieval analyses assumed a uniform (1D) composition and temperature across the terminator. Transmission spectra of planets exhibiting nonuniform compositions in the terminator region can lead to biased temperatures when subject to a 1D atmospheric retrieval (MacDonald et al. 2020; Pluriel et al. 2020). Given the strong temperature dependence of equilibrium H⁻ abundances, one would expect large H⁻ compositional gradients between different regions of the terminator. MacDonald et al. (2020) demonstrated that transmission spectra of planets for which H⁻ is only present on the warmer evening terminator result in retrieved 1D temperatures biased ∼1000 K colder than the terminator average temperature. An inter-terminator H⁻ abundance gradient therefore provides a potential explanation for the discrepancy between our retrieved temperatures, self-consistent models, and the existence of H⁻ opacity.

A further possibility is that the ∼10⁻⁸ H⁻ abundance we infer results from disequilibrium photochemistry. Lewis et al. (2020) recently showed that a similar H⁻ abundance can be produced in HAT-P-41b's atmosphere (T_eq ∼ 1900 K) by photochemical production of free electrons followed by dissociative electron attachment of H₂ (H₂ + e⁻ → H + H⁻).

An essential goal of atmospheric studies is to relate retrieved properties to planetary formation histories and environments. The abundance enhancements of elements relative to hydrogen, or metallicity, provide a crucial link to planetary formation mechanisms (e.g., Öberg et al. 2011; Madhusudhan et al. 2014; Mordasini et al. 2016). The solar system giant planets exhibit an inverse correlation between planet mass and metallicity (from C/H measurements), commonly interpreted as evidence of formation by core accretion (Pollack et al. 1996). On the other hand, many hot Jupiter exoplanets are consistent with substellar O/H ratios below the solar system mass–metallicity trend (Barstow et al. 2017; Pinhas et al. 2018; Welbanks et al. 2019). With a roughly Jovian mass (0.9 M_J), WASP-79b provides an opportunity to benchmark a hot Jupiter metallicity against elemental abundance measurements of Jupiter.

Our retrieved H₂O abundances generally suggest somewhat superstellar O/H ratios for the atmosphere of WASP-79b (∼0.3–34× stellar). This is consistent with the C/H abundance of Jupiter of ∼4× solar (Atreya et al. 2018) and recent preliminary measurements of the equatorial O/H abundance of Jupiter from JUNO (Li et al. 2020). However, our median H₂O abundances are ∼100× higher than those derived from transmission spectra of the similar-mass hot Jupiters HD 209458b and HD 189733b (Barstow et al. 2017; MacDonald & Madhusudhan 2017; Pinhas et al. 2018; Welbanks et al. 2019). We note that our retrieved H₂O abundances, and hence metallicities, may be biased by ≲1 dex toward higher abundances if a H⁻ abundance gradient exists between the morning and evening terminators (see MacDonald et al. 2020, their Figure 3). Even accounting for a factor of 10 H₂O bias, the maximum likelihood H₂O abundance for WASP-79b would remain ∼10× higher than those of HD 209458b and HD 189733b. This suggests the formation of WASP-79b may be more analogous to Jupiter than to other similar-mass hot Jupiters, indicative of a diversity of formation avenues at play across the hot Jupiter population.

Our retrievals favor models including stellar contamination, arising from unocculted faculae ≈500 K hotter than the photosphere of the host star and covering ≈15% of the stellar surface (Table 5). A similar contamination effect was also observed in the transmission spectrum of GJ 1214b (Rackham et al. 2017), though that study found ΔT ≈ 350 K and a faculae coverage fraction about five times smaller (≈3%). However, it is hard to reconcile the presence of faculae on the surface of WASP-79 with other available observations: although our XMM-Newton observations indicate a moderate level of chromospheric activity, the stellar photosphere is not expected to present large active regions—as indicated by the low photometric variations in the TESS light curves described in Section 4. However, certain geometrical configurations of the system can resolve this apparent discrepancy. As WASP-79b follows a nearly polar orbit (Addison et al. 2013), we do not know the inclination of WASP-79. This opens the possibility that we are observing the star pole on. We do not expect to see high levels of variability in the TESS data in this case, as effectively few active regions would rotate into or out of view. This scenario would explain the low-level photometric variability, while still allowing a high coverage fraction of unocculted faculae. While this remains speculative, we also note that the posteriors for our retrieved stellar parameters are broad. Consequently, a lesser degree of heterogeneity is still consistent with our transmission spectrum (e.g., ΔT ≈ 300 K temperature contrast with ≈7% coverage fraction lies within 1σ). A promising avenue for future work would be to better constrain the degree of stellar heterogeneity WASP-79 exhibits via out-of-transit stellar decomposition (e.g., Wakeford et al. 2019; Iyer & Line 2020).

WASP-79b is a shortlisted target for the JWST Transiting Exoplanet Community ERS Program (PI: Batalha, ERS 1366; Bean et al. 2018). If observed, a complete transmission spectrum will be constructed from 0.6 to 5.3 μm via observations with four instrument modes: NIRISS SOSS, NIRSpec G235H, NIRSpec G395H, and NIRCam F322W2. Our results inform the potential science return of such observations.

Consistent with previous studies (Skaf et al. 2020; Sotzen et al. 2020), we find an atmosphere with a high H₂O abundance (∼1%) and negligible cloud opacity. Our best-fitting models, therefore, predict a prominent 3 μm H₂O feature spanning ∼4 scale heights (see Figure 8), which will be readily detectable by all four JWST observations. Consequently, precise H₂O abundance and metallicity determinations (≲0.2 dex; Sotzen et al. 2020) will be possible. Spectrally resolved CO and CO₂ features around 4.5 μm with NIRSpec G395H will further allow a precise C/O ratio determination.

Our tentative inference of H⁻ opacity offers intriguing possibilities for the potential JWST observations of WASP-79b. First, NIRISS SOSS can readily assess the presence of H⁻ via precision measurements of the characteristic bound–free opacity in the optical and near-infrared. If confirmed, a high-significance H⁻ detection and abundance constraint would result. Second, the longest wavelength observations (>4 μm) with NIRSpec G395H may detect free–free H⁻ opacity, enabling one to measure the atmospheric electron mixing ratio (Lothringer & Barman 2020). JWST observations of WASP-79b, therefore, have the potential to open a window into ionic chemistry in hot Jupiter atmospheres.