Detection of Na, K, and H₂O in the hazy atmosphere of WASP-6b

ABSTRACT

We present new observations of the transmission spectrum of the hot Jupiter WASP-6b both from the ground with the Very Large Telescope FOcal Reducer and Spectrograph (FORS2) from 0.45 to 0.83 μm, and space with the Transiting Exoplanet Survey Satellite from 0.6 to 1.0 μm and the Hubble Space Telescope (HST) Wide Field Camera 3 from 1.12 to 1.65 μm. Archival data from the HST Space Telescope Imaging Spectrograph (STIS) and Spitzer are also re-analysed on a common Gaussian process framework, of which the STIS data show a good overall agreement with the overlapping FORS2 data. We also explore the effects of stellar heterogeneity on our observations and its resulting implications towards determining the atmospheric characteristics of WASP-6b. Independent of our assumptions for the level of stellar heterogeneity we detect Na i, K i, and H₂O absorption features and constrain the elemental oxygen abundance to a value of [O/H] ≃ −0.9 ± 0.3 relative to solar. In contrast, we find that the stellar heterogeneity correction can have significant effects on the retrieved distributions of the [Na/H] and [K/H] abundances, primarily through its degeneracy with the sloping optical opacity of scattering haze species within the atmosphere. Our results also show that despite this presence of haze, WASP-6b remains a favourable object for future atmospheric characterization with upcoming missions such as the James Webb Space Telescope.

1 INTRODUCTION

Transiting exoplanets currently present one of the best options towards studying the atmospheres of planets outside of the Solar system through observations of wavelength-dependent variations in their apparent radii as they occult their host star. These variations are intrinsically linked to the composition and structure of an exoplanetary atmosphere, as the starlight transmitted through the planetary limb is strongly modulated by the wavelength-dependent opacities of its constituent molecular species (Seager & Sasselov 2000). Tracing these variations as a function of wavelength, known as transmission spectroscopy, has already been successfully applied across a range of both ground- and space-based observatories, unveiling a host of atomic and molecular species in the atmospheres of exoplanets (e.g. Charbonneau et al. 2002; Redfield et al. 2008; Snellen et al. 2008; Sing et al. 2011b; Deming et al. 2013; Spake et al. 2018; Evans et al. 2018) as well as providing strong insights into their bulk atmospheric properties (e.g. Madhusudhan et al. 2011; Evans et al. 2017; Wakeford et al. 2018). In particular, Sing et al. (2016) show a large diversity in the atmospheres of a sample of 10 hot Jupiter exoplanets, revealing a continuum in the obscuring effects of haze and clouds on molecular absorption features present in their transmission spectra. Of the 10 exoplanets displayed by Sing et al. (2016), WASP-6b and WASP-39b were lacking in near-infrared observations between 1 and 2 μm, a region abundant in potential water absorption features. Wakeford et al. (2018) reported such observations for WASP-39b, providing a strong constraint on the water abundance in its atmosphere. In this study, we present these observations for WASP-6b, completing the search for water absorption features across this sample of exoplanets.

Space-based observations, such as those performed with the Hubble Space Telescope (HST) and Spitzer, have thus far proven to be the most prolific method towards the broad spectrophotometric characterization of exoplanet atmospheres (e.g. Charbonneau et al. 2002; Deming et al. 2013; Sing et al. 2016). However, ground-based characterization through multi-object differential spectrophotometry with the Very Large Telescope (VLT) FOcal Reducer and Spectrograph (FORS2) (Appenzeller et al. 1998), has recently been able to produce HST-quality transmission spectra for a variety of exoplanets (Bean et al. 2011; Nikolov et al. 2016; Gibson et al. 2017; Sedaghati et al. 2017; Nikolov et al. 2018). As part of a small survey to test the performance of FORS2 and assess the validity of previously observed spectroscopic features with HST, the optical spectra of WASP-31b, WASP-39b, and WASP-6b have been observed. In the case of WASP-39b and WASP-31b, these results have already been reported in Nikolov et al. (2016) and Gibson et al. (2017), respectively. In this study, we report the results for WASP-6b, the final target from our ground-based comparative program.

WASP-6b is an inflated hot Jupiter with a mass of |$0.485\, M_{\rm {Jup}}$|⁠, a radius of |$1.230\, R_{\rm {Jup}}$|⁠, and an equilibrium temperature of 1184 K (Tregloan-Reed et al. 2015) discovered by the Wide Angle Search for Planets (WASP) ground-based transit survey (Pollacco et al. 2006; Gillon et al. 2009). WASP-6b orbits with a period of P ≃ 3.36 d at a separation a ≃ 0.041 au around a mildly metal-poor G8V star (Gillon et al. 2009; Tregloan-Reed et al. 2015). Ibgui, Burrows & Spiegel (2010) demonstrate that the planet’s inflated radius could be due to tidal heating brought on by a non-zero eccentricity reported in Gillon et al. (2009). Whilst further radial velocity data from Husnoo et al. (2012) demonstrated that this eccentricity is not significantly non-zero, as initially inferred, it does not necessitate a circular orbit and as such the true cause of the inflation has yet to be definitively determined. Doyle et al. (2013) refine the bulk properties of the host star WASP-6 through spectroscopy, providing measurements of T_eff  = 5375 ± 65, log(g) = 4.61 ± 0.07, and [Fe/H]  = −0.15 ± 0.09. Finally, Tregloan-Reed et al. (2015) demonstrated that fluctuations in multiple transit light curves of archival photometry of WASP-6b could be attributed to a single star-spot anomaly. This enabled a more precise measurement on the sky projected spin–orbit alignment of λ = 7.2° ± 3.7° in agreement with Gillon et al. (2009).

The atmosphere of WASP-6b was initially probed spectrophotmetrically in the optical with the ground-based IMACS instrument on the 6.5-m Magellan Telescope by Jordán et al. (2013) who observed a decrease in transit depth as a function of wavelength, characteristic of a scattering haze, and no evidence of the Na i and K i absorption lines. Subsequent observations performed in the optical with HST’s Space Telescope Imaging Spectrograph (STIS) and Spitzer’s InfraRed Array Camera (IRAC, Nikolov et al. 2015) also demonstrated evidence of a scattering haze, however the Na i and K i lines were resolved in this case with significance levels of 1.2σ and 2.7σ, respectively. WASP-6b’s atmosphere has also been observed at secondary eclipse as the planet passes behind its host star from our point of view with Spitzer IRAC, providing dayside temperature estimates of |$1235_{-77}^{+70}$| and |$1118_{-74}^{68}$| K for the 3.6 and 4.5 μm channels, respectively (Kammer et al. 2015).

We present new spectrophotometric observations from 1.1 to 1.7 μm using the HST Wide Field Camera 3 (WFC3) instrument with the G141 grism for the exoplanet WASP-6b, the final object in the Sing et al. (2016) study without observations in this wavelength range. Additionally, we present new spectrophotometric observations from 0.4 to 0.8 μm performed from the ground using VLT FORS2. Recent photometric observations of WASP-6b performed from space with the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) are also included in our study. These data sets were analysed in tandem with a re-analysis of the archival STIS and Spitzer data sets on a common Gaussian process (GP) framework (Gibson et al. 2012a). We also perform light-curve corrections to account for the effects of stellar heterogeneity on the perceived transmission spectrum of WASP-6b, the presence of which can act to mimic the signatures of scattering hazes (McCullough et al. 2014; Rackham, Apai & Giampapa 2018; Pinhas et al. 2018; Alam et al. 2018; Rackham, Apai & Giampapa 2019).

Descriptions of our observations and the necessary data reduction are shown in Section 2. All light-curve fitting and analysis is presented in Section 3. An accounting of the effects of stellar heterogeneity is shown in Section 4. The resultant transmission spectra and the conclusions drawn from them using both forward- and retrieval-based models are described in Section 5. Finally, we summarize our results in Section 6.

2 OBSERVATIONS AND DATA REDUCTION

2.1 VLT FORS2

We obtained observations of two primary transits of WASP-6b using the VLT FORS2 GRIS600B (G600B) and GRIS600RI (G600RI) grisms in multi-object spectroscopy mode on 2015 October 3 and November 9, respectively, as part of program 096.C-0765 (PI: Nikolov). These observations utilize a mask with broad slits centred on WASP-6 and a nearby reference star (2MASS J23124095–2243232), all slits had a width of 25 arcmin, the slit lengths used in the G600B and G600RI observations were 31 and 90 arcmin, respectively. On the night of the G600B observations conditions began clear (less than 10 per cent of the sky covered in clouds, transparency variations under 10 per cent) and moved to photometric (no clouds, transparency variations under 2 per cent) approximately half-way through the observations. The exposure time was set at 100 s per exposure for a total of 152 exposures. During this night observations were halted for ∼30 min during transit ingress as the target passed through the zenith and was outside the observable region of the telescope. On the night of the G600RI observations, conditions began clear but moved to photometric for the bulk of the observation and the exposure time was set to 60 s per exposure for a total of 184 exposures. Towards the end of the transit an earthquake caused a guide star loss and as such observations were halted for ∼15 min.

We begin the data reduction by performing bias- and flat-field corrections on the raw data frames, followed by cosmic ray correction using two iterations of the l.a.cosmicalgorithm (van Dokkum 2001). Background flux subtraction for each spectrum was conducted using the median of a box of pixels outside of each spectral trace. Spectra were then extracted using the apall procedure within the iraf package (Tody 1993). Aperture widths for the spectral extraction were varied and values of 14 and 15 pixels were selected as they minimized the dispersion in the out-of-transit flux for the G600B and G600RI white light curves, respectively. We produce a wavelength solution for both observations using the spectra of an emission lamp taken with the calibration mask following each observation. In particular, a low-order Chebyshev polynomial was fit to a multitude of emission lines, the centres of which were determined through individual Gaussian fits. This wavelength solution was then applied to a single data frame to produce a reference spectrum for each observation. Finally, each extracted spectrum was then cross-correlated against its respective reference in order to account for subpixel shifts in the dispersion direction, the maximum resultant shifts were ∼1.2 and ∼0.3 pixels for the G600B and G600RI data sets, respectively. Representative spectra of both WASP-6 and the reference star are shown in Fig. 1 for both the G600B and G600RI observations.

2.2 HST WFC3

A primary transit of WASP-6b was also observed using the HST WFC3 G141 grism on 2017 May 6 as part of General Observer (GO) program 14767 (PI: Sing and López-Morales). All exposures were taken in sequence across five HST orbits, with 13 exposures per orbit, except for the first orbit which only consisted of 10 exposures. Each exposure was performed in forward spatial scanning mode (McCullough & MacKenty 2012), where the telescope slews in the cross-dispersion axis during the exposure, allowing for longer exposure times whilst avoiding saturation on the detector. For the first orbit the exposure times were set to ∼184 s, whilst the remaining orbits had exposure times of ∼138 s. All exposures employed the SPARS25 readout mode and used a scan rate of ∼0.46 pixels s⁻¹.

Reduction of the spectra began with the .ima files output from the CALWF3 pipeline. Each .ima file contains multiple reads for each individual spatial scan, up to the final full scan image. We do not however perform spectral extraction on the final frame of each scan but rather the sum of differenced frames, following Deming et al. (2013). This has the advantage of reducing the impact of cosmic rays and hot pixels, whilst also reducing the overall sky background. For each differenced read, pixels beyond a mask of 35 pixels above and below the centre of the spectral trace were zeroed before extraction of the differenced frame following (Evans et al. 2016). Finally, we then sum all of the differenced frames for each spatial scan to produce a final differenced frame scan.

To perform cosmic ray correction, these frames were stacked into a single cube so that the variation of each pixel could be tracked as a function of time. Each pixel was smoothed temporally with a Gaussian filter and pixel deviations between this and the initial data cube larger than 8σ were flagged as cosmic rays. Static bad pixels were also flagged by searching for deviations greater than 10σ between each individual unsmoothed pixel and the median of a span of 5 pixels in the cross-dispersion direction, centred on the initial pixel. These cosmic rays and static pixels were then replaced by a linear interpolation of the pixel to the point spread function of the same median span. Using a second mask of 50 pixels above and below the centre of the final scans, the 2D spectra were summed along the cross-dispersion axis to produce a 1D spectrum for each scan. This mask width was selected as it provided the minimal white light curve out-of-transit scatter across a range of 30–80 pixels in steps of 5 pixels. The background was subtracted from each spectrum using the median of a box of pixels in a region of the detector unpolluted by the diffuse light from the edges of the spatial scan.

Wavelength solutions were obtained by cross-correlating each individual spectrum with an atlas¹ (Kurucz 1976, 1993) stellar spectrum, with parameters similar to WASP-6 (T_eff = 5500K, log(g)=4.5, and [M/H] = −0.2), convolved with the throughput of the G141 grism. Before cross-correlation, both spectra were smoothed with a Gaussian filter to inhibit the effects of spectral lines and focus the correlation on the steep edges of the G141 throughput. This process revealed shifts in the dispersion direction across the course of observation within ∼0.12 pixels. An example 1D spectrum from the G141 observations is shown in Fig. 1.

2.3 TESS

The TESS is currently performing an all sky search for transiting exoplanets in a single broad-band filter from 0.6 to 1.0 μm (Ricker et al. 2014). Due to the broad 24° × 96° field of view, TESS holds enormous potential not only for discovering new exoplanets, but also observing transits of already known transiting systems. With the public release of the TESS Sector 2 data, 7 clear transits of WASP-6b can be readily identified from 2018 August 23 to September 19.

To obtain the TESS light curve spanning, this time period we initially used the pre-calibrated and extracted light curve held in the lc.fits file. However, on closer inspection we found indications of a non-optimal pipeline correction and as such choose to perform our own correction on the uncorrected light curve in the same file. We follow a Pixel Level Decorrelation (PLD) systematics removal method on the raw data as implemented by the lightkurvepython package (Lightkurve Collaboration et al. 2018). PLD has already been used successfully as a systematics correction technique on both Spitzer (Deming et al. 2015) and K2 data (Luger et al. 2016, 2018) and we refer the reader to these references for further information on the PLD technique itself. Finally, to prepare for the transit light-curve analysis, we extract seven separate portions from the complete light curve, each centred on one of the observed transits. Each individual extracted light curve spans from roughly 5 h pre-transit to 5 h post-transit, in order to facilitate an effective out-of-transit baseline determination.

2.4 Archival data

In order to fully exploit the data that are available to us we opt to perform a re-analysis of the previously reported HST STIS and Spitzer IRAC data (Nikolov et al. 2015). Specifically, there were two spectroscopic transit observations with the STIS G430L grism from 0.33 to 0.57 μm, one spectroscopic transit using the STIS G750L grism from 0.55 to 1.03 μm, and one photometric transit for each of the Spitzer IRAC 3.6 and 4.5 μm bandpasses. Performing such a re-analysis can account for transit depth baseline offsets between these data sets and those in this study by fitting all light curves under a common set of prior system parameters. Furthermore, the implementation of a stellar heterogeneity correction, and its changes to the system parameters (Section 4) necessitates further light-curve fitting. A complete re-analysis ensures that any comparisons between the spot-corrected and uncorrected data sets are not influenced by the differing light-curve fitting methodologies of this study and that of Nikolov et al. (2015).

With respect to the data reduction of the observations themselves, all light curves were extracted following the same methodology outlined in Nikolov et al. (2015). For the STIS data this involves spectral extraction following the APALL procedure in iraf (Tody 1993), and photometry is performed for the Spitzer data through time-variable aperture extraction. For the Spitzer IRAC light curves, there are thousands of independent photometric measurements throughout each observation and to reduce the computational intensity of the light-curve fitting procedure described in Section 3 we bin each light curve into 1000 bins, corresponding to a cadence of ∼15 and ∼16 s for the 3.6 and 4.5 μm bands, respectively.

3 LIGHT-CURVE ANALYSIS

White light curves for the G600B, G600RI, and G141 data sets were produced by summing the flux for each individual spectrum along the dispersion axis from 0.449 to 0.617, 0.529 to 0.833, and 1.0 to 1.8 μm, respectively. Spectrophotometric light curves were produced for the G600B, G600RI, and G141 data sets by summing the flux within 12, 34, and 28 respective bins across the wavelength ranges displayed in Fig. 1.

Below ∼0.45 μm, the G600B flux levels are the lowest of both of the FORS2 data sets and inherently contain a limited amount of information due to the higher photon error. Whilst using a larger bin size could alleviate this, the contribution of differential extinction due to a spectral-type mismatch between the target and reference must also be considered. In the case of our observations such a mismatch is evident in the different spectral profiles of the target and reference star in Fig. 1. The flux of the reference star from 0.40 to 0.45 μm is  50 per cent that of the target, whereas at 0.6 μm this value is  80 per cent. Therefore, the data below 0.45 μm not only contain the lowest flux levels of our FORS2 observations, but their accuracy is impacted the most by the differential extinction. Furthermore, including such a wavelength range would also impart further differential extinction effects on every other spectrophotometric bin in the G600B data set due to the nature of the common-mode, white-light correction performed during the light-curve fitting. In an effort to mitigate the impact of differential extinction on our final transmission spectra we therefore choose to exclude the G600B data below ∼0.45 μm.

In the case of the G600B and G600RI observations, all light curves were also produced for the reference star. Before fitting any of the G600B or G600RI light curves, we first correct for dominant atmospheric effects by dividing the raw flux of the target by that of the corresponding wavelength range reference. The spectrophotometric bins for all observations are displayed in Fig. 1. As the TESS observations are photometric they hold no spectral information and were treated as white light curves in terms of fitting. Finally, we obtain the archival STIS and Spitzer light curves across identical wavelength ranges as described in Nikolov et al. (2015).

During both the G600B and G600RI observations, the target needed to be re-acquired and as such all light curves suffer from incomplete phase coverage, this also results in the separate pieces of each light curve exhibiting differing systematics effects. Throughout our analyses we were unable to accurately and effectively account for these systematic offsets due to the significant, or complete, absence of in transit observations for one piece of each light curve. As such, in the analysis presented here we exclude the pre-ingress data for the G600B observation and the post-egress data for the G600RI observation. The first orbit, and first spectrum of all other orbits, of the G141 observation exhibit much stronger systematics than the other obtained spectra due to charge trapping in the detector (Zhou et al. 2017). We therefore opt to remove these data from our analysis in line with many other studies (e.g. Knutson et al. 2014; Sing et al. 2016; Wakeford et al. 2018; Mikal-Evans et al. 2019) that have been performed since the first spatial scanning WFC3 transit observations were made (Deming et al. 2013).

3.1 White light curves

To perform all light-curve fitting we follow Gibson et al. (2012a), accounting for the transit and instrumental signals simultaneously by treating the data for each light curve as a GP using the python library George (Ambikasaran et al. 2014). GP fitting methodologies have been successfully applied to a range of transit observations (Gibson et al. 2012a,b, 2017; Evans et al. 2013, 2015, 2016, 2017, 2018; Cartier et al. 2017; Kirk et al. 2017, 2018, 2019; Louden et al. 2017; Sedaghati et al. 2017; Mikal-Evans et al. 2019) from both the ground and space thus far and enable the measurement of the systematic signal without assuming any prior knowledge on its functional form. We obtain the best-fitting model to each light curve by marginalizing over the constructed GP likelihoods using Markov Chain Monte Carlo (MCMC) as implemented by the python library emcee (Foreman-Mackey et al. 2013). When executing each MCMC, we first initialized a group of 150 walkers near the maximum likelihood solution, identified using a Nelder–Mead simplex algorithm as implemented by the fmin function within the scipy library. We run a group for 500 samples and then use the best run to initialize a second group of 150 walkers in a narrow space of this solution. This second group was then run for 3000 samples, with the first 500 samples being discarded as burn-in.

We list the individual subtleties for each data set throughout the GP fitting procedure below, however there are some aspects which remained unchanged regardless of the data set. For the GP covariance amplitudes of all data sets, we utilize gamma priors of the form p(A) ∝ e^−100A as in Evans et al. (2018) in order to favour smaller correlation amplitudes and reduce the effects of outliers. Additionally, we follow previous studies and fit for the natural logarithm of the inverse length-scale hyperparameters (e.g. Evans et al. 2017; Gibson et al. 2017; Evans et al. 2018), but limit these quantities with a uniform prior ranging between the cadence of the observation and twice the length of the observation. This prescription encourages the GP to fit the broader systematic variations that occur during the transit, with shorter variations described by white noise and longer variations accounted for by the linear baseline trend. Finally, in all cases the orbital period was held fixed to the value of P  = 3.36100239 d from Nikolov et al. (2015) and the eccentricity was held fixed to the value of e  = 0.041 from Husnoo et al. (2012).

3.1.1 G600B and G600RI

To describe the mean function of the GP we use the model transit light curves of Mandel & Agol (2002) generated using the batmanpython library (Kreidberg 2015) multiplied by a linear airmass baseline trend. We initially tested a time baseline trend however found that this restricted the final GP fitting of shorter frequency variations within the light curves, by utilizing a linear airmass baseline trend the non-linear sloping of the light curves was better matched and the GP had more freedom to fit these shorter frequency variations. Whilst the observed airmass trend can be included in the GP directly as a decorrelation parameter we found this necessitated stricter priors on the length-scale hyperparameters and did not measurably improve the fitting. As such, we opt to include this term through the baseline trend. To construct the covariance matrix of the GP, we use the Matérn ν  = 3/2 kernel, with time as the decorrelation parameter. Other decorrelation parameters were also tested both individually and in combination such as: spectral dispersion drift, cross-dispersion drift, full width at half-maximum, ambient temperature, ambient pressure, and telescope rotation angle. Despite this, no clear correlations were observed and therefore we excluded these parameters from the final analysis.

Unlike the other data sets, for the FORS2 analysis we account for limb-darkening following the two-parameter quadratic law. The treatment is different as these observations were performed from the ground where the Earth’s atmosphere acts as an filter for the incoming light. Crucially, the response of the atmosphere is a function of wavelength and varies as a function of the zenith distance, which varies throughout the observations. Instead of making explicit assumptions about this atmospheric transmission and including it directly in our determination of the pre-computed limb-darkening coefficients we choose to fit for the coefficients themselves. We select the quadratic limb-darkening law in order to improve computational efficiency by reducing the number of fit parameters whilst still providing an accurate description for the true limb darkening of WASP-6 given its temperature (Espinoza & Jordán 2016).

For the G600RI observation, we allow the transit depth R_p/R_*, inclination i, normalized semimajor axis a/R_*, transit central time T₀, linear trend parameters, and quadratic limb-darkening parameters c₁ and c₂ to vary throughout the fit. However, in the case of the G600B observation, we found that the paucity of transit coverage provided imprecise determinations of i and a/R_* and as such perform a simpler fit after retrieving the weighted average best-fitting parameters, see Section 3.1.5.

The presence of high-frequency variations from ∼2 to 3 and ∼0 to 1 h after mid-transit for the G600B and G600RI light curves, respectively, strongly constrain the hyperparameters of the GP fit which leads to overfitting of other variations within the light curve. In order to assess the impact on the fit transit parameters we restricted the priors on these hyperparameters such that the high-frequency variations could no longer bias the GP fitting. Whilst this significantly reduced the perceived overfitting, we find that all fit transit parameters are unaffected by this change and lie within 1σ of the original fit. Therefore, and in addition to the lack of prior knowledge on these hyperparameters, we opt not to perform such a restriction for any of the final white light-curve fits.

3.1.2 STIS and G141

The mean function of the GP is described identically to the G600B and G600RI mean functions, except using a linear time baseline trend. To construct the covariance matrix of the GP, we use the Matérn ν  = 3/2 kernel, with HST orbital phase, dispersion shift and cross dispersion shift identified as the optimal decorrelation parameters. Limb darkening was accounted for through the four-parameter non-linear law. During the fitting we allow the transit depth R_p/R_*, inclination i, normalized semimajor axis a/R_*, transit central time T₀, and linear trend parameters to vary throughout the fit and we fixed all four non-linear limb-darkening values to values calculated from the atlas model described in Section 2, following Sing (2010). Finally, as there are two independent light curves in the STIS G430L observations we performed a joint fit between them, only allowing the transit central time for each light curve to vary independently.

3.1.3 Spitzer

The mean function of the GP is described identically to the G600B and G600RI mean functions, except using a linear time baseline trend. We construct the covariance matrix following Evans et al. (2015). Specifically, we construct a kernel k = k_xy + k_t where k_xy is a squared exponential kernel, with the photometric centroid x and y coordinates as the decorrelation parameters, and k_t is a Matérn s  = 3/2 kernel, with time as the decorrelation parameter. Constructing such a kernel allows us to account for the smooth variations in pixel sensitivities as well as residual correlated noise in the light curve. Limb darkening was accounted for through the four-parameter non-linear law. During the fitting we allow the transit depth R_p/R_*, inclination i, normalized semimajor axis a/R_*, transit central time T₀, and linear trend parameters to vary throughout the fit and we fixed all four non-linear limb-darkening values similarly to the STIS and G141 observations.

3.1.4 TESS

The mean function of the GP is described identically to the G600B and G600RI mean functions, except using a linear time baseline trend. To construct the covariance matrix of the GP we use the Matérn ν  = 3/2 kernel, with time as the decorrelation parameter. Limb darkening was accounted for through the four-parameter non-linear law. During the fitting we allow the transit depth R_p/R_*, inclination i, normalized semimajor axis a/R_*, transit central time T₀, and linear trend parameters to vary throughout the fit and we fixed all four non-linear limb-darkening values similarly to the STIS and G141 observations.

3.1.5 Best-fitting models

In order to obtain the best-fitting model to each data set we determine the weighted average values of the orbital inclination and the normalized semimajor axis (Table 1). Using these values we performed the fit to the G600B data set, where we allowed the transit depth R_p/R_*, transit central time T₀, linear trend parameters, and the quadratic limb-darkening parameters u₁ and u₂ to vary. In addition, we repeat the fit for each other light curve, with the orbital inclination and normalized semimajor axis fixed to the weighted average values and the transit central time to that of its respective original fit. The G600B, G600RI, and G141 light curves, alongside the systematics corrected light curves are displayed in Fig. 2, all TESS light curves are displayed in Fig. 3, all STIS light curves are displayed in Fig. A1, all Spitzer light curves are displayed in Fig. A2, and all relevant MCMC results are displayed in Table 2.

3.2 Spectrophotometric light curves

Prior to the full spectrophotometric fits, we correct all of the spectrophotometric light curves for wavelength independent (common-mode) systematics. In the case of the G600B and G600RI data sets, we follow Nikolov et al. (2016) and determine a common-mode correction by dividing each uncorrected transit white light curve by its final best-fitting transit model. To apply the correction we divide all spectrophotometric light curves by the common-mode calculated from their parent white light curve. For the G141 data set, we correct for common-mode systematics following the shift-and-fit method of Deming et al. (2013). In this case, a reference spectrum was first produced by averaging all of the out-of-transit spectra. Each individual spectrum was then matched against this reference through stretching vertically in flux and shifting horizontally in wavelength following a linear least-squares optimization. We then separate the spectral residuals of the previous fit into 28 wavelength bins spanning 1.13–1.65 μm. Each spectrophotometric residual was then added to a transit model constructed using the best-fitting parameters from the white light-curve fit and limb-darkening calculated for the relative wavelength bin to produce the spectrophotometric light curves. All corrections can be seen under each systematics corrected light curve in Fig. 2.

All spectrophotometric light curves were then fit following the same process as their corresponding white light curves. In each case however, the inclination and normalized semimajor axis were fixed to the weighted average values calculated from the white light-curve fits and the transit central time was fixed to that of each respective white light-curve fit. Additionally, for the G600B and G600RI light curves the quadratic limb-darkening parameter u₂ was fixed to a value calculated from the atlas model described in Section 2 for each individual wavelength bin. The results for all best-fitting transit depths are displayed in Tables B1 and B2 and all spectrophotometric light curves for the G600B, G600RI, G141, and STIS data sets are displayed in Figs 4–6, and A3, respectively.

The initial transmission spectrum of these spectrophotometric light curves revealed an offset in transit depth between the G600B and G600RI data sets. Whilst activity of the host star can lead to such offsets, the stellar variability monitoring performed in Nikolov et al. (2015) shows that potential offsets are of a magnitude ΔR_p/R_* ≃ 0.00022, much too small to account for the observed offset of ΔR_p/R_* ∼ 0.002. Furthermore, the very good agreement of the G600RI data set with the STIS measurements (Section 5.1) of Nikolov et al. (2015) demonstrates that the cause of this offset most likely lies with the G600B data set. Due to the poor phase coverage of the G600B data set, there are almost no observations during ingress, this produces a large uncertainty in the transit central time and subsequently the absolute transit depth, which may be responsible for the offset we see. Therefore, to account for this offset we apply a vertical shift to the G600B data set by performing a weighted least-squares minimization on the difference between the spectrophotometric bins in the overlapping region between the G600B and G600RI data sets, leaving the relative vertical shift of the G600B data set as a free parameter in the minimization. This results in a shift of ΔR_p/R_* = 0.00248, equivalent to ∼1.5σ of the error on the transit depth of the G600B white light curve. A full transmission spectrum with this offset included is shown in Fig. 7.

4 CORRECTING FOR STELLAR HETEROGENEITY

Stellar activity leads to the presence of heterogeneities on stellar surfaces through the magnetically driven formation of cooler regions known as star-spots and hotter regions known as faculae. The presence of spots (or faculae) on the surface of a star results in a wavelength dependent variation in the stellar baseline flux due to the respective differences in the emission profiles of the relatively cool spot (or relatively hot faculae) and the stellar surface itself. As the stellar baseline flux is crucial in determining transit depth, the presence of an unocculted star-spot during a transit observation will necessarily produce a wavelength-dependent variation in the measured transit depth (Rackham et al. 2018, 2019). If significant enough, this variation can produce an artificial slope in the optical region of the final measured transmission spectrum, potentially mimicking the effects of haze in the atmosphere (Pont et al. 2008; Sing et al. 2011a; McCullough et al. 2014; Alam et al. 2018; Pinhas et al. 2018). These wavelength-dependent variations can also impact individual spectral features due to the differential emission of specific stellar lines. Previous studies have displayed small decreases in the amplitude of Na i absorption following a stellar heterogeneity correction (e.g. Sing et al. 2011a; Alam et al. 2018), however this effect is typically secondary to the artificially induced optical slope.

To estimate the impact surface stellar heterogeneities may have on our observations we obtained a proxy of the magnetic activity level of WASP-6 using a measurement of |$\log (R^\prime _{HK})$|⁠. This value has been previously quoted without uncertainties as −4.741 in Sing et al. (2016), however analysis of the emission cores of the Ca ii H and K lines in the HARPS spectra of Gillon et al. (2009) results in a direct measurement of −4.511 ± 0.037, indicating that WASP-6 is a moderately active star compared to the broader population of cool stars (Boro Saikia et al. 2018). We therefore endeavour to account for the effects of unocculted star-spots following the methodology of Alam et al. (2018).

4.1 Photometric monitoring of WASP-6

We estimate the long baseline variability of WASP-6 by considering all 18 317 images from the TESS observations previously described in Section 2.3 in addition to 435 R-band images from the Tennessee State University 14-inch Celestron Automated Imaging Telescope (AIT) taken from 2011 September to 2019 January 2019 (Fig. 8). Initially, we also incorporated 738 V-band images taken from 2013 November to 2018 July as part of The Ohio State University All-Sky Automated Survey for Supernovae (ASAS-SN) (Shappee et al. 2014; Jayasinghe et al. 2018) into our photometric monitoring data set as in Alam et al. (2018). However, on comparing the contemporaneous ASAS-SN and TESS data we find a ∼4 times larger photometric scatter in the ASAS-SN data set compared to the more precise TESS sample and, as such, exclude it from our analysis in order to avoid influencing the variability amplitude estimation with such a noise-dominated data set.

4.2 The stellar rotation period

In order to perform an accurate fit to the photometric monitoring data, it is necessary to have a measurement of the stellar rotation period. However, a range of rotation periods have been reported for WASP-6. In particular, Jordán et al. (2013) find a period of 16 ± 3 d based on the vsinI  = 2.4 ± 0.5 km |$\rm {s}^{-1}$| measurement from Doyle et al. (2013), Nikolov et al. (2015) determine a period of 23.6 ± 0.5 d from a portion of their AIT photometric monitoring, and by tracking transit star-spot crossings Tregloan-Reed et al. (2015) find a period of 23.80 ± 0.15 d, assuming the star had rotated only once between successive observed crossings.

We also perform a measurement of this rotation period through virtue of the very high cadence TESS observations. Even from an initial inspection of the light curve shown in Fig. 8, a clear sinusoidal variation can be seen. In order to determine that this variation is not due to an instrumental effect we inspect the light curves and background flux of the four closest neighbouring stars to WASP-6 with TESS light-curve observations. We find that none of the stars exhibit the same sinusoidal variation as WASP-6, and they all exhibit similar variations in their background flux. To determine the rotation period itself, we perform a least-squares minimization using a simplistic sinusoidal model on the data with all transit events removed. This resulted in an inferred period of 12.18 ± 0.05 d.

Even though this method of model fitting is quite rudimentary, the determined period is clearly in contradiction to current estimates of the stellar rotation period. This contradiction suggests that the variability observed is likely not that of a single spot feature rotating with a period equal to that of the stellar rotation period. Alternatively, the perceived TESS period can be explained by the spot coverage during the TESS epoch being concentrated on opposite hemispheres of the star, rather than one single hemisphere. During a period of AIT photometry performed shortly after the TESS observations from 2018 September to 2019 January, we find a standard deviation of 3.8 mmag, in contrast to previous seasons where this reached up to 8.1 mmag. This reduced variability is further justification of the measured TESS period being a result of hemispherically varying star-spot coverage and not intrinsic to the TESS instrument itself. Further high-quality photometric monitoring will likely be necessary to fully resolve the discrepancy between these observations. For subsequent analysis however, we adopt the stellar rotation period of 23.6 ± 0.5 d from Nikolov et al. (2015) as this estimate was made over much longer timescales compared to the estimates of Jordán et al. (2013) and Tregloan-Reed et al. (2015).

4.3 Modelling and correction of unocculted star-spots

The variability of WASP-6 was modelled following the methodology of Alam et al. (2018). We perform a GP regression model fit to the photometric monitoring data constructed with a three-component kernel which models: the quasi-periodicity of the data, irregularities in the amplitude, and stellar noise. A gradient-based optimization routine was used to locate the best-fitting hyperparameters and a uniform prior was placed on the stellar rotation period, centred on the value of 23.6 ± 0.5 d from Nikolov et al. (2015) with a width three times that of the standard deviation. The TESS bandpass ranges from 0.6 to 1.0 μm and is less susceptible to active photometric variations compared to the AIT R-band observations. This should not affect the wavelength dependence of our determined spot correction however, as the estimated variability amplitude is ultimately used as a reference to normalize the true model wavelength-dependent correction factor (equation 1). Despite this, the discrepancy of the measured TESS period from the measured period in other studies (Jordán et al. 2013; Nikolov et al. 2015; Tregloan-Reed et al. 2015), and the reduced variation in a subset of AIT data described in Section 4.2, does indicate that the variability of the star as a whole was also lower during this epoch. Because the variability amplitude is crucial in determining the spot correction, we opt to perform separate fits to the TESS and AIT data sets. To avoid influencing the GP fitting with the lower variance AIT data, we exclude 41 measurements obtained shortly after the TESS epoch which correspond to the subset described in Section 4.2. Due to the large size of the TESS data set (∼18 000 data points), we bin the data down by a factor of 10 in order to make the GP fitting computationally tractable.

Whilst the TESS data is well sampled and more precise than the AIT data, we may be perceiving a lower level of variability due to the TESS bandpass or the lower intrinsic variability of WASP-6 during the TESS epoch (Section 4.2). Comparatively, the AIT data have a much broader temporal coverage and could therefore be more indicative of the longer term variability of WASP-6, though as there are no contemporaneous measurements with the TESS data set their accuracy is not guaranteed. The TESS and AIT model fits therefore provide respectively more conservative or realistic estimates of the true stellar variability. All such fits to the photometric monitoring data are displayed in Fig. 8.

We are then able to correct for the unocculted spots in the transit light curves following Huitson et al. (2013). Under the assumption that there is always some level of spot coverage on the stellar surface, the maximum observed stellar flux does not correspond to the flux emitted by an entirely unspotted surface. Using the amplitude of the GP fit to both the TESS and AIT photometric monitoring data, we determine different estimates for the unspotted stellar flux F^′  = max(F) + kσ, where F is the observed photometric monitoring data, σ is the dispersion of these photometric measurements, and k is a value fixed to unity. Whilst an accurate value of k can be difficult to determine a k  = 1 has been shown to be suitable for active stars (Aigrain, Pont & Zucker 2012). Furthermore, varying the chosen value of k does not significantly influence the wavelength dependence of the correction and mainly influences the offset of the transmission spectrum baseline (Alam et al. 2018). For each estimate, the fractional dimming due to stellar spots was then calculated as |$f_{\rm norm} = \overline{F/F^{\prime }}$|⁠, giving the amplitude of the spot correction at the variability monitoring wavelength as Δf₀  = 1 – f_norm.

5 DISCUSSION

The observed transmission spectrum of WASP-6b reveals a variety of spectroscopic features present both in the uncorrected and spot-corrected analyses (Fig. 10). In particular, the broad absorption feature at 1.4 μm indicates the presence of H₂O in the atmosphere. Additionally, narrow-band absorption features at 0.589 and 0.767 μm due to Na i and K i are also evident in the optical. Finally, a distinct increase in transit depth across optical wavelengths is seen, indicative of a scattering haze and in agreement with Nikolov et al. (2015). The primary difference between the uncorrected and spot-corrected data sets is the presence of a vertical offset across the full wavelength range. This offset is not wavelength independent however and the spot correction has acted to slightly reduce the gradient across the optical slope. This wavelength dependence is clearly identified by the difference in transit depth between the uncorrected and AIT-corrected data sets at the shortest wavelength bin compared to that of the longest wavelength (Fig. 10).

5.1 Archival data comparisons

The transmission spectrum of WASP-6b had already been measured using the available HST STIS and Spitzer IRAC data sets (Nikolov et al. 2015). In order to compare our independent reduction against these results we overplot both the uncorrected transit depths from this study, with those from this prior published study (Fig. 11). The different reductions agree quite well, with all measurements within 1σ of one another. A minor discrepancy in transit depth is seen for the longest wavelength STIS bins and the Spitzer photometry. These discrepancies are likely due to the slightly different measured system parameters which were held fixed during the independent fittings in addition to slight differences in the adopted stellar limb-darkening parameters. The error bars for the reduction performed in this study are larger than those of those from the original reduction, primarily due to the difference between the model marginalization and GP approaches towards light-curve fitting.

As the STIS and VLT FORS2 data sets have a broad overlapping wavelength range we reproduce the VLT FORS2 transmission spectrum using an identical wavelength binning as the HST STIS measurements to facilitate a comparison between the results (Fig. 11). It is evident from this comparison that whilst our results agree very well at the shortest and longest wavelengths, there is a small disparity in the measurements centred around the Na i absorption line. We calculate a weighted average transit depth across five wavelengths bins centred on the Na i absorption line for the G600RI and the STIS data sets, resulting in R_p/R_*’s of 0.14628 ± 0.00031 and 0.14520 ± 0.00043, respectively. We exclude the G600B data set from the calculation to avoid any bias due to the applied vertical shift as described in Section 3.

As the offset reduces proportionally with separation from the Na i line centre, this signal could be indicative of an observation of the pressure-broadened wings from the full Na i feature in the FORS2 data sets. Such wings have recently been definitively observed in the atmosphere of the hot Jupiter WASP-96b (Nikolov et al. 2018). Given these wings are not present in the STIS data set, this could suggest we are observing variability in the atmosphere of WASP-6b. However, this offset being of an instrumental or systematic origin cannot be excluded, particularly as the FORS2 observations are taken from the ground where systematic variations are not as well understood and harder to model. The possibility that this discrepancy has been caused by the STIS observations in particular also cannot be excluded as there exists robust evidence that systematics in STIS observations resulted in a spurious detection of K in WASP-31b (Gibson et al. 2017, 2019). The true cause of the discrepancy, be it physical or systematic, can not be determined with these data and additional observations at higher signal-to-noise ratio and over long time-scales will be required to investigate this further.

5.2 Goyal forward models

In order to explore the bulk properties of WASP-6b, we fit the observed transmission spectrum to a grid of forward models (Goyal et al. 2018, 2019). These models are generated using the 1D radiative-convective equilibrium code atmo. Initially we opted to use the more recent generic model grid (Goyal et al. 2019) in our analysis as it allowed for a broader coverage of the parameter space than the WASP-6b specific grid from Goyal et al. (2018). However, as subsolar metallicity forward models have yet to be implemented into the generic grid our ability to accurately fit the observed data was ultimately restricted. As such, we used the WASP-6b specific grid (Goyal et al. 2018) in order cover the subsolar metallicity range of parameter space.

With the arrival of the Gaia Data Release 2 (Gaia Collaboration et al. 2018), the distance to WASP-6 has been more accurately determined as |$d = 197.1^{+0.4}_{-1.6}$| pc (Bailer-Jones et al. 2018), significantly different to the prior measurement of 307 pc. This re-estimation has significant effects on the inferred stellar radius of WASP-6 which in turn affects the estimation of planetary radius from the observed transit depths. A mismeasurement of the planetary radius naturally leads to a mismeasurement of the planetary gravity, a currently fixed parameter for the planet-specific forward model grid of Goyal et al. (2018). Following the methodology of Morrell & Naylor (2019), we performed spectral energy distribution (SED) fitting on WASP-6 using near-ultraviolet (NUV), optical, and near-infrared broad-band photometry. The fitted integrated flux allows us to measure its luminosity, and the shape of the SED determines its so-called T_SED (see Morrell & Naylor 2019, for details). By combining this with the revised distance measurement, we obtained an updated estimate of the radius of WASP-6, and subsequently the radius of WASP-6b. This radius results in a new value for the planetary gravity of |$g = 10.55^{+0.19}_{-0.39}$| ms⁻², notably different from the previous estimate of g = 8.71 ± 0.55 ms⁻² (Gillon et al. 2009). Changes in gravity can have significant effects on the computed forward models (Goyal et al. 2018, 2019) and therefore to fit our observed data we use a more updated forward model grid for WASP-6, identical to the original shown in Goyal et al. (2018) except recomputed for a value of g = 10.5.

The model grid used consists of 3920 different transmission spectra varying in temperature, metallicity, C/O ratio, scattering haze, and uniform cloud. The scattering haze is implemented through the use of a haze enhancement factor α_haze which simulates an increase in the total scattering of small aerosol particles in the atmosphere. Similarly, the uniform cloud is implemented through a variable cloudiness factor α_cloud, which produces the effects of a cloud deck through a modification to the wavelength-dependent scattering using the strength of grey scattering due to H₂ at 350 nm. Irrespective of the true cloud composition, implementing a grey cloud is appropriate for our observations as at the observed wavelengths Mie scattering predicts essentially grey scattering profiles (Wakeford & Sing 2015). Further details on the grid parameters, including their ranges and implementations, can be found in Goyal et al. (2018).

Each model spectrum was fit in turn by producing a binned version of the spectrum which matches the selected spectrophotometric bands from the data reduction and then averaged to produce a single value of transit depth in each bin. A χ² measurement between the observed and model data was then computed following a least-squares minimization scheme with a varying wavelength-independent vertical offset. These fits were performed for both the uncorrected and both spot-corrected transmission spectra and the best-fitting models for each are presented in Fig. 10.

For the uncorrected and the TESS-corrected transmission spectra, we find a best-fitting model of |$T = 1334\,$| K, subsolar metallcity [M/H] = −1.0, slightly supersolar C/O ratio of [C/O]  = 0.70, moderate hazes α_haze = 10, and no evidence of clouds α_cloud = 0 corresponding to a |$\chi ^2_\nu$|  = 1.10 and 0.98, respectively. For the AIT-corrected transmission spectrum however, we find a best-fitting model of |$T = 1334\,$| K, subsolar metallcity [M/H] = −1.0, solar C/O ratio of [C/O]  = 0.56, moderate hazes α_haze = 10, and no evidence of clouds α_cloud = 0 corresponding to a |$\chi ^2_\nu$|  = 0.99. To explore the discrepancies and commonalties between the grid fits to the uncorrected and corrected data sets we produce χ² contour maps (Madhusudhan & Seager 2009) as shown in Fig. 12. We begin by constructing 2D grids of every possible pair of model parameters. In each separate grid, and at every individual grid point dictated by the resolution of the model parameters, we vary all the remaining model parameters in turn and determine the model with the smallest χ². Across these new χ² spaces we determine contours which correspond to models in the parameter space which are Nσ from the overall best-fitting model following Goyal et al. (2018).

The primary differences between the data sets are the existence of subsets of model fits more favoured by the lowest metallicities and the highest haze enhancement factors for only the uncorrected data set. These subsets are present because the wavelength dependence of stellar heterogeneity acts to increase the gradient of the optical slope in the observed data, an effect that is somewhat degenerate with lower metallicity and hazy atmospheres (Goyal et al. 2018). Whilst both the lowest metallicities and highest haze enhancements factors are not as favoured in tandem, they both correspond to model fits favouring a lower level of C/O ratio. This is because both low metallicity and high haze enhancement factor act to suppress the H₂O absorption features beyond the constraints set by the G141 data set and as such the C/O ratio must be reduced in order to re-inflate the H₂O features to match the observations. In summary, the χ² contour map for even the conservative TESS-corrected data set indicates that these highest haze enhancement factors, lowest metallicities, and lowest C/O ratios are likely effects of stellar heterogeneity on the transmission spectrum of WASP-6b and not truly symptomatic of its atmosphere. However, a moderate haze enhancement of at least α_haze = 10 is strongly constrained, and a preference towards subsolar metallicities is still evident, independent of the addition of a spot correction.

Whether or not a spot correction is used, temperatures of |$1334 \,$| K are primarily preferred for each grid fit. Comparatively, the measured dayside temperatures for WASP-6b are |$1235^{+70}_{-77}$| and |$1118^{+68}_{-74}$| from the 3.6 and 4.5 μm Spitzer IRAC channels, respectively (Kammer et al. 2015). As these values are within ∼1σ they do not suggest a disagreement, however, it is worthwhile assessing the source of the slight preference of the grid model fits towards limb temperatures higher than that measured from the dayside. As the model grid varies in temperature steps of |$150 \,$| K, the model cannot settle on a precise temperature estimate and is therefore likely to be somewhat discrepant from the true value. However, there are models at a temperature |$T = 1184\,$| K which should in theory match the true temperature of WASP-6b’s limb more accurately. Looking to Fig. 12, the preferred temperature is strongly constrained below the |$1484 \,$| K grid models, as at approximately this temperature absorption features due to TiO and VO start to become significant in the optical (Fortney et al. 2008) and are strongly disfavoured by the observed FORS2 and STIS data sets. As temperature acts to increase the gradient of the optical slope (Goyal et al. 2018), it is also degenerate with the effects of stellar heterogeneity. Therefore the models at |$1334 \,$| K are the most favoured as it is the highest temperature, and thus steepest slope, that the model grid can produce without generating conflicting TiO and VO features. Fig. 12 demonstrates this as the model preferences for the highest temperatures are slightly reduced upon application of the spot corrections, with the most significant difference being for the AIT-corrected data set. As the best-fitting temperature for the AIT correction is still beyond what we would expect given the dayside temperatures already reported it could even suggest that the spot correction used has been underestimated. However, a subset of |$1184\,$| K models are comfortably within the 2σ region for every data set and therefore conclusively determining the true effect of stellar heterogeneity on the best-fitting model temperature will require further investigation with observations at a higher signal-to-noise ratio.

To determine the significance of the perceived detections of the Na i and K i features we begin by performing a quadratic interpolation of the baseline of the best-fitting model to each data set from 0.4 to 0.9 μm using regions of the optical slope with no clear absorption features as anchors for the interpolation. The interpolation then served as a comparison against the weighted mean value of the G600B, G600RI, STIS 430, and STIS 750 data contained with the Na i and K i lines. Detection significances are summarized in Table 3, these values indicate at least a 3σ detection of the Na i and K i narrow-line signatures in the atmosphere WASP-6b, irrespective of an applied spot correction.

5.3 ATMO retrieval modelling

The previously available transmission spectra of WASP-6b has been the subject of multiple retrieval-based model analyses thus far. First by Barstow et al. (2017) who utilize the nemesis retrieval code to demonstrate that the atmosphere of WASP-6b is best described by Rayleigh scattering clouds at high altitudes. In addition, Pinhas et al. (2018) perform a retrieval using the aura code, demonstrating that the atmosphere of WASP-6b is best described as a combination of the effects of stellar heterogeneity and atmospheric hazes. However, in an effort to fit the widely disparate STIS and Spitzer points this retrieval predicts a very low H₂O abundance, a claim that has not been possible to verify or refute until the recent acquisition of HST WFC3 data from this study.

Due to the wealth of new data available with the addition of the FORS2, WFC3, and TESS observations, we perform our own atmospheric retrieval on the uncorrected and spot-corrected data sets using the atmo Retrieval Code (arc) which has already been used for a variety of transmission spectra to date (Wakeford et al. 2017a, 2018; Nikolov et al. 2018; Spake et al. 2018; Evans et al. 2018). For the retrieval model, the relative elemental abundances for each model were calculated in equilibrium. For each model, equilibrium chemistry was calculated on the fly, using input elemental abundances fit by assuming solar values and we allowed for non-solar scaled elemental compositions by fitting the carbon, oxygen, sodium, and potassium elemental abundances ([C/H], [O/H], [Na/H], and [K/H]), which can potentially all be constrained by the transmission spectrum. We fit all remaining species by varying a single quantity for the trace metallicity, [M_trace/H]. Throughout this study, all abundances are quoted as [X/H] which is logarithmic relative to the Sun, with all solar abundances taken from Asplund et al. (2009). The resulting chemical network consisted of 175 neutral gas phase species, 93 condensates, and the ionized species e⁻, H⁺, H⁻, He⁺, Na⁺, K⁺, C⁺, Ca⁺, and Si⁺. By varying both C and O separately, we mitigate several important modelling deficiencies and assumptions compared to varying the C/O ratio as a single parameter (Drummond et al. 2019). For the spectral synthesis, we included the spectroscopically active molecules of H₂, He, H₂O, CO₂, CO, CH₄, NH₃, Na, K, Li, TiO, VO, FeH, and Fe. The temperature was assumed to be isothermal, fit with one parameter, and we also included a uniform haze fit with the enhancement factor. A differential-evolution MCMC was used to infer the posterior probability distribution which was then marginalized (Eastman, Gaudi & Agol 2013), we ran 12 chains each for 30 000 steps, discarding the burn-in before combining them into a single chain. Uniform priors were adopted, with the log₁₀ abundances allowed to vary between −12 and −1.3.

The resulting best-fitting retrieval models for the uncorrected, TESS-corrected, and AIT-corrected data sets all provide good fits to the data, with χ²  = 75, 71, and 73, respectively, for 86 degrees of freedom. We show a visual representation of the retrieval for the AIT-corrected data set in Fig. 13 and the mean values for each individual retrieval are shown in Table 4. To facilitate comparisons between the uncorrected and corrected data sets, we plot the retrieval posteriors for each together in Fig. 14. As with the forward model grid fits shown in Section 5.2, there are clear differences between the uncorrected and spot-corrected data sets, particularly for the temperature, radius, and haze opacity. The difference in radius is a natural result of performing the spot correction, as this results is a wavelength-dependent shift in the transmission baseline to lower transit depths. Given the square root of the transit depth δ  = R_p/R_*, and that the stellar radius is fixed during the retrieval, any decrease in the transit depth will subsequently produce a decrease in the estimated planetary radius. In a similar fashion to the forward model grid fits, the highest temperatures and highest levels of haze opacity are favoured by the uncorrected data set, the cause of which being the degeneracy between these properties and the effects of stellar heterogeneity on the uncorrected transmission spectrum. Upon performing a spot correction, the best-fitting temperature and haze opacity falls as the gradient of the optical slope has been reduced. However at least a moderate amount of haze is still required irrespective of spot correction.

Due to the freedom of the retrieval analyses, we were also able to investigate the specific elemental abundances inferred from the measured transmission spectra. First, as the C, O, Na, and K abundances were fit independently throughout the retrieval analysis the measured metallicity only encompasses the other elemental constituents of the atmosphere. The subsolar metallicity measured across all retrieval analyses therefore show that no other substantial absorber is required to fit the measured transmission spectra. The Spitzer data points are the only observations sensitive to carbon-bearing species in the atmosphere such as CH₄, CO, and CO₂, however, given their non-negligible uncertainties and minimal relative offset the retrieved carbon abundance is largely unconstrained and merely represents an upper limit. This is true across all data sets as the addition of a stellar heterogeneity correction has a marginal effect towards the infrared. We constrain the carbon abundance to subsolar at 3σ for the uncorrected and AIT data sets, and at 2σ for the TESS data set. Our limit on the carbon abundance suggests that H₂O is the primary oxygen-bearing species, and from the observed feature we constrain the oxygen abundance to a subsolar value, irrespective of a spot correction. For the best-fitting retrieval model to the AIT-corrected data set our oxygen abundance corresponds to a water abundance of log(H₂O)  = −4.87. Given the lack of WFC3 data available to previous studies of WASP-6b, this water abundance is the first to be informed by an observed water absorption feature in transmission. Furthermore, given the extensive optical data from FORS2 and STIS, this result is robust to previously observed degeneracies of water abundance and reference pressure (Griffith 2014; Pinhas et al. 2018). Contrasting to oxygen, the Na and K abundances are relaxed to lower values following the application of a spot correction as the lone Na i and K i absorption features lie in the optical region where stellar heterogeneity has a significant effect on the observed slope. Upon a reduction in the slope opacity, these abundances must necessarily drop to fit the observed data. Specifically for the AIT correction, we see variations in sodium of supersolar, [Na/H]  = |$1.33^{+0.42}_{-0.67}$|⁠, to solar/supersolar, [Na/H]  = |$0.83^{+0.67}_{-0.80}$|⁠, and potassium of solar/supersolar, [K/H]  = |$0.22^{+0.65}_{-0.74}$|⁠, to subsolar/solar, [K/H]  = |$-0.12^{+0.71}_{-0.74}$|⁠. Given the measurement precision, we cannot explicitly quantify the impact of the correction as both the [Na/H] and [K/H] abundances lie within 1σ of their inferred uncorrected abundances. Despite this, the broader shifts of their full retrieved distributions (Fig. 14) indicate that neglecting to account for the affects of stellar heterogeneity in future, higher precision, observations may lead to strictly incorrect determinations of their abundances.

As the metallicity we retrieve excludes C, O, Na, and K, we cannot perform a comparison to the [M/H] distributions obtained as part of the forward model analysis in Section 5.2. However, comparing the retrieved [O/H] to the forward model [M/H] we see similar distributions indicating a subsolar metallicity. Additionally, whilst the slightly supersolar abundances of [Na/H] and [K/H] do not completely agree with the subsolar [M/H], the large uncertainties of these distributions indicate that an overall subsolar metallicity cannot be ruled out.

5.4 WASP-6b in context

Our determined, spot-corrected, oxygen abundance of [O/H]  = |$-0.84^{+0.40}_{-0.39}$| and sodium abundance of [Na/H]  = |$0.83^{+0.67}_{-0.80}$| are slightly disparate to the determined subsolar metallicity of the host star of [Fe/H]  = −0.15 ± 0.09 (Doyle et al. 2013), whilst the potassium abundance is in good agreement at [K/H]  = |$-0.12^{+0.71}_{-0.74}$|⁠. Variations in these elemental abundances relative to the host star could be indicative of formation history (e.g. Öberg, Murray-Clay & Bergin 2011), however in the case of WASP-6b, the current uncertainties are not sufficiently constrained to make such determinations, with all values lying within 2σ of the host star metallicity. Further observations of the atmosphere of WASP-6b will be necessary to provide more detailed constraints on these elemental abundances. In particular, due to the presence of carbon-bearing molecular features beyond 2 μm such as CO, CO₂, and CH₄, spectroscopic observations with the upcoming James Webb Space Telescope (JWST) will provide stronger constraints on the carbon abundance, of which this study could only provide an upper limit. This in turn will enable robust constraints on the C/O ratio and progress our understanding of the formation history of WASP-6b.

Irrespective of the application of a stellar heterogeneity correction, both the forward and retrieval models require some level of haze opacity enhancement in order to describe the steep optical slope of the transmission spectrum. In the context of hot Jupiter atmospheres, this haze is often thought of as either photochemically produced, or condensate dust, scattering species within the atmosphere (Marley et al. 2013). In the case of the condensate species it is thought that the lofting of particles from deeper atmospheric cloud decks can serve to populate the upper atmosphere and lead to the observed scattering we see (e.g. Parmentier, Showman & Lian 2013). Despite this, the most recent simulations of condensate particle formation in the atmosphere of the hot Jupiter HD 189733b (Bouchy et al. 2005) fail to fully reproduce its observed scattering slope (Lee et al. 2017; Powell et al. 2018). At the temperature of WASP-6b, generation of hydrocarbons through photochemistry was initially thought to be inhibited (Liang et al. 2004) and whilst sulfur photochemistry may play a role (Zahnle et al. 2009), it primarily induces a scattering slope below 0.45 μm, whereas the observed slope of WASP-6b extends further into the optical. However, recent laboratory experiments have shown that hydrocarbons may form not just in cool exoplanet atmospheres (Hörst et al. 2018; He et al. 2018), but also in hot atmospheres beyond 1000 K with a sufficiently high [C/O]  = 1 (Fleury et al. 2019), a possibility our observations cannot definitively rule out. Additionally, the effects of wind-driven chemistry act to homogenize the atmospheres of tidally locked hot Jupiters such as WASP-6b and can lead to significant increases in the abundance of CH₄ compared to standard equilibrium models (Drummond et al. 2018a,b). Given photolysis of CH₄ can drive the formation of haze precursors (Lavvas, Coustenis & Vardavas 2008), this increase in abundance may naturally lead to their more efficient production. Furthermore, of the well-characterized hot Jupiter atmospheres, WASP-6b and HD 189733b present an interesting comparison as they have similar temperatures, both orbit active stars (⁠|$\log (R^\prime _{HK})$|  = −4.511 and −4.501, respectively), and both exhibit strong haze scattering slopes across the optical (Sing et al. 2016). Recent simulations of HD 189733b by Lavvas & Koskinen (2017) have shown that the formation of photochemical haze ‘soots’ higher in the atmosphere are not excluded and can match its observed transmission spectrum. Moreover, the increased UV flux that these two planets are subject to due to their large host star activity levels is likely acting to enhance the rate of photochemical haze production in their atmospheres (Kawashima & Ikoma 2019). Possible evidence to this conclusion is seen in the potential trend towards stronger scattering haze signatures with reducing |$\log (R^\prime _{HK})$| (increasing activity) observed in the hot Jupiter population study of Sing et al. (2016). An exact determination of whether the haze produced in the atmosphere of WASP-6b is of photochemical origin, condensate dust origin, or a combination of the two, was not possible as part of this study due to their similar opacities at the wavelengths of these observations (e.g. Nikolov et al. 2015). In future analyses however, the relative contributions of both photochemical and condensate haze components should be considered in order to describe this observed scattering.

Amongst the population of spectroscopically studied exoplanets, the atmosphere of WASP-6b is one of the haziest. Previous studies of its atmosphere predicted a small (Nikolov et al. 2015; Sing et al. 2016) amplitude H₂O feature at 1.4 μm, however the feature observed as part of this study is slightly larger than anticipated. This increase is likely due to the seemingly small Spitzer transit depths biasing the model estimates prior to the acquisition of the FORS2 and WFC3 data sets. To quantify the size of the H₂O feature relative to an assumed clear atmosphere for WASP-6b, we determine the scaled amplitude of the water feature following Wakeford et al. (2019). Specifically, we begin by taking a clear atmosphere forward model from the grid used throughout this paper (Goyal et al. 2018) with: the equilibrium temperature of WASP-6b, solar metallicity, solar C/O ratio, and no haze or cloud opacity components. We then scale this model to fit the data using a model defined as S1 = (S0 × p₀) + p₁, where S0 is the clear atmosphere model, p₀ is the model amplitude scale factor and p₁ is a baseline offset. For the AIT-corrected data set, we determine p₀  = 64 ± 12 per cent, in contrast to the median amplitude across the observed population of p₀  = 33 ± 24 per cent (Wakeford et al. 2019). These new observations indicate that despite the presence of haze, WASP-6b remains a favourable target for atmospheric characterization, particularly with JWST. This potential for JWST to characterize hazy hot Jupiters such as WASP-6b is in contrast to those who exhibit flat, cloudy spectra such as WASP-31b (Gibson et al. 2017) and WASP-101b (Wakeford et al. 2017b).

6 CONCLUSIONS

We present the most complete optical to infrared transmission spectrum of the hot Jupiter WASP-6b to date utilizing new observations performed with HST WFC3, VLT FORS2, and TESS in addition to re-analysed existing HST STIS and Spitzer IRAC data. The impact of host star heterogeneity on the transmission spectrum was investigated and we correct the observed light curves to account for these effects under different assumptions for the level of stellar activity. All reduced transmission spectra then undergo a retrieval analysis fitting in addition to being fit to a grid of forward atmospheric models.

Across all data sets we find clear evidence for Na i, K i, and H₂O within the atmosphere of WASP-6b in addition to a steep increase in transit depth towards the optical. After applying both forward model and retrieval analyses, we find that at least a moderate haze enhancement is required to describe the optical slope, however when neglecting even a conservative stellar heterogeneity correction, higher and potentially erroneous haze enhancement factors are more preferred. An analogous effect is also seen in the estimated temperature, where higher and potentially unphysical temperatures are preferred when there is no stellar heterogeneity correction. Both of these effects likely stem from the degeneracy of these properties and the impact of stellar heterogeneity towards increasing the optical slope of the transmission spectrum.

Whilst the precision of current observations is not sufficient to definitively estimate the impact of stellar heterogeneity on the transmission spectrum of WASP-6b, the parameter differences observed upon the application of a stellar heterogeneity correction indicate that its effect should not be neglected for future observations of exoplanetary atmospheres around moderately active stars. Despite the presence of haze in its atmosphere, WASP-6b remains a favourable target for further characterization. Contemporaneous and broader wavelength measurements of its transmission spectrum with missions such as JWST will enable a more detailed characterization of its atmosphere in addition to the precisely determining the effects stellar heterogeneity has on its appearance.

ACKNOWLEDGEMENTS

This work is based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere under European Southern Observatory programme 196.C-0765 in addition to observations associated with program GO-14767 made with the NASA/ESA HST that were obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc, under NASA contract NAS 5-2655. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. Based on observations made with the NASA Galaxy Evolution Explorer, GALEX is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034. ALC and SM are funded by a UK Science and Technology Facilities Council (STFC) studentship. MKA acknowledges support by the National Science Foundation through a Graduate Research Fellowship. HRW acknowledges support from the Giacconi Prize Fellowship at STScI, operated by AURA. This research has made use of NASAs Astrophysics Data System and the python modules numpy, matplotlib, and scipy.

Footnotes

REFERENCES

APPENDIX A: ARCHIVAL LIGHT-CURVE FITS

APPENDIX B: SPECTROPHOTOMETRIC LIGHT-CURVE FITS