The Very Low Albedo of WASP-12b from Spectral Eclipse Observations with Hubble

Thermal measurements of hot Jupiters suggest that these gas giant exoplanets often have moderate Bond albedos (${A}_{B}\approx 0.4$, the fraction of incident energy reflected to space; Schwartz et al. 2017). However, many previous searches for reflected light from hot Jupiters have found little-to-none at optical wavelengths where the host star emits most of its energy (geometric albedo ${A}_{g}\lt 0.1;$ e.g., Rowe et al. 2008; Kipping & Spiegel 2011; Heng & Demory 2013; Dai et al. 2017). It is unclear what is causing this apparent contradiction between constraints from thermal emission and optical reflection. Previous Hubble Space Telescope (HST) eclipse observations of HD 189733b with the Space Telescope Imaging Spectrograph (STIS) showed an increase in reflectivity toward bluer wavelengths which may, at least in part, explain the discrepancies between these two techniques (Evans et al. 2013).

A direct way to probe the back scattering efficiency of a hot Jupiter's atmosphere is observing the planet at optical wavelengths (where thermal emission is negligible) during eclipse, when the planet is near full phase and passes behind its host star. This method requires at least an order of magnitude higher photometric precision than transit observations of the same planet because the planet will be fainter than its host star, while the occulted area remains the same.

Observing an atmosphere at different orbital phases can provide further information about the scattering particles (e.g., Demory et al. 2013; Esteves et al. 2013; Heng & Demory 2013; Garcia Munoz & Isaak 2015; Shporer & Hu 2015; Oreshenko et al. 2016). Parmentier et al. (2016) suggested a connection between reflected light phase curve measurements and a sequence of condensate cloud models, but this only covered temperatures up to ${T}_{\mathrm{eq}}\sim 2200$ K: well below the equilibrium temperature of WASP-12b (${T}_{\mathrm{eq}}=2580$ K; Collins et al. 2017).

WASP-12b orbits a G0V star with an orbital period of 1.09 days (Hebb et al. 2009). While the host star is fairly faint (V = 12), WASP-12b's close semimajor axis and large radius (a = 0.0234 au, ${R}_{p}=1.90\,{R}_{J}$, ${R}_{p}=0.19\,{R}_{* }$; Collins et al. 2017) make it an excellent target for detailed study. Transit observations of WASP-12b range from 0.3 to 4.5 μm, and eclipse observations range from 0.9 to 8.0 μm (e.g., Hebb et al. 2009; López-Morales et al. 2010; Campo et al. 2011; Madhusudhan et al. 2011; Cowan et al. 2012; Crossfield et al. 2012; Copperwheat et al. 2013; Föhring et al. 2013; Sing et al. 2013; Swain et al. 2013; Stevenson et al. 2014a, 2014b; Croll et al. 2015; Sing et al. 2016). This work presents the first optical eclipse measurement of WASP-12b.

The atmospheric composition of WASP-12b has been extensively studied (e.g., Madhusudhan et al. 2011; Crossfield et al. 2012; Swain et al. 2013; Stevenson et al. 2014b), with initial claims of a C/O ratio greater than unity. This was first challenged by Crossfield et al. (2012) and Cowan et al. (2012), who instead reported an isothermal photosphere for WASP-12b. The recent detection of water in the planet's atmosphere has now firmly refuted the carbon-rich hypothesis (Kreidberg et al. 2015). Sing et al. (2013) found that the best-fit model for WASP-12b transmission spectroscopy was Mie scattering by an aluminum-oxide (Al₂O₃) haze. Barstow et al. (2017) found that an optically thick Rayleigh scattering aerosol with a 0.01 mbar top pressure best described the transmission observations, but the model poorly described the steep increase in transit depth at optical wavelengths. Schwartz et al. (2017) used thermal phase variations and eclipse depths to determine a Bond albedo of ${A}_{B}={0.2}_{-0.12}^{+0.1}$ and a dayside effective temperature of ${T}_{\mathrm{day}}=2864\pm 15$ K.

On 2016 October 19, a single eclipse of WASP-12b was observed with five HST orbits, using the STIS G430L grating (290–570 nm). The first HST orbit has significantly worse systematics than the four later orbits as a result of the repointing of the telescope, so these data were removed from the subsequent analysis. This left two HST orbits out of eclipse (one before and one after) when the planet and host star were both visible with the planet near full phase, as well as two HST orbits during eclipse when the planet was behind its host star, leaving only the star's light visible. These observations were granted as a part of programme GO-14797 (PI: Crossfield).

We used the same data collection method as previously used for similar observations (Sing et al. 2011, 2013, 2016; Evans et al. 2013). The subarray readout mode with a wide $52^{\prime\prime} \times 2^{\prime\prime}$ slit was used to minimize time-varying slit losses; this produced 1024 × 128 pixel images. In previous HST/STIS observations, the first frame from each HST orbit had systematically lower counts, so a 1 s dummy exposure was obtained at the beginning of each orbit, which successfully mitigated this systematic effect. This dummy exposure was then followed by 10 science exposures lasting 279 s each (the maximum recommended duration to avoid excessive cosmic-ray hits). Our final, analyzed data set thus contains 40 exposures collected over 331 minutes.

The raw STIS data were reduced (bias-, dark-, and flat-corrected) using the latest version of the CALSTIS1 pipeline and the relevant up-to-date calibration frames. Cosmic-ray events were identified and removed following Nikolov et al. (2014), as were all pixels identified as "bad" by CALSTIS. Overall, ∼9% of the pixels in each 2D spectrum were affected by cosmic rays with another ∼5% identified as "bad," resulting in a total of ∼14% interpolated pixels.

Next, the IRAF procedure apall was used to extract spectra from the calibrated .flt science files. We tested apertures between 9.0 and 17.0 pixels in intervals of 2 pixels and found that an 11.0 pixel aperture resulted in the lowest light curve residual scatter after fitting the white light data. However, the difference between apertures was minute (∼1 ppm). We then used cross-correlation to correct for subpixel shifts along the dispersion axis. The x1d files from CALSTIS were then used to calibrate the wavelength axis. Finally, both "white light" and six spectral channel light curves were produced by integrating the appropriate flux from each bandpass.

WASP-12b's host star WASP-12A is also orbited by two M-dwarf companions bound in a binary system 106 away from WASP-12A (Bergfors et al. 2011; Sing et al. 2013; Bechter et al. 2014). For our observations, the spectrograph slit orientation was chosen to be perpendicular to the line connecting WASP-12A and WASP-12(B,C) to allow maximal separation in the spatial direction of the resulting FITS files. The spectrum of the stellar companions is visually distinguishable from WASP-12A in the raw spectra and does not fall within our small spatial-axis aperture.

The top panel of Figure 1 shows the raw light curve binned across the entire STIS G430L bandpass ("white light"). There is a strong, repeated trend in flux, with exposures from each orbit appearing to follow a roughly polynomial trend. This systematic is well known and is believed to be the result of the thermal cycle of HST throughout its orbit as well as the movement of the spectral trace on the detector (e.g., Brown et al. 2001; Sing et al. 2011; Huitson et al. 2012; Evans et al. 2013).

Zoom In Zoom Out Reset image size

These systematic trends are also observed during HST/STIS observations of planetary transits, and a standard approach to remove them is assuming polynomial variations as a function of auxiliary variables (e.g., Sing et al. 2011; Huitson et al. 2012). More recently, Gibson et al. (2011) used Gaussian processes (GPs) to model the HST systematics, as the choice of polynomial model can potentially bias the results. For this reason, we modeled the systematics with the GP library george (Foreman-Mackey 2015) and used the same method as Evans et al. (2013). A detailed discussion of modeling systematics with GPs can be found in Gibson et al. (2012a, 2012b, 2013). We also attempted to fit the systematic variations with a polynomial model, which gave results consistent with our GP model.

3.1. Gaussian Process Model

The likelihood of a GP model is described as a multivariate normal distribution with

$$\begin{eqnarray}&&p({\boldsymbol{f}}| {\boldsymbol{X}},{\boldsymbol{\theta }},{\boldsymbol{\Omega }},{\boldsymbol{t}})={ \mathcal N }({\boldsymbol{E}}({\boldsymbol{\Omega }},{\boldsymbol{t}}),{\boldsymbol{\Sigma }}({\boldsymbol{X}},{\boldsymbol{\theta }})),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{f}}$ is the 40 measured fluxes, ${\boldsymbol{E}}$ is the eclipse function, and ${\boldsymbol{\Sigma }}$ is the kernel (covariance matrix). The time at the midpoint of each exposure is represented by ${\boldsymbol{t}}$. Further, ${\boldsymbol{X}}={[{\boldsymbol{\phi }},{\boldsymbol{\psi }}]}^{T}$ is the matrix of covariates, where ${\boldsymbol{\phi }}$ is the orbital phase of HST, and ${\boldsymbol{\psi }}$ is the slope of the spectral trace on the detector (computed using IRAF's apall procedure). These two covariates were selected as they provided the lowest scatter in the residuals after calibration. We also tested the inclusion of two additional covariates: the y-intercept of the spectral trace on the detector and the measured shifts of the spectral trace along the dispersion axis. However, the inclusion of these additional covariates did not significantly impact our results or uncertainties, likely because the covariates themselves are significantly correlated with the other covariates.

Our eclipse parameters are given by ${\boldsymbol{\Omega }}={[\alpha ,\delta ,\beta ]}^{T}$, where α is the baseline flux consisting of light emitted from both the planet and star, δ is the fractional eclipse depth ($\delta \,={F}_{\mathrm{planet}}/{F}_{\mathrm{star}}$), and β describes a constant rate of change in the baseline flux over time. Since we did not observe during eclipse ingress or egress, we used a boxcar function to describe the eclipse signal, with

$$\begin{eqnarray}{E}_{i} & = & \alpha (1-\delta {B}_{i})(1+\beta ({t}_{i}-{t}_{0}))\\ {B}_{i} & = & \left\{\begin{array}{ll}0 & i\in \mathrm{Orbit}\ 2\ \mathrm{or}\ 5\\ 1 & i\in \mathrm{Orbit}\ 3\ \mathrm{or}\ 4,\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where t₀ is the time of the first exposure.

Our GP parameters are given by ${\boldsymbol{\theta }}={[C,{L}_{\phi },{L}_{\psi },{\sigma }_{w}]}^{T}$, where C² is the maximum covariance, L_ϕ and L_ψ are covariance lengthscales, and ${\sigma }_{w}$ is the white noise level. We adopted the squared-exponential kernel:

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{nm}}={C}^{2}\exp \left[-\displaystyle \sum _{i=0}^{1}\displaystyle \frac{{({{\boldsymbol{X}}}_{{in}}-{{\boldsymbol{X}}}_{{im}})}^{2}}{{L}_{i}^{2}}\right]+{\delta }_{{nm}}{\sigma }_{w}^{2},\end{eqnarray} \tag{ 3 }$$

where ${L}_{i}={[{L}_{\phi },{L}_{\psi }]}_{i}$ and ${\delta }_{{nm}}$ is the Kronecker delta function. This kernel can be simply understood as requiring that observations be strongly correlated if they have similar spectral trace slope and HST orbital phase, while observations further from each other in covariate space are more weakly correlated. This then describes a smoothly varying function of the covariates, with the addition of white noise.

The final model is then given by

$$\begin{eqnarray}&&{f}^{* }={\boldsymbol{\mu }}({\boldsymbol{\phi }},{\boldsymbol{\psi }})+{\boldsymbol{E}}({\boldsymbol{\Psi }}),\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{\mu }}({\boldsymbol{\phi }},{\boldsymbol{\psi }})$ is the GP model mean.

The Markov Chain Monte Carlo (MCMC) ensemble sampler software emcee (Foreman-Mackey et al. 2013) was used to explore the seven parameters, determining the most likely parameter values and their uncertainties. For computational reasons, the variables used in this MCMC were {δ, $\mathrm{ln}(\alpha )$, β, $\mathrm{ln}({L}_{\phi })$, $\mathrm{ln}({L}_{\psi })$, $\mathrm{ln}({\sigma }_{w}^{2})$, and $\mathrm{ln}({C}^{2})$}. Using logarithms removes the need to use a prior to obtain strictly positive values. While the eclipse depth, δ, should be strictly positive, we allowed for negative values to ensure an unbiased estimate. A uniform prior was used so $\mathrm{ln}({L}_{\phi })\lt 0$ and $\mathrm{ln}({L}_{\psi })\lt 0$, which has the effect of ensuring that these lengthscales are within a few orders of magnitude of the variations in the covariates. For un-normalized white light data, the best-fit values from a 10⁶ step MCMC chain were {$\delta =(-5.3\pm 7.4)\times {10}^{-5}$, $\alpha =({4.198}_{-0.008}^{+0.012})\times {10}^{7}$, $\beta =0.0056\pm 0.0015$, ${L}_{\phi }={0.05}_{-0.03}^{+0.19}$, ${L}_{\psi }={0.03}_{-0.03}^{+0.27}$, ${\sigma }_{w}^{2}\,=({5.2}_{-1.3}^{+1.7})\times {10}^{7}$, and ${C}^{2}=({8}_{-5}^{+14})\times {10}^{7}$}.

The top panel of Figure 1 shows the median model and uncertainty from a 10⁶ step MCMC chain overplotted on the raw white light flux measurements. The bottom panel of Figure 1 shows the light curve produced by dividing the raw spectra by the median model (excluding the change in flux during eclipse), with the median eclipse model overplotted. The clear linear trend in the calibrated flux of HST orbit #5 (bottom panel of Figure 1) shows that there is still substantial correlated noise in the data that could not be described by any of the four considered covariates.

The STIS G430L spectra were binned into six spectral channels to allow moderate wavelength resolution while keeping uncertainties on each channel sufficiently small to be able to test atmospheric models. Each spectral channel was modeled independently using the GP method described above. Light curves after GP calibration are shown for each spectral channel in Figure 2, and the relevant results are tabulated in Table 1. Eclipse depths were found using the median value from a 10⁶ step MCMC chain, while the 84 and 97.5 percentiles were used to determine upper limits. The larger uncertainties in eclipse depth at shorter wavelengths are due to lower stellar flux and detector sensitivity.

Zoom In Zoom Out Reset image size

Table 1.  Eclipse Depths and Geometric Albedos

Wavelengths	Eclipse Depth, δ (ppm)		${\delta }_{\mathrm{thermal}}$	Geometric Albedo, Ag
(nm)	Best Fit	97.5% Upper Limit	(ppm)	Best Fit	97.5% Upper Limit
290–570	−53 ± 74	96	56	−0.035 ± 0.050	0.064
290–336	−60 ± 540	1020	10	−0.04 ± 0.36	0.68
336–383	90 ± 290	670	20	0.06 ± 0.20	0.45
383–430	−30 ± 180	330	40	−0.02 ± 0.12	0.22
430–476	−60 ± 130	210	60	−0.039 ± 0.089	0.14
476–523	−70 ± 130	190	100	−0.045 ± 0.087	0.13
523–570	−50 ± 150	240	160	−0.036 ± 0.098	0.16

Download table as:  ASCIITypeset image

Because WASP-12b is so strongly irradiated, the peak of its thermal emission is expected to be at ∼1 μm for a ∼2800 K dayside temperature (Schwartz et al. 2017). For this reason, we calculated the predicted eclipse depths due to thermal radiation from WASP-12b, assuming a T = 3000 K blackbody for WASP-12b (hotter than inferred from infrared observations due to the greater depth of the optical photosphere; Cowan & Agol 2011) and a standard G0V spectrum from Pickles (1998) for WASP-12. These depths (${\delta }_{\mathrm{thermal}}$) are summarized in Table 1 and are all within our 97.5% confidence interval upper limits.

If interpreted as solely due to reflected light, eclipse depths can be converted to geometric albedo using

$$\begin{eqnarray}&&{A}_{g}=\delta {\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{-2},\end{eqnarray} \tag{ 5 }$$

where R_p = 1.90 R_J is the radius of the planet, and $a\,=0.0234\,\mathrm{au}$ is its orbital semimajor axis (Collins et al. 2017). Applying Equation (5) to the best-fit eclipse depths and their corresponding upper limits gives constraints on the geometric albedo across the STIS G430L wavelength range (summarized in Table 1).

Our reported uncertainties on eclipse depths are ∼40% higher than the photon limit. The increased scatter in our data may be the result of incomplete modeling of the systematic noise. Alternatively, our uncertainties may be limited by intrinsic stellar variability. Given the slow rotation period of WASP-12 compared to the observing window (${P}_{\mathrm{rot}}\gtrsim 23$ days given $v\sin i\lt 2.2\,\mathrm{km}\,{{\rm{s}}}^{-1}$; Hebb et al. 2009), variability due to stellar rotation (e.g., starspots passing in and out of view) should not significantly affect our observations. However, our Sun's total irradiance (spatially and spectrally integrated) is known to vary at the 10⁻⁴ level on timescales of minutes to hours as a result of solar convection and oscillations (Kopp 2016). Given the G0V spectral class of WASP-12, similar variations may also be present and may explain the greater than Poisson limit uncertainties as well as the residual correlated noise in the calibrated time-series spectra.

We use the NEMESIS spectral retrieval tool (Irwin et al. 2008; Barstow et al. 2014) to produce predicted model spectra given two previously proposed models for WASP-12b: an Al₂O₃ haze and a cloud-free atmosphere. NEMESIS is not a radiative equilibrium code; rather, it takes an atmospheric model and calculates incident and emergent flux and will not take into account heating from incoming stellar radiation. The limits from our HST/STIS eclipse observations firmly reject both models; we find ${\chi }^{2}$ per datum (${\chi }^{2}/{N}_{\mathrm{obs}}$, ${N}_{\mathrm{obs}}=6$) of 41 and 10 for the Al₂O₃ haze and cloud-free models.

Given its exceedingly high equilibrium temperature (${T}_{\mathrm{eq}}=2580$ K; Collins et al. 2017), WASP-12b would technically lie within Sudarsky et al.'s (2000) Class V (${T}_{\mathrm{eff}}\gt 1500$ K) but is far hotter than any planet they considered. On the planet's dayside, WASP-12b is far too hot for condensates to form (Wakeford et al. 2017). However, temperatures near the planet's day–night terminator, and across the planet's nightside, may be cool enough to allow for the formation of condensates that could affect transmission spectroscopy without significantly affecting dayside eclipse spectroscopy.

Also, it is expected that Na i absorption (which is important at lower temperatures) will not contribute much to the low albedo of WASP-12b as most of the sodium will be ionized on the hot dayside. Instead, it is expected that the atmosphere will be dominated by Rayleigh scattering from atomic hydrogen and helium, with a small contribution from electron scattering. The red line in Figures 3 and 4 shows the predicted eclipse depth (binned to a resolution of 1 point per 5 nm) from Crossfield et al. (2012) made with the PHOENIX atmosphere code adapted for hot Jupiters as described in Barman et al. (2001, 2005). In this model, reflected light makes up $\lesssim 10$% of the eclipse depth (${A}_{g}\lesssim 0.002$) at the shortest wavelengths and $\ll 1$% of the eclipse depth at infrared wavelengths; the remainder of the eclipse depth is due to thermal emission. This model gives a ${\chi }^{2}$ per datum of 0.9 for our HST/STIS data (${N}_{\mathrm{obs}}=6$), but a worse ${\chi }^{2}$ per datum of 3 for all of the data plotted on Figure 4 (${N}_{\mathrm{obs}}=21$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are significant differences between the PHOENIX model and the cloud-free model produced by NEMISIS, including but not limited to the inclusion of atomic hydrogen opacities (lines and bound-free opacities), as well as the typical opacities more commonly associated with cool stellar photospheres. Also, the PHOENIX model results from a self-consistent calculation of the thermal structure, chemistry, line-by-line opacities (as well as scattering), and irradiation, thereby accounting for important changes that occur in the hot upper layers of WASP-12b (for example, the transition from H₂ to H at low pressures and the thermal ionization of Na and K).

Our observations cover the blackbody peak of WASP-12 (∼450 nm) and show that little of the incident radiation at these wavelengths is reflected by the planet. Geometric albedo is related to spherical albedo through a phase integral q such that A_s = qA_g, and Bond albedo is equal to the flux-weighted, wavelength-averaged spherical albedo. If we assume diffuse scattering (q = 1.5), our "white light" 97.5% confidence upper limit on the geometric albedo across the STIS bandpass (${A}_{g}\lt 0.064$) suggests $\lt 10 \%$ of the energy received at these wavelengths is reflected. However, since the wavelengths observed cover only 36% of the incident stellar energy, the Bond albedo is not well constrained by these measurements and is consistent with Schwartz et al.'s (2017) measurement of ${A}_{B}={0.2}_{-0.12}^{+0.1}$.

Our results are in stark contrast with those for the much cooler HD 189733b, the only other hot Jupiter with spectrally resolved reflected light observations (Evans et al. 2013); those data showed an increase in albedo with decreasing wavelength. The fact that the first two exoplanets with optical albedo spectra exhibit significant differences demonstrates the importance of spectrally resolved reflected light observations and highlights the great diversity among hot Jupiters.

T.J.B. acknowledges support from the McGill Space Institute Graduate Fellowship and from the FRQNT through the Centre de recherche en astrophysique du Québec. J.K.B. acknowledges support from the Royal Astronomical Society Research Fellowship. I.J.M.C. was supported under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement No. 336792. This work is based on observations made with the NASA/ESA Hubble Space Telescope, 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-26555. The Al₂O₃ and cloud-free models tested in this work were made with the NEMESIS code developed by Patrick Irwin. We have also made use of free and open-source software provided by the Matplotlib, Python, and SciPy communities.

Facility: HST(STIS) - Hubble Space Telescope satellite.