Near-IR transmission spectrum of HAT-P-32 b using HST/WFC3
M. Damiano G. Morello A. Tsiaras T. Zingales G. Tinetti

	
	
	
	Mario Damiano
	mario.damiano.15@ucl.ac.uk
	
	0000-0002-1830-8260]M. Damiano
	Department of Physics & Astronomy, University College London, Gower Street, WC1E6BT London, United Kingdom
	INAF - Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, I-90134 Palermo, Italy
	
	
	Department of Physics & Astronomy, University College London, Gower Street, WC1E6BT London, United Kingdom
	
	
	Department of Physics & Astronomy, University College London, Gower Street, WC1E6BT London, United Kingdom
	
	
	Department of Physics & Astronomy, University College London, Gower Street, WC1E6BT London, United Kingdom
	INAF - Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, I-90134 Palermo, Italy
	
	
	Department of Physics & Astronomy, University College London, Gower Street, WC1E6BT London, United Kingdom
	
	
	
		We report here the analysis of the near-infrared transit spectrum of the hot-Jupiter  which was recorded with the Wide Field Camera 3 (WFC3) on-board the Hubble Space Telescope (HST).  is one of the most inflated exoplanets discovered, making it an excellent candidate for transit spectroscopic measurements. To obtain the  transit spectrum, we have adopted  different analysis methods, both parametric and non parametric (Independent Component Analysis, ICA), and compared the results. The final spectra are all consistent within 0.5σ. The uncertainties obtained with ICA are larger than those obtained with the parametric method by a factor ∼1.6 - 1.8. This difference is 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 therefore worst-case estimates. To interpret the spectrum of   we used , our fully Bayesian spectral retrieval code. As for other hot-Jupiters, the results are consistent with the presence of water vapor (logH_2O = -3.45_-1.65^+1.83), clouds (top pressure between 5.16 and 1.73 bar). Spectroscopic data over a broader wavelength range will be needed to de-correlate the mixing ratio of water vapor from clouds and identify other possible molecular species in the atmosphere of .
	
	
	
	
	
	
	

§ INTRODUCTION
 
	In the past decade the Hubble Space Telescope has been an invaluable observatory to study the properties of exoplanetary atmospheres. The majority of the planets observed to date are hot and gaseous as they are the easiest targets to probe. Transit observations in the UV, VIS and IR have started to provide important insights into the chemical composition and structure of the atmospheres of gas-giants orbiting very close to their star. Many of these atmospheres appear to be in the hydrodynamic escape regime given their vicinity to the stellar host <cit.>. Common atmospheric components detected include alkali metals <cit.> and water vapor <cit.>. Condensates or hazes have also been identified <cit.>. Some of the data also suggest that carbon-bearing or more exotic species, such as TiO and VO <cit.>, are present in some of these atmospheres. Finally, eclipse and phase curve observations have enabled to glimpse into the atmospheric thermal properties and global circulation of a few of these objects <cit.>.
	
	In this work we analyze the near-infrared transit spectrum of the hot-Jupiter  (T_eq = 1786 K) <cit.> obtained with the WFC3 camera on-board the HST.  is one of the most inflated exoplanets discovered, being less massive than Jupiter (M_p = 0.79   M_Jup) but having almost twice its radius (R_p = 1.789   R_Jup). The atmosphere of  has been observed with ground-based instruments in the optical wavelengths, revealing a featureless transmission spectrum <cit.>. In addition, <cit.> suggested the presence of a thermal inversion in the atmosphere of  to interpret eclipse observations.
	
	We used our dedicated WFC3 pipeline <cit.> to extract the transit light-curves per wavelength channel and obtain the planetary spectrum (Section <ref>). We used in parallel Independent Component Analysis to correct for the instrumental systematics, and investigate the effect of different analysis techniques on the same data set (Section <ref>). The final spectrum was analyzed using our fully Bayesian spectral retrieval code,  <cit.>.
	
	
	
	
	

§ DATA ANALYSIS
 
	
	
	

 §.§ Observations
 
	
	The spatially scanned spectroscopic images of  were obtained with the G141 grism and are available from the MAST archive[<https://archive.stsci.edu/>] (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× 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 <cit.>, 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.
	
	

 §.§ Extraction of light-curves
 
	
	Before extracting the light-curves (white and spectral), all frames were reduced using the routines described in <cit.>. HAT-P-32 A has an M1.5 stellar companion, HAT-P-32 B <cit.>. The dispersed signals from HAT-P-32 A and B are blended when using the scanning mode. However, these two stars are separated enough <cit.> 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 <cit.> 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.
	
	
	
	

 §.§ Parametric fitting

	
	It is known that instrumental systematics (known as “ramps”) affect the WFC3 infrared detector both in staring <cit.> and scanning modes <cit.>. The brighter is the star, the stronger are the ramps. In the case of , the host star is relatively faint (K_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 <cit.>, i.e. we adopted an analytic function with two different types of ramps, short-term and long-term, to correct the data:
	
	
    R(t) = (1-r_a(t-t_v))(1-r_b1e^-r_b2(t-t_0))

	
	where, t is the mid-time of each exposure, t_v is the time when the visit starts, t_0 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[<https://github.com/ucl-exoplanets/pylightcurve>], which returns the flux as a function of time using the non-linear limb darkening law <cit.>. The limb darkening coefficients were fitted on the profile of a star similar to HAT-P-32 A (T_* = 6207 K, [Fe/H] = -0.04 dex, log(g_*) = 4.33 [cgs]), using a modified version of the ATLAS stellar model described in <cit.>. 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 <ref>. We also assumed a circular orbit. The results are shown in Figure <ref> and reported in Table <ref>. As we can see 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 while 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).
	
	
	
	Finally, for each wavelength bin we divided the spectral light-curve by the white light-curve <cit.> and fitted a linear trend simultaneously with a relative transit model:
	
	
    n_λ(1+χ_λ)(F_λ/F_W)

	
	where n_λ is the normalization factor that needs to be calculated for each bin, χ_λ is the wavelength-dependent linear ramp <cit.>, (F_λ/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 <ref> 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 <ref>). 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 is 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:
	
		
  * (stacked) using a unique reference light-curve for each spectral bin, obtained by summing the relative stripe light-curves;
		
  * (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 <ref> (top panel). Both methods give the same modulation with the exception of a few bins where the differences are within 0.3 σ.
	
	
	

§ ICA
 
	Independent Component Analysis (ICA) is a blind signal-source separation (BSS) technique which is able to separate the source signals in a set of observations without any prior knowledge about the signals themselves or their mixing ratios. In many applications, observations are well-represented as linear combinations of certain (unknown) source signals:
	
    x = As

	where x = (x_1,x_2,...,x_n)^T is the column vector of observed signals, s = (s_1,s_2,...,s_n)^T is the column vector of source signals, and A is the so-called mixing matrix. The original source signals are retrieved through a linear transformation that maximizes their mutual independence, according to one or more statistical estimators <cit.>:
	
    s = Wx

	ICA has been used to remove instrument systematics and other astrophysical signals in exoplanetary light-curves obtained with Kepler, HST/NICMOS <cit.>, Spitzer/IRS <cit.>, and Spitzer/IRAC <cit.> with excellent results. We refer the reader to those publications and the relevant cited literature for more technical details about ICA and the different implementations. In this paper, we discuss a similar approach to the analysis of spectroscopic time series obtained with HST/WFC3 using the scanning-mode technique.
	The main steps of the algorithm are:
	
		
  * ICA decomposition;
		
  * Fitting;
		
  * Finalizing the parameter error bars.
	
	

 §.§ ICA decomposition

	After the preliminary reduction, we obtained 12 stripe light-curves for each of the 20 spectral bins, as described in Section <ref>. We performed ICA for all bins separately, by using the corresponding 12 stripe light-curves as input time series (vector x in Equation <ref>). 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.
	

 §.§ Fitting

	Following a standard ICA algorithm (e.g. <cit.>), we simultaneously fitted a transit model (with the same parameters as in Section <ref>) 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 <ref>.
	
	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 <ref>), 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, σ_0^2, is initialized to the variance of the residuals, then updated at any iteration. The best fitting parameters are estimated as μ_par±σ_par,0, where μ_par and σ_par,0 are the mean value and standard deviation of the relevant parameter chain, respectively.
	
	

 §.§ Final error bars

	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 σ_ICA^2 term to the likelihood's variance:
	
    σ_par = √(σ_ICA^2 + σ_0^2/σ_0^2)σ_par,0

	The σ_ICA^2 term is calculated as:
	
    σ_ICA^2 = ∑_j o_j^2 ISR_j

	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. In plain words, the σ_ICA term is the weighted sum of the errors attributed to the independent components extracted with ICA. We refer the reader to <cit.> 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 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 lead to final error bars smaller than photon noise limited.
	
	
	

§ ATMOSPHERIC RETRIEVAL

	To interpret the spectrum of , we use  <cit.>, a Bayesian spectral retrieval code which uses line lists provided by ExoMol <cit.>, HITRAN <cit.> and HITEMP <cit.>.
	We assumed an atmosphere dominated by  molecular Hydrogen and Helium, with a mean molecular weight of 2.3 amu. We considered as candidate trace gases a broad range of molecules, including H_2O, C_2H_2, CH_4, CO_2, CO, HCN, NH_3, VO and TiO. The RobERt (Robotic Exoplanet Recognition, ) module restricts the list of detectable molecules, based on the observed spectral pattern, to H_2O, TiO and VO. Given the relatively narrow spectral range probed, we assumed an isothermal profile and molecular abundances constant with pressure. In addition, we set uniform priors to the fitted parameters, which were: the mixing ratios of the molecules (10^-12 - 10^-1), the effective temperature of the planet (1400 - 2100 K), the radius of the planet (1.56 - 2.10 R_Jup) and the cloud top pressure (10^-3 - 10^6 Pa).
	
	In addition to our best-fit model, we fitted for a fully cloudy atmosphere (straight line), and an atmosphere that contains only H_2O as an active gas (no TiO and VO). We calculated the Bayes factors relative to the cloudy model, as follows:
	
	
    B_m = logE_m/E_c

	
	where E_m is the Bayesian evidence of the test model and  E_c is the one of the cloudy model.
	
	We found B_m = 12.2 for the pure-water model and B_m = 12.3 for the model that includes water, TiO and VO. These values correspond to a 5.3 σ detection of water, while they are inconclusive about the presence of TiO and VO <cit.>. 
	
	
	
	
	
	

§ RESULTS

	Figure <ref> reports the stacked spectra obtained with the parametric pipeline and with stripe-ICA. Table <ref> reports the numerical results. The four spectra, i.e. stacked and weighted obtained with the two detrending algorithms, are all consistent within 0.5 σ (maximum discrepancy for a wavelength-bin). Also, the error bars for the corresponding stacked and weighted spectra are similar, within less than 10% in average.
	
	
	
	
	
	
	
	

 §.§ Retrieval Results

	The transmission spectrum of HAT-P-32 b and the best fit to it, retrieved with , is shown in Figure <ref>. The best fitting values and the posterior distributions are shown in Tab <ref> and Figure <ref>.
	With the  exception of water vapor, the fitted values for all the other molecular mixing ratios are smaller than 10^-7. This result means that they are not detectable from this data-set. The water vapor mixing ratio oscillates, instead, between  logH_2O = -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 water vapor mixing ratio, clouds' top pressure, planetary radius at 10 bar and temperature is noticeable in Figure <ref>, indicating there is a degeneracy of solutions.
	
	
	
	
	
	
	
	
	

§ DISCUSSION
 
	

 §.§ Comparison between detrending algorithms
 
	As mentioned in the previous section, the error-bars obtained with ICA are larger by a factor ∼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.
	
	

 §.§ Comparison with other observations
 
	
	Previous ground-based observations of the transit of  in the optical wavelengths <cit.> 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.
	
	

 §.§ Strong water feature
 
	
	Water vapour has been detected, to date, in the atmospheres of about 10 hot Jupiters <cit.>. <cit.> identify two classes of hot Jupiters, essentially mostly cloudy or with a strong water signature. The observed trend suggests that hotter (T_eq >  700 K) and more inflated (logg > 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  (T_eq = 1786 K; logg > 2.8) has one of the strongest water features so far detected (∼ 500 ppm, 5.3 σ).
	
	
	

§ CONCLUSION
 
	We have reported here the analysis of the near-infrared transit spectrum of the hot-Jupiter HAT-P-32 b which was recorded with the  Wide Field Camera 3 on-board the Hubble Space Telescope. 
	
	To obtain the  transit spectrum, we have adopted different analysis methods, which include  parametric and non parametric techniques (Independent Component Analysis, ICA), and compared the results. The final spectra are all consistent within 0.5σ. The uncertainties obtained with ICA are larger than the ones obtained with the parametric method by a factor ∼1.6 - 1.8. The larger uncertainties 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 uncertainties are therefore the worst-case estimates. 
	
	To interpret the spectrum of HAT-P-32 b, we used , a fully Bayesian spectral retrieval code.
	
	
	
	As for other hot-Jupiters, the results are consistent with the presence of water vapor (logH_2O = -4.66_-1.93^+1.66) and probably clouds (top pressure between 5.16 and 1.73 bar). Spectroscopic data over a broader wavelength range will be needed to de-correlate water vapour's mixing ratio from clouds and identify other possible molecular species in HAT-P-32 b atmosphere. 
	
	

§ ACKOWLEDGEMENTS
 
	The authors are supported by the European Research Council (ERC) grant ExoLights. M. D. and T. Z. are also supported by the Istituto Nazionale di AstroFisica - Osservatorio Astronomico di Palermo (INAF-OAPa). G. T. is supported by a Royal Society URF. This work is also supported by STFC (ST/P000282/1).
	
	The authors wish to thank Ingo Waldmann, Marco Rocchetto, Jonathan Tennyson and Sergey Yurchenko for their input.
	
	
	
		apj
	
	
	

§ DILUTION FACTOR FOR THE COMPANION STAR

	
	In this particular data set the nearby companion HAT-P-32B can be separated from the host star HAT-P-32A. This gives us the opportunity to calculate the dilution factor between the two stars , so that it can be used by other studies where the two stars cannot be separated. To calculate the dilution factor we extracted the light-curve of the companion by adjusting the wavelength calibration and using an aperture expanding five pixels above and below the spatially scanned spectrum. Given the undispersed image of the system taken in the beginning of the observation, we found the companion to be shifted by -3.5487 and -23.9709 pixels along the horizontal and vertical axes of the detector, respectively. For this calculation we considered only the spectra obtained during the last HST orbit as we noticed that the dilution factor was varying linearly during the two orbits before the transit. The most possible explanation for this behavior is the long-term, linear, ramp.
	
	
	
	
