Near-IR Transmission Spectrum of HAT-P-32b using HST/WFC3

The spatially scanned spectroscopic images of HAT-P-32b were obtained with the G141 grism and are available from the MAST archive³ (ID:14260, PI: Deming Drake). The data set contains five consecutive HST orbits and each exposure is the result of 14 non-destructive reads, with a size of 256 × 256 pixels in the SPARS10 mode (exposure time = 88.435623 s). With this configuration the maximum signal level is $2.6\times {10}^{4}$ electrons per pixel and the total scan length is approximately 40 pixels.

During the light curve analysis, the first of the five orbits was discarded. This is a standard practice for exoplanet transit observations (e.g., Deming et al. 2013; Huitson et al. 2013; Haynes et al. 2015; Tsiaras et al. 2016a), as the telescope needs to stabilize into its new position. Of the remaining four HST orbits, the first and the fourth provide the out-of-transit baseline, while the second and the third capture the transit. The data set contains, for calibration purposes, a non-dispersed (direct) image of the target, obtained using the F139N filter.

Before extracting the light curves (white and spectral), all frames were reduced using the routines described in Tsiaras et al. (2016b). HAT-P-32A has an M1.5 stellar companion, HAT-P-32B (${T}_{\mathrm{eff}}=3565\pm 82\,{\rm{K}}$; Zhao et al. 2014). The dispersed signals from HAT-P-32A and HAT-P-32B are blended when using the scanning mode. However, these two stars are separated enough ($2\buildrel{\prime\prime}\over{.} 923\pm 0\buildrel{\prime\prime}\over{.} 004$, Zhao et al. 2014) to avoid blending when the differential reads (the difference between two consecutive non-destructive reads, or "stripes") are considered. For each stripe, we determined the photometric aperture, taking into account the wavelength-dependent photon trajectories (Tsiaras et al. 2016b) and obtained a set of 12 white light curves. The same criterion was used to extract the spectral light curves, obtaining a set of 12 time series for each one of the 20 spectral bins. The wavelength range of each bin was chosen in order to have a similar flux level across all bins.

It is known that instrumental systematics (known as "ramps") affect the WFC3 infrared detector both in staring (Berta et al. 2012; Swain et al. 2013; Wilkins et al. 2014) and scanning modes (Deming et al. 2013; Knutson et al. 2014; Kreidberg et al. 2014; Tsiaras et al. 2016a, 2016b). The brighter the star, the stronger the ramps. In the case of HAT-P-32b, the host star is relatively faint (${K}_{\mathrm{mag}}=9.99$) so we did not expect very strong ramps.

We fitted the ramps on the white light curve using a similar approach to Kreidberg et al. (2014); i.e., we adopted an analytic function with two different types of ramps, short-term and long-term, to correct the data:

$$\begin{eqnarray}&&R(t)=(1-{r}_{a}(t-{t}_{v}))(1-{r}_{b1}{e}^{-{r}_{b2}(t-{t}_{0})})\end{eqnarray} \tag{ 1 }$$

where, t is the mid-time of each exposure, t_v is the time when the visit starts, t₀ is the time when each orbit starts, r_a is related to the long-term ramp, and ${r}_{b1},{r}_{b2}$ are related to the short-term ramp.

To model the transit light curve we used our Python package, PyLightcurve,⁴ which returns the flux as a function of time using the nonlinear limb-darkening law (Claret 2000). The limb-darkening coefficients were fitted on the profile of a star similar to HAT-P-32A (T_* = 6207 K, [Fe/H] = −0.04 dex, $\mathrm{log}({g}_{* })=4.33$ [cgs]), using a modified version of the ATLAS stellar model described in Howarth (2011). In this fit we took into account the variable sensitivity of the G141 grism across its wavelength range. The observations do not cover both the ingress and the egress of the transit, hence we could not fit for the semi-major axis and inclination, which have been fixed to the values reported in Table 1. We also assumed a circular orbit. The results are shown in Figure 1 and reported in Table 2. As can be seen in the residuals near the egress, there are a few points that appear to deviate significantly from the model fitted to the white light curve. The root mean square (rms) of the residuals is 180 ppm, significantly higher than the photon-noise limited rms of 100 ppm. The error bar in transit depth is accordingly higher than the photon-noise limited case. This behavior could be caused by star-spots, which notoriously might generate a wavelength-dependent astrophysical signal and therefore a distortion on the spectrum. However, we repeated the analysis excluding these points and the spectrum was not affected. We also tested changing the orbital parameters up to 1σ, in both directions, from their reference values. However, this did not improve the white light curve residuals. The main effect was a shift of ∼120 ppm in the white light curve transit depth—i.e., up to 1.3σ from our reported uncertainty—but again no detectable effect on the spectrum (differential transit depths vary by less than 0.25σ on average).

Zoom In Zoom Out Reset image size

Table 1.  Parameters of the HAT-P-32b System (Hartman et al. 2011)

Stellar Parameters
[Fe/H] [dex]	−0.04 ± 0.08
Teff [K]	6207 ± 88
${M}_{* }\,[{M}_{\odot }]$	1.160 ± 0.041
${R}_{* }\,[{R}_{\odot }]$	1.219 ± 0.016
$\mathrm{log}({g}_{* })$ [cgs]	4.33 ± 0.01
Planetary Parameters
Teq [K]	1786 ± 26
${M}_{{\rm{p}}}\,[{M}_{\mathrm{Jup}}]$	0.860 ± 0.164
${R}_{{\rm{p}}}\,[{R}_{\mathrm{Jup}}]$	1.789 ± 0.025
a [AU]	0.0343 ± 0.0004
Transit Parameters
T0 [BJD]	2454420.44637 ± 0.00009
Period [days]	2.150008 ± 0.000001
${R}_{{\rm{p}}}/{R}_{* }$	0.1508 ± 0.0004
$a/{R}_{* }$	${6.05}_{-0.04}^{+0.03}$
i [deg]	88.9 ± 0.4

Download table as:  ASCIITypeset image

Table 2.  White Light Curve Fitting Results

Limb-darkening Coefficients (1.125–1.650 μm)
a1	0.603336
a2	−0.223032
a3	0.281379
a4	−0.13988
Fitting Results
${T}_{0}\,(\mathrm{HJD})$	2457408.95783 ± 0.00004
${R}_{{\rm{p}}}/{R}_{* }$	0.1521 ± 0.0003

Download table as:  ASCIITypeset image

Finally, for each wavelength bin we divided the spectral light curve by the white light curve (Kreidberg et al. 2014) and fitted a linear trend simultaneously with a relative transit model:

$$\begin{eqnarray}&&{n}_{\lambda }(1+{\chi }_{\lambda })({F}_{\lambda }/{F}_{W})\end{eqnarray} \tag{ 2 }$$

where ${n}_{\lambda }$ is the normalization factor that needs to be calculated for each bin, ${\chi }_{\lambda }$ is the wavelength-dependent linear ramp (Tsiaras et al. 2016b, 2016a), and $({F}_{\lambda }/{F}_{W})$ is the ratio between the spectral light curve and the white light curve. We fitted this model using the same orbital parameters listed in Table 1 and the white ${R}_{P}/{R}_{* }$ ratio obtained from the white light curve fitting. The limb-darkening coefficients were calculated for each bin using the same method as for the white light curve (see Table 4). The rms of the residuals for the spectral light curves (on average 474 ppm) is close to the photon-noise limited case (on average 443 ppm). The corresponding error bars in relative transit depths are also ∼10% above the photon noise limit. This proves that the deviation from the model seen for the white light curve is not wavelength dependent.

The planetary spectrum was fitted using MCMC and following two different approaches, leading to a "stacked" and a "weighted" spectrum:

• 1.  
  (stacked) using a unique reference light curve for each spectral bin, obtained by summing the relative stripe light curves;
• 2.  
  (weighted) fitting each stripe light curve alone, then taking the weighted mean for each spectral bin.

Following the first method, we obtained the white light curve shown in Figure 1 (top panel). Both methods give the same modulation, with the exception of a few bins where the differences are within 0.3σ.

After the preliminary reduction, we obtained 12 stripe light curves for each of the 20 spectral bins, as described in Section 2.2. We performed ICA for all bins separately, by using the corresponding 12 stripe light curves as input time series (vector ${\boldsymbol{x}}$ in Equation (3)). Similarly, the light curves integrated over the 20 spectral bins for each stripe were used as input white light curves. Thus, we obtained one set of components for each spectral bin and an additional set for the whole spectral range. The transit signal is mainly contained in the first components of all sets, while the other components are predominantly instrument systematics and noise.

Following a standard ICA algorithm (e.g., Morello 2015; Morello et al. 2016), we simultaneously fitted a transit model (with the same parameters as in Section 2.3) and a linear combination of the non-transit components to the relevant raw light curves. We computed the "stacked" and a "weighted" spectra, as described in Section 2.3.

The residuals obtained for the stacked white light curve have been included as an additional component in the spectral fits. This step is equivalent to dividing by the white light curve as is done in the parametric fitting (see Section 2.3), in order to remove possible undetrended systematics common to all wavelengths.

The fitting process is as follows. First, we run a Nelder-Mead optimization algorithm to find the parameter values minimizing the fitting residuals, then we use them as starting values for a Markov Chain Monte Carlo (MCMC) calculation with 300,000 iterations. The likelihood's variance, ${\sigma }_{0}^{2}$, is initialized to the variance of the residuals, then updated at any iteration. The best-fitting parameters are estimated as ${\mu }_{\mathrm{par}}\pm {\sigma }_{\mathrm{par},0}$, where ${\mu }_{\mathrm{par}}$ and ${\sigma }_{\mathrm{par},0}$ are the mean value and standard deviation of the relevant parameter chain, respectively.

To fully account for the potential bias associated with the detrending technique, the final error bars are re-scaled with respect to the MCMC error bars inferred from the residuals only, by adding a ${\sigma }_{\mathrm{ICA}}^{2}$ term to the likelihood's variance:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{par}}=\sqrt{\displaystyle \frac{{\sigma }_{\mathrm{ICA}}^{2}+{\sigma }_{0}^{2}}{{\sigma }_{0}^{2}}}{\sigma }_{\mathrm{par},0}.\end{eqnarray} \tag{ 5 }$$

The ${\sigma }_{\mathrm{ICA}}^{2}$ term is calculated as:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ICA}}^{2}=\displaystyle \sum _{j}{o}_{j}^{2}{\mathrm{ISR}}_{j}\end{eqnarray} \tag{ 6 }$$

where ISR is the so-called interference-to-signal-ratio matrix computed with ICA, and o_j are the coefficients of the non-transit components. Simply put, the ${\sigma }_{\mathrm{ICA}}$ term is the weighted sum of the errors attributed to the independent components extracted with ICA. We refer the reader to Morello et al. (2015, 2016) for additional details.

The error bars for the weighted spectrum are calculated as the simple arithmetic means of the error bars derived from fitting the independent components to the single stripes. These are worst-case estimates, as they do not scale when combining the results from the stripes. Scaling the error bars would not be theoretically correct, as the individual fits are not independent, given that they adopt the same components, which are estimated using the information contained in all of the stripes. The error bars obtained with ICA are larger than the ones obtained with the parametric approach by a factor ∼1.6 (weighted) and ∼1.8 (stacked). Note that scaling the error bars in the weighted approach would have led to final error bars smaller than photon-noise limited.

The transmission spectrum of HAT-P-32b and the best fit to it, retrieved with ${ \mathcal T }$-REx, is shown in Figure 3. The best-fitting values and the posterior distributions are shown in Table 4 and Figure 4. With the exception of water vapor, the fitted values for all of the other molecular mixing ratios are smaller than 10⁻⁷. This result means that they are not detectable from this data set. The water vapor mixing ratio oscillates, instead, between ${\mathrm{log}{\rm{H}}}_{2}{\rm{O}}=-{3.45}_{-1.65}^{+1.83}$ depending on the clouds' top pressure, which could occur between 5.16 and 1.73 bar. A strong correlation between the water vapor mixing ratio, the clouds' top pressure, planetary radius at 10 bar, and temperature is noticeable in Figure 4, indicating there is a degeneracy of solutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Limb-darkening Coefficients ${a}_{1-4}$ and Transit Depth ${({R}_{p}/{R}_{* })}^{2}$ for the Wavelength Channels

${\lambda }_{1}-{\lambda }_{2}\,(\mu {\rm{m}})$		a1	a2	a3	a4	${({R}_{{\rm{p}}}/{R}_{* })}^{2}\,(\mathrm{ppm})$	${({R}_{{\rm{p}}}/{R}_{* })}^{2}\,(\mathrm{ppm})$
						UCL Pipeline	Stripe-ICA
1.1250	1.1511	0.632741	−0.481904	0.701108	−0.306091	22940 ± 112	22961 ± 184
1.1511	1.1767	0.619205	−0.434713	0.64011	−0.282483	22862 ± 100	22890 ± 184
1.1767	1.2011	0.614294	−0.41589	0.610565	−0.272242	23091 ± 105	23057 ± 181
1.2011	1.2247	0.599151	−0.360648	0.544934	−0.247917	23083 ± 105	23163 ± 186
1.2247	1.2480	0.584001	−0.29953	0.465487	−0.216442	22893 ± 102	22926 ± 179
1.2480	1.2716	0.581928	−0.282551	0.441745	−0.210655	22878 ± 102	22863 ± 242
1.2716	1.2955	0.58946	−0.229732	0.322997	−0.169253	22951 ± 110	22966 ± 189
1.2955	1.3188	0.57237	−0.227002	0.362724	−0.181489	22864 ± 123	22911 ± 200
1.3188	1.3421	0.569522	−0.202303	0.325228	−0.166816	23176 ± 94	23188 ± 203
1.3421	1.3657	0.564634	−0.163366	0.265035	−0.14235	23335 ± 129	23381 ± 189
1.3657	1.3901	0.561817	−0.127278	0.200548	−0.113503	23255 ± 103	23285 ± 236
1.3901	1.4152	0.561832	−0.0979712	0.148201	−0.0914278	23122 ± 111	23064 ± 182
1.4152	1.4406	0.572262	−0.100901	0.133369	−0.0848254	23382 ± 119	23396 ± 195
1.4406	1.4667	0.58462	−0.111943	0.124656	−0.0799948	23202 ± 115	23129 ± 196
1.4667	1.4939	0.600205	−0.136878	0.140204	−0.0874595	23181 ± 130	23170 ± 294
1.4939	1.5219	0.609784	−0.134319	0.11158	−0.0721681	23041 ± 160	22987 ± 204
1.5219	1.5510	0.626375	−0.139701	0.0839621	−0.0555132	22890 ± 114	22959 ± 201
1.5510	1.5819	0.647904	−0.193435	0.120068	−0.0635888	23076 ± 131	22978 ± 230
1.5819	1.6145	0.663831	−0.223633	0.124246	−0.0583813	22871 ± 102	22886 ± 171
1.6145	1.6500	0.686226	−0.267069	0.137329	−0.0557593	22611 ± 111	22680 ± 198

Download table as:  ASCIITypeset image

Table 4.  Fitting Results for HAT-P-32b Atmosphere

Parameter	Value
${\mathrm{log}{\rm{H}}}_{2}{\rm{O}}$	$-{4.66}_{-1.93}^{+1.66}$
${T}_{\mathrm{eff}}\,[{\rm{K}}]$	${1553}_{-91}^{+174}$
${R}_{p}\,[{R}_{\mathrm{jup}}]$	$-{1.76}_{-0.04}^{+0.05}$
${P}_{\mathrm{cld},\mathrm{top}}\,[\mathrm{bar}]$	${3.39}_{-1.66}^{+1.77}$

Download table as:  ASCIITypeset image

As mentioned in the previous section, the error bars obtained with ICA are larger by a factor of ∼1.6–1.8 compared to the ones obtained with the parametric fitting. The larger error bars obtained with ICA are the trade-off for higher objectivity, due to the lack of any assumption about the instrument systematics compared to the parametric approach. The ICA error bars are worst-case estimates. It is worth noting that the discrepancies between the spectra obtained with the different methods are smaller than the parametric error bars, suggesting that, in this case, the ICA error bars might be overly conservative.

Previous ground-based observations of the transit of HAT-P-32b in the optical wavelengths (Gibson et al. 2013; Zhao et al. 2014; Mallonn & Strassmeier 2016; Nortmann et al. 2016) did not find evidence of spectral modulations due to molecules. Our cloud top pressure is consistent with their measurements within 1σ, hence the water detection in the infrared is not controversial.

Water vapor has been detected, to date, in the atmospheres of about 10 hot Jupiters (Iyer et al. 2016). Stevenson (2016) identified two classes of hot Jupiters, essentially mostly cloudy or with a strong water signature. The observed trend suggests that hotter (${T}_{\mathrm{eq}}\gt 700$ K) and more inflated ($\mathrm{log}g\gt 2.8$) planets are more likely to have a strong water signature than cooler and smaller ones, but the current sample is not statistically significant. In agreement with this scenario, we find that HAT-P-32b (${T}_{\mathrm{eq}}=1786$ K; $\mathrm{log}g\gt 2.8$) has one of the strongest water features so far detected (∼500 ppm, 5.3σ).