© ESO, 2016

1. Introduction

Two decades after the first detection of an exoplanet around a solar-type star by Mayor & Queloz (1995) the field of exoplanet science is fast moving and has expanded into many subfields. A new main focus is the characterization of exoplanet atmospheres. The most successful approach for studying the atmospheric properties of planets has been the measurement of their transmission and emission spectra from multicolor observations of the occultation events in transiting planetary systems. While the emission from a planet can be inferred from the drop in flux during the secondary eclipse, the transmission spectrum of a planet can be obtained during the primary eclipse. This is possible since the atmosphere of planets are opaque at wavelengths, where the atmospheric constituents absorb light causing larger effective planet radii and, thus, deeper transits. Many successful measurements of wavelength-dependent planet radii have been obtained from space using the Hubble Space Telescope (HST; e.g., Charbonneau et al. 2002; Pont et al. 2008; Berta et al. 2012; Sing et al. 2015). Moreover, in the last four years, ground-based observations have also yielded promising results (e.g., Bean et al. 2010; Murgas et al. 2014; Jordán et al. 2013; Gibson et al. 2013a). However, both space-based and ground-based data often are affected by systematic noise signals, which need to be addressed before a high quality transmission spectrum can be extracted. In the past, the correct treatment of these noise signals was subject of scientific debate and led to disagreements between the conclusions of several groups studying the same data sets (e.g., Tinetti et al. 2007; Ehrenreich et al. 2007; Désert et al. 2009; Gibson et al. 2011, for HD 189733b). As a general consequence, this has created doubts concerning the robustness of presented results. In this paper, we aim to demonstrate with the case of HAT-P-32b that reliable results for the transmission spectrum of a planet can be obtained from the ground. The hot Jupiter with a mass of M = 0.860 ± 0.16 M_Jup and a radius of R = 1.789 ± 0.025R_Jup was discovered by Hartman et al. (2011) around a late-type F dwarf star (Vmag = 11.44) at an 2.15 day orbit. The dayside temperature of the planet was measured at T_eq = 2042 ± 50 K by Zhao et al. (2014) from secondary eclipse observations in the H, K_S, 3.6, and 4.5μm bands. An optically close companion was discovered in 2013 by Adams et al. The stellar companion was recently studied in more detail and concluded to be an M-dwarf bound to the HAT-P-32 system from proper motion and AO measurements (Ngo et al. 2015; Zhao et al. 2014). Both studies place the effective temperature of the companion at about T_eff ≈ 3500 K. Following the notation used in these works, in the following we refer to this stellar companion as HAT-P-32B and to the planet host star and the planet as HAT-P-32A and HAT-P-32Ab, respectively. Knutson et al. (2014) observed the HAT-P-32 system among several other planet host stars for radial velocity (RV) trends that could indicate additional companions. They found a long trend signal for HAT-P-32A pointing to the existence of yet another body in the system. A transit timing variation (TTV) study of 45 transit events by Seeliger et al. (2014) looking for evidence of an additional body found no evidence for variations larger than 1.5 min.

Gibson et al. (2013b) obtained a ground-based optical transmission spectrum of HAT-P-32Ab using the Gemini Multi-Object Spectrograph (GMOS; Davies et al. 1997) at Gemini North. Their results show a flat transmission spectrum. We probed a very similar wavelength range using the long-slit method at the OSIRIS instrument at the 10-meter class telescope Gran Telescopio CANARIAS (GTC) aiming to verify the nature of the transmission spectrum and further demonstrate the potential of GTC/OSIRIS as a reliable survey instrument for observations of this kind.

The paper is organized as follows. We first describe the observing setup and data reduction in Sect. 2. This is followed by a description of the white light curve analysis in Sect. 3 and a discussion of the white light curve results in Sect. 4. Here we also address systematic noise signals we found in both data sets. In Sect. 5 we describe the source of the largest noise signal and motivate its correction for the narrowband light curves followed by a description of the extraction of the transmission spectrum during the analysis of the narrowband light curves in Sect 6. We present and discuss our results for the transmission spectrum in Sect. 7 and draw our conclusions in Sect. 8. The study of photometric and spectroscopic data of the companion HAT-P-32B to derive its stellar parameters and extract values necessary for the correction of its diluting effect on the transit depth of HAT-P-32Ab can be found in Appendix A.

2. Observations and data

We observed HAT-P-32Ab twice during transit on 2012 September 15 (JD 2 456 185.5, hereafter referred to as Run 1) and on 2012 September 30 (JD 2 456 200.5, hereafter referred to as Run 2) with the OSIRIS instrument (Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy; Sánchez et al. 2012) mounted at the Spanish 10.4 m telescope GTC. We chose the method of long-slit spectroscopy, in which the planet host star and a suitable reference star are both placed inside one long slit. The grism R1000R was used to disperse the light over the range from 518 to 918 nm. With an exposure time of 10 s (Run 1) and 7 s (Run 2) continuous time series of 321 (Run 1) and 700 (Run 2) optical spectra were obtained each night covering the duration of the whole transit event in both cases. We used slightly different observing setups each night, working with different regions of the CCD detector. The OSIRIS detector consists of a mosaic of two 2048 × 4096 pixel Marconi CCD42-82 chips.

In Run 1 the chosen reference star Ref1 was considerably fainter (ΔVmag = 2.15) than HAT-P-32A and located at a distance of 56.0″ = 0.93′. A custom made 12″ wide slit was used and both stars were placed on CCD1. In Run 2 we chose a brighter reference star, Ref2, (ΔVmag =−0.467) located at a 191.0″ = 3.18′ separation from HAT-P-32A. As a result of the larger distance between the stars, they could not both be placed in CCD1. In order to have both stars on the same CCD and, thus, avoid possible complications from variations in the detector properties we placed both stars on CCD2, which has a larger unvignetted field of view through the slit than CCD1. We also exchanged the custom made 12″ wide slit for a 10″ wide slit, since the latter is slightly longer (extending 0.567′ further into CCD2) giving both stars more room in spatial direction. The setup for both observing runs is illustrated in Fig. 1 and the coordinates for both reference stars are given in Table 1. The observing conditions during both nights were good with an average seeing of 1.06″ in Run 1 and 1.12″ in Run 2. The seeing was not stable in either run, varying between 0.78″ and 1.82″ during the first and between 0.69″ and 2.42″ during the second night.

Owing to complications during the observation almost no out-of-transit data was obtained during Run 1 and part of the data was rendered useless by a light reflection passing over the detector contaminating the red part of the spectrum of HAT-P-32A for approximately 20 min (35 frames) during the second half of the transit. An example of a contaminated frame is shown in Fig. 2.

After each run about 50 bias frames and 100 sky flats were taken. For wavelength calibration, spectra of the HgAr, Xe, and Ne lamps were taken with a 1″ wide slit.

2.1. Data reduction and spectral extraction

We employed standard data reduction procedures, subtracting the median averaged bias and dividing by the median averaged flat field. We then applied a wavelength calibration to every column of the image rebinning the data to a homogeneous wavelength grid using the IDL routine rebinw from the PINTofALE package (Kashyap & Drake 2000), guaranteeing flux conservation. This step ensured that every pixel in an image row corresponded to the same wavelength. The extraction of the stellar spectra from the images was then performed using the optimal extraction algorithm (Horne 1986). The algorithm performed well on our data sets, yielding lower noise levels than other, more simple approaches. During the extraction a median averaged spatial profile was used to identify and mask cosmic ray strikes. Due to the close projected distance of HAT-P-32A and the fainter M-dwarf companion (which was not yet discovered at the time of the observations) it was not possible to reliably exclude HAT-P-32B by choosing a narrow extraction aperture. Instead, we choose a very wide aperture (80 pixel = 20.32′′) ensuring that HAT-P-32B was fully within this aperture at all times. As a consequence, the diluting effect of its additional flux on the transit depth had to be corrected during transit modeling (see Sect. 3).

During the observations the stars drifted slightly in spatial and also in dispersion direction. We monitored the drift in spatial direction by fitting a Gaussian function to the stellar profile tracing the position of the peaks. In this step, we also retrieved the full width at half maximum (FWHM) of the fitted Gaussian profile to monitor the seeing variations. The drift of the stars in dispersion direction caused small shifts in the wavelength solution with respect to that obtained for a star perfectly centered within the slit. We monitored and corrected these wavelength shifts by calculating the cross correlation of each spectrum with the first spectrum of the respective run.

2.2. Light curves

Since any telluric variations during the observing runs will have affected both planet host star and reference star in the same manner, these effects can be neglected if only the relative light curves of these two objects are considered. We created white light curves from the data by dividing the total summed flux over all wavelengths for HAT-P-32A by the total sum of the reference star spectra for each measurement. The resulting white light curves are shown in Fig. 3. We further created 20 narrowband channel light curves from the data of Run 2 by dividing the wavelength range into intervals of 20 nm. The spectra of HAT-P-32A and reference star Ref2 are shown in Fig. 4 together with the channel limits. Channel #13 encompassed both the telluric oxygen bands and the potassium resonance lines (K i at 766.5 and 769.9 nm) predicted for exoplanet atmospheres at moderate temperatures. Akin to what has been reported by Parviainen et al. (2016), we found that the noise level, estimated from the standard deviation of the out-of-transit light curve scatter in each individual wavelength point, rises significantly in the deeper of the two telluric oxygen bands, negatively affecting the signal-to-noise of the entire channel #13 light curve. We therefore constructed an additional channel (channel #13b) of only 15 nm width, spanning from 763 to 788 nm. This channel was still encompassing the expected potassium lines but avoided the high noise region as sketched in Fig. 5.

3. Analysis of white light curves

The relative white light curves for both runs are shown in Fig. 3. We found both white light curves to be affected by red noise. Previous works dealing with data obtained with GTC/OSIRIS have reported on similar noise signals (Sing et al. 2012; Murgas et al. 2014; Mallonn et al. 2015). The three works explored different systematics models that included terms depending on the seeing and air mass. We also found indicators for a possible correlation of the red noise with auxiliary parameters and explored these possible correlations by including a systematic noise model in our light curve fitting. We found the minimal χ² using the IDL implementation of the Levenberg-Markward algorithm mpfit by Markwardt (2009), while modeling the light curves with a model of the form (1)where is the analytical transit model described by Mandel & Agol (2002), is a systematic noise model and c_w is the relative flux contribution of the stellar companion HAT-P-32B in white light (i.e., c_w = f_{HAT - P - 32B}/f_{HAT - P - 32A}, where f is the total flux integrated over the white light wavelength range). In the following we refer to c_w and its equivalents for the narrowband channels as the dilution factor. The Mandel & Agol transit model parameterizes the transit using the radius ratio between the planet and star R_p/R_⋆, the semimajor axis of the planet orbit in units of the stellar radius a_p/R_⋆, the inclination of the planet orbit i, a reference time for the midtransit T_C1 and two coefficients, u₁ and u₂, describing the stellar limb darkening with a quadratic limb darkening law. Also needed is the period of the planet orbit P, which we kept fixed to the value given by Hartman et al. (2011), i.e., P = 2.150008 days, and the eccentricity of the orbit e, which we kept fixed at e = 0. The code implementation of Eastman et al. (2013) was used for the calculation of the transit model. The relative flux contribution of the stellar M-dwarf companion c_w and its uncertainty were determined from our own data of Run 2, as described in detail in Appendix A, and this value was kept fixed during this optimization step. We tested 135 different systematic noise models where each model was a different combination of polynomial functions depending on the four auxiliary parameters: position in spatial direction xpos, position in dispersion direction ypos, seeing/FWHM of the stellar profile fwhm, and air mass airm. All combinations of different polynomial orders between 0 and 2 (between 0 and 4 for the FWHM) for all four auxiliary parameters were explored. The most complex form of we tested, consequently, was of the form (2)where n₀ is the normalization and x_1,2, y_1,2, a_1,2 and f_1,2,3,4 are model parameters scaling the influence of the respective auxiliary parameter.

To avoid overfitting because the minimal χ² generally decreases with a rising number of free parameters, a compromise between the number of free parameters and minimization of the χ² should be found by means of an objective criterion. We used the Bayesian information criterion for the case of an unknown variance of the data (3)where RSS is the sum of squared residuals, k is the number of free parameters, and N is the number of data points, to determine the best choice of . We compared the BIC of all tested models and found that using a combination of a first-order polynomial of the position drift in dispersion direction (ypos), a first-order polynomial of the air mass, and a third-order polynomial function of the FWHM yields the lowest BIC for the data of Run 1, i.e., a model of the form (4)An abbreviated compilation of the results of the BIC comparison can be found in Table 2, where, for practicality, only the models with an ΔBIC = BIC − BIC_min< 10 are listed. The best-fit results for the planet-to-star radius ratio derived using these model approaches differ only slightly, and all lie well within the 1σ uncertainty interval (derived in Sect. 3.1) of the best-fit result obtained with the lowest BIC model. In Table 2 the models are identified by a short code, which is to be read as follows: “xpos” is the key word for the position drift in spatial direction, “ypos” for the drift in dispersion direction, “airm” for the air mass, and “fwhm” for the seeing/FWHM of the stellar profile. The number following each of these key words indicates the highest polynomial order that was allowed to be non-zero and free in the model fit. Using the systematic noise model given in Eq. (4) we achieved a good fit to the data of Run 1 with an almost Gaussian distribution of the residuals (with a normalized standard deviation of 463 ppm). The remaining correlation of the residuals is explored in the Sect. 3.1, which focuses on the error estimates of the results. A plot of the white light curve with the best-fit model can be found in Fig. 6.

After repeating this process for the data of Run 2, we found the minimal BIC for a systematics model consisting of a first-order polynomial of the position drift in spatial direction (xpos), a first-order polynomial of the air mass, and a third-order polynomial function of the FWHM. However, neither this model nor any other tested combination yielded a good fit to the whole curve leaving the residuals still strongly correlated (see Fig. 7). This lead us to conclude that the systematic noise, which is strongly distorting the first part of the light curve of Run 2, has a different origin. We investigate this origin in Sect. 5.

3.1. Errors

For the white light curve fit of Run 1 we probe the probability distributions of the model parameters with a Markov chain Monte Carlo (MCMC) sampling using the emcee package (Foreman-Mackey et al. 2013). We used noninformative priors for all parameters except the dilution coefficient c_w, which we allowed to vary within its uncertainties using a Gaussian prior. The MCMC sampling was run with an ensemble of 600 walkers each starting with slightly different parameter populations as seeds for the chains. We let the chains run for 20 000 accepted steps and defined a burn-in phase, in which the chains are not yet fully converged, to be over after 5000 steps omitting all chain-steps prior to this mark. The results can be considered robust if all the chains converge to the same probability distribution, which they did in our case. We determined a thinning factor for each chain as the largest autocorrelation length of any parameter within the chain. The average thinning factor was ≈200, which resulted in a total of 45 891 accepted sample points for the merged distribution of an accepted chain steps. The correlation plots of the posterior distributions for all parameters of the white light curve fit of Run 1 are shown in Fig. 8. The 1σ uncertainties of each parameter were calculated as the limits encompassing 68.27% of all sampled points.

3.2. Red noise estimation

We calculated the red noise factor β to determine the level of red noise remaining in the light curve residuals and consequent underestimation of the parameter uncertainties. This factor was introduced by Winn et al. (2008) and is based on a comparison of the progression of the standard deviation of the time-binned light curve residuals with the behavior expected for a light curve purely affected by Gaussian (i.e., white) noise. The latter is described by Eq. (5), (5)where σ₁ is the amplitude of the Gaussian white noise, N the number of adjacent points binned together, and M the total number of bins. The difference between theoretical and actual progression of the standard deviation of the binned residuals is caused by red noise. Winn et al. define the ratio between the two curves as the “red noise factor” , which can be used to inflate the MCMC chain derived error bars of all parameters. We determined the maximal value of this ratio to be β_{Run 1} = 1.345 for the white light curve data from the first run when the systematics model with the lowest BIC was used and β_{Run 2} ≥ 6.5 for the data from Run 2, regardless of which systematics model was used. The β factor close to unity for Run 1 demonstrates that the systematics model (Eq. (4)) used was appropriate to correct all noise sources. In contrast, the high value derived for Run 2 clearly indicates a remaining noise source (see Sect. 5), suggesting the model approach used for Run 1 to be insufficient for the analysis of the data of Run 2.

4. Results and discussion of the white light curves

The results for the white light curve fit to the data from Run 1 are listed in Table 3 together with the literature parameters of recent high precision studies of HAT-P-32Ab. The error bars of these results were inflated with the red noise factor of β = 1.345 derived in Sect. 3.1. We find our results for the planet parameters to be consistent with the literature values. The difference in radius ratio compared to the results of Hartman et al. (2011) and Seeliger et al. (2014) might be explained by the different wavelength region probed in this study but is more likely caused by the transit depth dilution from additional flux of the M-dwarf companion HAT-P-32B, which was not accounted for in these two studies.

A rough dilution correction of the given radius ratios using (6)with the dilution factor values , , provided by Zhao et al. (2014) in similar broadband filters as those used in the two studies mentioned above, yields comparable radius ratios to that derived by us and Gibson et al. (2013b).

4.1. Transit timing comparison

In order to compare our best-fit transit time, T_C1,2456185 with literature predictions, we use the ephemeris given in each paper to calculate the expected transit time and its uncertainty using Eq. (7) as follows: (7)where T₀ is a reference transit time, E the epoch, and P the orbital period given in the respective paper. We find that our result for the transit time best agrees with the prediction given by the discovery paper (Hartman et al. 2011). Using their ephemeris information and our measured transit time to refine the period we yield P_new = 2.15000806(24) days.

We also calculated the predicted transit time for Run 2 T_C2,2456200, which we could not reliably measure from a white light curve transit fit. To calculate the prediction for T_C2,2456200 from our Run 1 results, we used Eq. (7) with T₀ = T_C1,2456185, E = 7 and P = P_new. The results are given in Table 3.

4.2. Comparison with theoretical limb darkening

As the limb darkening coefficients are wavelength dependent and cannot be directly compared with literature results obtained in different filters, we compared these coefficients to theoretical values. To derive these theoretical coefficients, we calculated the wavelength-dependent theoretical limb darkening profiles for a star with the basic stellar properties of HAT-P-32A ([Fe/H] =−0.04 ± 0.08, log g (cgs) = 4.33 ± 0.01, Hartman et al. 2011, and T_eff = 6269 ± 64 K, Zhao et al. 2014) via interpolation from the PHOENIX specific intensity spectra library by Husser et al. (2013). We then weighted each wavelength with its actual contribution to the measured stellar flux during observation, taking into account the instrument response function and the telluric absorption and summed all contributing limb darkening profiles to derive the theoretical white light profile. We then renormalized the model information so that μ = 0 actually corresponds to the outer edge of the star, in other words, the region where the mean optical depth corresponds to unity. This is not the case for the raw model data due to the spherical symmetry assumption used in the PHOENIX code. We repeated this process varying the adopted stellar parameters for HAT-P-32A within their reported errors to explore the uncertainties of the intensity profile. Finally, the theoretical limb darkening coefficients and their errors were derived by fitting the intensity profiles with a quadratic limb darkening law. We found that our best-fit limb darkening coefficients were lower than the values predicted by the PHOENIX stellar models (u_1,theory = 0.340 ± 0.056 and u_2,theory = 0.245 ± 0.073). This difference might be caused by insufficiencies of the PHOENIX models, errors in the assumed stellar parameters, or undetected systematic noise in the white light curve of Run 1.

4.3. Consequences for the retrieval of the transmission spectrum of the planet

Despite the good agreement of our Run 1 white light curve results with the literature data, the lack of extensive out-of-transit data led us to deem the data set unsuitable for a transmission spectroscopy study, where a reliable measurement of very small changes in the transit depth is essential.

Further, we considered the results obtained for the white light curve fit of the Run 2 as unreliable since a clear and dominant systematic noise signal in the data remained uncorrected (see Sect. 3.2). We further investigate this noise signal in Sect. 5 and motivate a correction for the narrowband channel data of Run 2, which then can be used to derive a transmission spectrum of HAT-P-32Ab.

5. GTC/OSIRIS instrument specific systematic noise affecting the data of Run 2

In Sect. 3 we determined that the systematics affecting the data of Run 2 cannot be sufficiently corrected by solely using a similarly simple decorrelation function of auxiliary parameters as was sufficient for the data of Run 1. The largest nonastrophysical systematic present in the white light curve of Run 2 is the large distortion of the first half of the light curve. It is present in all narrowband channel light curves (see Fig. 15, left panel). When studying the raw white light curves of Run 2 we found that the distortion in the relative light curve coincides with a sinusoidal-like feature found in both the raw light curve of HAT-P-32A and the reference star (see Fig. 9). This particular and slightly curious feature in the first half of the observation does not appear to be caused by telluric atmosphere variations. A closer investigation of the two bumps showed that while they are present in the data of both stars, they are slightly shifted in phase (and possibly have a different amplitude). As a consequence, they do not cancel out when the ratio of both curves is taken to correct any telluric effects. We found that the bumps reappear with a frequency of roughly 60° of the instrument rotation angle. To better visualize our findings, we divided the raw light curve of HAT-P-32A by a theoretical transit light curve using our best-fit transit parameters from Run 1, cleaned both raw light curves from an air mass trend, and plotted them over the instrument rotation angle (see Fig. 10). We fitted all four bumps with a Gaussian function to derive an estimate for the values of rotation angle at which the peak flux occurs. Since the instrument rotation angle changed during the first half of the observing run and then stayed nearly constant (see Fig. 11) this lead to the bumps only manifesting in the first half of the run. The 60° symmetry of the bump feature prompted us to suspect these flux amplitude variations to be caused by vignetting in pupil space. A vignetting of the pupil would reduce the overall amount of captured light that reaches the detector. Because of the hexagonal shape of the GTC primary mirror, the pupil is also not invariant under rotation but exhibits a 60° rotational symmetry. As a consequence, the amount of vignetted pupil area would not only depend on the distance of the projected pupil but also on its effective rotation angle to the source of the vignetting. Vignetting could, for example, occur at the secondary and/or tertiary mirror for off-axis rays. The GTC conceptual design document (Alvarez et al. 1997, Fig. 4.4) indicates that this should only have a significantly impact on targets located at separations larger than ≈ 4′ from the telescope optical axis (which includes Ref2 but not HAT-P-32A in Run 2; see Fig. 1). Potential vignetting caused by other objects located in the optical path in pupil space might, however, also affect targets located closer to the telescope pointing.

We constructed a simple toy model to simulate the expected footprint of pupil vignetting at a hypothetical source. The effect predicted by this toy model exhibited a 60° periodicity with rotation angle very similar to the signal found in our data of Run 2 (see Fig. 12). The toy model further showed a dependency of the signal on the initial target position relative to the telescope optical axis, which was regulating the phase, amplitude, and precise shape of the effect. The true signal affecting the data could be more complex in case the source of vignetting were located off-center from the telescope optical axis (divergent from the toy model assumption) or if the signal were a superposition of several vignetting signals.

The difference in systematic noise signal amplitude between the two stars of Run 2, which are separated by more than 3′ (≡1/3 of the OSIRIS field of view) appears to be small, suggesting that for any given observing setup no signals with significantly larger amplitudes (i.e., larger than 1%) are to be expected. Systematic noise signals of amplitudes slightly smaller than 1% might be negligible for many other science cases. In our case, however, noise signals of this amplitude are more than two orders larger than the astrophysical signals we are aiming to detect and, therefore, need to be addressed. Since the amount of vignetting is only dependent on the projected pupil position and effective pupil rotation angle, it should be wavelength independent as long as the vignetting occurs before the light has passed any dispersing optical elements. In this case the produced systematic noise signal affects all of the narrowband channel light curves and the white light curve in exactly the same way and should divide out if the ratio between any color light curve and the white light curve is taken.

6. Analysis of the color light curves of Run 2

We proceeded to divide all narrowband light curves of Run 2 by the white light curve. This served to cancel out all wavelength-independent systematic noise signals, that is, signals that affected all curves in an identical manner. Such common-mode corrections are usual practice when dealing with spectrophotometric transit data (Gibson et al. 2013b; Sing et al. 2015). In accordance with expectation, the correction appears to fully dispose of the light curve distortion that we suspect is caused by wavelength-independent pupil vignetting (see Sect. 5). All remaining (i.e., wavelength-dependent) noise signals in the narrowband light curves are found to be linked with auxiliary parameters of the observation in a similar manner as for the white light curve of Run 1 (see Sect. 3). We fitted all differential light curves with a model of the form (8)where is the analytical transit model for the narrowband light curve, c_n the dilution factor (i.e., relative flux contribution of the stellar companion HAT-P-32B) in this particular narrowband, the model for the white light transit, c_w dilution factor in white light, and a systematic noise model for all remaining wavelength-dependent noise sources. The latter was of the same form as the systematic model used during the white light curve fit, that is, a combination of polynomial functions of the auxiliary parameters detector position, air mass, and FWHM.

Alternatively the original relative light curves can be fitted by a model of the form (9)where (10)is the common mode systematic noise signal derived from the residuals of the white light curve data of Run 2. This approach is mathematically identical to that given in Eq. (8).

In both cases all wavelength-independent transit model parameters (i.e., a_p/R_⋆, i, T_C2, P, the white light planet-to-star radius ratio, the white light limb darkening coefficients, and the dilution parameters c_w and c_n), were kept fixed in this optimization. The values for c_n were determined as explained in Appendix A. We chose to fix the white light and wavelength-independent transit parameters to the values we derived from Run 1 in Sect. 3, as they represent a measurement taken with the same instrument and at the same wavelength interval. Not fixing the values but letting them vary in their uncertainty intervals did not change our final results for the relative wavelength-dependent change in radius ratio, and only offset the exact value around which the narrowband radius ratios and limb darkening coefficients varied. The same held true when we fixed the wavelength-independent and white light transit parameters to the literature values provided by Gibson et al. (2013b), Hartman et al. (2011), or Seeliger et al. (2014). Between these tested configurations adopting our best-fit Run 1 white light parameters did yield the lowest overall χ². We, therefore, adopted the values resulting from these white light curves values as our final results. To allow an independent use of our derived transmission spectrum, we also provide the relative change in transit depth with respect to the white light curve transit depth, i.e., , together with our results for the absolute values of the radius ratios in Table 5. These differential transit depths are independent of the white light curve parameters chosen during the analysis and, therefore, free of any possible systematic errors caused by the uncertainties inherent to these parameters.

We again tested all different combinations of systematics models. This time we allowed the highest polynomial order for the air mass and spatial and dispersion position to be 2 and for the FWHM to be 3. The BIC was again calculated using Eq. (3). As we were aiming to determine the model that best explained all 20 channels simultaneously, the joined sum of squared residuals of all 20 light curves (i.e., channel #1–#20, excluding the subchannel #13b) was used in the BIC calculation. The overall lowest BIC was reached with the model (11)where n₀, x₁, y₁, a₁, and f_1,2 are wavelength-dependent parameters and, therefore, different for every channel. An overview of the BIC comparison can be found in Table 4 in which all models with a ΔBIC < 1000 are listed.

When investigating the individual channels separately we found that most of the color light curves were satisfied with a less complex model without a spatial position (xpos) dependent term. The χ² of the three bluest channels, however, was significantly minimized by introducing this additional parameter, leading to an overall favored BIC without affecting the final results for the other channels.

The results from the different models listed in Table 4 are plotted in Fig. 13 together with the values we adopted as our final results. It can be seen that at the blue end of the spectrum a few models give notably different results. The models yielding such deviant results are, however, not those favored when determining the best model for each channel individually. The resulting radius ratios for the individually determined best models (red squares) are in good agreement with the adopted results derived from the homogeneous analysis.

Errors.

We repeated the exploration of the posterior parameter distributions with MCMC as described for the white light curve analysis in Sect. 3.1, now letting the values for the wavelength-dependent dilution factors c_n vary within their uncertainty intervals using Gaussian priors. This time we used 300 walkers and let them run for 10 000 accepted steps. We discarded all points of the burn-in period, which was over after 1000 steps. The typical thinning factors for the chains were ≈82, leaving us with a final distribution of an average of 33 300 points for each wavelength channel light curve fit. An example of the resulting posterior parameter distributions is shown for channel #9 in Fig. 14. We then calculated the red noise factor for every narrowband channel light curve as described in Sect. 3.2 and inflated the MCMC derived error bars accordingly. The resulting values for β lie between 1.022 and 1.950 and the resulting uncertainties of the wavelength-dependent radius ratio after multiplication with β lie between 337 and 972 ppm. The light curves before and after the full systematic noise correction are shown in Fig. 15. The resulting radius ratios and limb darkening coefficients as well as β factors are given in Table 5.

7. Results and discussion of the transmission spectrum

The result for the wavelength-dependent, planet-to-star radius ratio of HAT-P-32Ab is shown in Fig. 16. We find that it shows little variation in the probed region between 518 and 918 nm. The variations are significantly smaller than two planetary atmospheric scale heights (H_p ≈ 1400 km). This estimate for H_P is derived with Eq. (12), (12)using the planet equilibrium temperature T_p = 2042 K given by Zhao et al. (2014), the planetary gravitational acceleration of g_p = 6.6069m/s² from Hartman et al. (2011), and the approximation of the mean molecular weight of the atmosphere m as that of a solar abundance hydrogen helium mixture. In Eq. (12), k_B stands for the Boltzmann constant. The best-fit results for the limb darkening coefficients are shown in Fig. 17. They vary smoothly with wavelength except for the coefficients for channel #13, the channel encompassing the telluric oxygen bands at ≈761 nm. We compare the results for the limb darkening coefficients with theoretical predictions. Our approach to the narrowband channel analysis entailed division by the white light curve. Therefore, effectively only relative changes of the transit depth and stellar limb darkening to the white light transit depth and white light limb darkening were measured. Consequently, the absolute values of our resulting color-dependent limb darkening values are affected by the assumed white light limb darkening values. Since we adopted the best-fit white light curve parameters from Run 1, where the limb darkening coefficients are divergent from theoretical expectations, the resulting color-dependent coefficients are by default divergent as well. To compare these coefficients to the theoretical expectations independently, we calculated modified theoretical values that correspond to the theoretical predictions for the color-dependent coefficients under assumption of the white light coefficients fixed to the best-fit results of Run 1. These modified theoretical values were derived from PHOENIX specific intensity spectra as described in Sect. 4.2, where instead of the whole white light wavelength region, now only the corresponding narrowband channels were summed to derive the limb darkening profile. The resulting modified theoretical prediction is shown together with the actual measured limb darkening coefficients in Fig. 17. It can be seen that the strong change in limb darkening in respect to the neighboring wavelength regions that we see for channel #13 is not expected from theory. Since the results for subchannel #13b do not show the same divergence from theory, we conclude that the results of channel #13 are affected by uncorrected noise due to the strong telluric absorption bands and are unreliable. Channel #13 does not carry any significant information that is not also represented by channel #13b. Therefore, we decide to exclude the results for the planet-to-star radius ratio of channel #13 in the following comparison to literature data and planet atmosphere model predictions.

7.1. Comparison to theoretical models and literature data

We compared the results to theoretical atmosphere models with various degrees of alkali metal abundance depletion (see Fig. 18). The atmospheric models were calculated with a line-by-line radiative code described in Iro et al. (2005) with updated opacities described in Agúndez et al. (2012) and Montalto et al. (2015). We assumed a clear atmosphere (no clouds) in thermochemical equilibrium with a solar abundance of elements for the nominal model. Furthermore, other models with various changes in titanium oxide (TiO) and alkali abundances were explored. The stellar heating was provided using a stellar synthetic spectrum¹ of a G-type star with an effective temperature of 6000 K from Castelli & Kurucz (2003). We found no agreement between the measured transmission spectrum and the expected absorption features of the sodium and potassium resonance doublets (Na i at 589.0 and 589.6 nm, K i at 766.5 and 769.9 nm, where the doublets are unresolved at our resolution) predicted by the nominal model (see top panel of Fig. 18). The radius ratio in channel #4 encompassing the sodium doublet is lower by ≈12.9σ than the predicted value and the radius ratio in channel #13b, encompassing the potassium doublet, is lower by ≈7.4σ. Further, there is no indication of a detection of the predicted broad wings of the sodium feature in the channels neighboring channel #4. Neither an enhancement of the titanium oxide abundance in the probed upper layers of the atmosphere (see top panel of Fig. 18) nor a depletion of either only sodium or only potassium (see middle panel of Fig. 18) improves the fit between the model and data. Reducing the abundance of both these alkali metals by a factor of 1000 yields a model that agrees fairly well with the data with a reduced χ² of 1.58 (20 degrees of freedom). This is a similar but slightly worse match than the match to a straight line representative of a gray atmosphere (see bottom panel of Fig. 18), which results in a reduced χ² of 1.50 (20 degrees of freedom). A gray atmosphere signal could be caused by a high altitude cloud layer masking the fingerprint of the atmosphere below. Alternatively, a significantly lower atmospheric scale height of the terminator region than the one predicted by the measured equilibrium temperature of the planet could be responsible. The amplitude of all expected atmospheric features would shrink accordingly and in an extreme case they would be hidden in the uncertainties of the measurement. To be in very good agreement with the data, the scale height would have to be lower than the current prediction by a factor of 6.8. Such a decrease in scale height could, for example, be caused by either a lower than predicted terminator temperature (≈300 K) or an increase of the assumed mean molecular weight by that factor (i.e., m ≈ 16 kg/kmole; for comparison the mean molecular weight of water vapor is m_H2O = 18.02 kg/kmole). A more plausible explanation, however, would be that both of these two regulating factors are divergent from the values assumed in the current scale height calculation.

Our results for the planet-to-star radius ratio are in good agreement with the study of Gibson et al. (2013b) who observed two transit events of HAT-P-32Ab with GMOS at Gemini North. They analyzed both data sets separately and then combined the final results for the planet-to-star radius ratio of both data sets. In Fig. 19 we compare our results with their combined results.

7.2. Potential of GTC/OSIRIS as a tool for transmission spectroscopy

We found the light curves to be heavily affected by systematics, which can be modeled, however. Only the channel containing the telluric oxygen absorption band at ≈761 nm was too heavily affected by noise to be sufficiently corrected. The data showed correlation with the position of the stars on the chip (and in the slit) with air mass, which is expected because of the difference in color between the planet host star and reference star, and with seeing. The latter correlation suggests slit losses that affect the planet host star and reference star differently. This problem can be overcome by choosing a larger slit width. We adopted this for the following runs (e.g., Parviainen et al. 2016). While the rotation dependent distortion cannot be overcome, it affects the data less when the changes in the rotation angle are small. It also has a smaller impact when the two stars are close to each other on the chip since vignetting should then be similar for both stars. In such cases, that is, when no obvious distortions can be seen in the data, disregarding their possible hidden existence could lead to systematically erroneous transit parameters, but should have negligible impact on the derived transmission spectrum as long as only relative variations in the radius ratio are considered.

8. Conclusion

We were able to derive a high precision transmission spectrum for the inflated hot Jupiter HAT-P-32Ab showing no prominent absorption features and, thus, supporting the results of the earlier work by Gibson et al. (2013b). The study allowed us to detect and understand the low-level instrument systematics affecting GTC/OSIRIS and will help to improve future measurements. The independent confirmation of ground-based results from a different ground-based facility affected by different systematic noise signals is a step toward re-establishing faith in the reliability of (ground-based) transmission spectroscopy measurements. We have confidence in the potential of ground-based facilities and GTC/OSIRIS, in particular, as an excellent tool for larger surveys.

Note added in proof. We would like to note that during a late stage of this paper’s referee process another study of HAT-P-32b by Mallonn & Strassmeier (2016) was uploaded to ArXiv. This study further confirms the measurement of a flat transmission spectrum for this planet’s atmosphere.

---

Acknowledgments

This work is based on observations made with the Gran Telescopio Canarias (GTC) installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias in the island of La Palma. It is partly financed by the Spanish Ministry of Economics and Competitiveness through projects ESP2013-48391-C4-2-R and ESP2014-57495-C2-1-R. L.N. acknowledges support from the DFG Graduiertenkolleg 1351 Extrasolar Planets and their Host Stars. F.M. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under the program ANR-12-BS05-0012 Exo-atmos. Many of the plots shown in this paper were made using Matplotlib (Hunter 2007).

References

Appendix A: Companion HAT-P-32B

In 2013Adams et al. discovered an optical companion to HAT-P-32A using adaptive optics (AO) with Arizona Infrared imager and Echelle Spectrograph (Aries, Principle Investigator: Don McCarthy), which is fed by the Multiple Mirror Telescope’s (MMT) AO beam. The optical companion was later confirmed to be bound to the HAT-P-32 system by Ngo et al. (2015) from proper motion measurements and AO imaging. HAT-P-32B was further characterized as an M-dwarf by Zhao et al. (2014) and Ngo et al. (2015) who used near-infrared broadband AO imaging to constrain its stellar parameters.

Because of the close proximity of the two stars, HAT-P-32A and B, the flux of the latter contributed to our measurements of the former. We needed to determine the wavelength-dependent flux ratio between the two stars as precisely as possible to correctly include this effect in our models. In the following, we describe how we extracted this information from our GTC/OSIRIS spectra. Subsequently, we use this data and additionally obtained near-infrared observations to derive improved stellar parameters for HAT-P-32B.

Appendix A.1: Optical spectrum − GTC/OSIRIS

Since HAT-P-32B was undiscovered prior to our observations the observing setup was not optimized to maximize the projected distance of the two spectra of HAT-P-32A and B on the chip. We still were able to detect HAT-P-32B as a separated object in Run 2 as a deformation of the spatial profile of HAT-P-32A.

Using the out-of-transit data available in Run 2 we determined the flux peak of the HAT-P-32A spectrum in every frame at every wavelength and then added frames in 10-frame time bins resulting in 25 images with increased signal to noise. We then fitted a double-peak model based on an empirical profile function to the spatial profiles of HAT-P-32A and B for every wavelength cut though the spectrum for all 25 images. We used an empirical profile based on a Moffat function and assumed symmetry of the stellar profile. It is described in detail in Sect. A.1.1. A sample fit to the spatial double profile can be found in Fig. A.1.

Appendix A.1.1: Empirical profile

We acquired a first approximation of the peak position and amplitude of the two stellar profiles by fitting each of them with a Moffat function enforcing identical values for the half widths and for the Moffat indices of both profiles. The result showed that the wings of the actual stellar profiles were underestimated by the Moffat approximation. Consequently, we moved on to an empirical approach. For this we assumed that the actual profile is symmetrical with respect to its central axis and that the point spread functions (PSFs) of both stars are the same. We constructed an empirical profile, , using the side of the stellar profile of HAT-P-32A, which is not distorted by the additional flux of HAT-P-32B, and mirroring it at the central axis of the HAT-P-32A profile. We then used a scaled down, shifted in pixel position version of to model the profile of HAT-P-32B. The final model for the superposition of both stellar profiles was of the form (A.1)In this model the scaling factor (s), the center of the stellar profile of HAT-P-32A (x), and the distance of the centers of both stars (x_A−B) were free parameters. From the optimization of this model we obtained the flux for HAT-P-32A (f_A) and HAT-P-32B (f_B) for every wavelength element in each of the 25 images by summing over the respective theoretical profiles.

Appendix A.1.2: Dilution of the optical transit depth measurements of HAT-P-32Ab

Since HAT-P-32B was within the aperture used in Sect. 2.1 to determine the flux of HAT-P-32A, its additional flux had a diluting effect on the transit depth. This effect is color dependent as HAT-P-32B contributes different amounts of flux in the individual wavelength channels defined in Sect. 2.2. We determined the contribution in every wavelength channel by first summing the measured flux of both stars separately within the channel limits in each of the 25 frames. For each channel we then divided the resulting total flux of both stars and adopted the mean of the 25 flux ratios as the final result of the dilution factor c_n = f_A/f_B and their standard deviation as its uncertainty. The resulting values for each channel are listed in Table 5 in Sect. 7.

Appendix A.1.3: Resulting relative optical spectrum

While for the correction of the multicolor channel light curves in this paper only low resolution information for HAT-P-32b was needed, higher resolution information was available. We used this data to determine the stellar properties of HAT-P-32B in Sect. A.4. In order to remove the OSIRIS instrument sensitivity function from the data and clean it of any telluric absorption, we used the relative spectrum of HAT-P-32B and HAT-P-32A for this analysis. Consequently, all wavelength-dependent telluric and instrumental effects are divided out. Since the slit alignment was optimized to center HAT-P-32A and the reference star (Ref2, see Table 1) within it, HAT-P-32B was not perfectly centered in the slit. As a consequence, the wavelength solution for its spectrum is slightly shifted compared to the wavelength solution of HAT-P-32A. We first applied this shift in wavelength and then calculated the flux ratio for every wavelength element in each of the 25 frames. We again adopted the mean of the 25 results as the final value and their standard deviation as the uncertainty. The resulting relative spectrum, covering the wavelength range 518−918 nm, is shown in Fig. A.3.

Appendix A.2: Infrared photometry - WHT/LIRIS

In addition to the optical spectrum we obtained J, H, and K_S band measurements with the William Herschel Telescope (WHT) in the night of 6 October 2012 using the Long-slit Intermediate Resolution Infrared Spectrograph (LIRIS; Acosta-Pulido et al. 2002) in imaging mode. We took 18 images and obtained 15 flats and darks 15 for each filter. A sample image of these observations for each filter is shown in Fig. A.2. We used PSF fitting to obtain the absolute flux for HAT-P-32B and HAT-P-32A. We then calculated the flux ratio f_B/f_A for every image and adopted the mean as the final value and the standard deviation as the uncertainty for each filter. The results are listed in the first column of Table A.1.

Appendix A.3: Comparison to literature

We compared our results of the flux ratio c_filter = f_A/f_B of HAT-P-32B and HAT-P-32A with the broadband measurements from the literature obtained in the optical (g′, r′, i′, z′) by Zhao et al. (2014) with the Robo-AO instrument (Baranec et al. 2014) on the 60 inch telescope at the Palomar Observatory and in the near-infrared (J, H, K′ and Ks) by Zhao et al. (2014) and Ngo et al. (2015) with NIRC2 a near-infrared imager (Principle Investigator: Keith Matthews) using the AO system of Keck-II (Wizinowich et al. 2000) and Adams et al. (2013) using MMT/Aries.

The optical wavelength region probed with our GTC/OSIRIS transit observations fully overlaps with the r′ and i′ bands. We folded our data with the respective filter curves and found the results (c_r′ = 0.0023(7), c_i′ = 0.0064(8)) to be consistent with the study of Zhao et al. (2014, c_Zhao,r′ = 0.003(1), c_Zhao,i′ = 0.006(2).

Since the WHT/LIRIS filter curves for the near-infrared broadband filters J, H, and K_S differ slightly from those used by the Keck-II/NIRC2 facility, an exact comparison between the respective measurements is not feasible. When neglecting these small differences in filter transmission, we found, however, that our results are consistent within 1σ with the near-infrared values derived by all three studies, i.e., Zhao et al. (2014), Ngo et al. (2015), and Adams et al. (2013). Both Zhao et al. (2014) and Ngo et al. (2015) independently analyzed the Keck-II/NIRC2 data obtained in the H and K_S band passes (while the J and K′ band data was only analyzed by Ngo et al. 2015) and arrived at different results, which are only consistent with each other within 2σ.

Appendix A.4: Comparison to theoretical models

Using their broadband measurements, Zhao et al. (2014) and Ngo et al. (2015) both determined physical properties of HAT-P-32B and arrived at similar results with effective temperatures around T_eff = 3550 K. The exact results are listed in Table A.2.

Both studies made use of the PHOENIX stellar atmosphere models but relied on slightly different approaches, stellar parameters for HAT-P-32A, and different fixed values for the metallicity.

We conducted our own analysis by fitting a theoretical model of the form (A.2)to our optical and near-infrared data and the literature broadband measurements, excluding the upper limit Zhao et al. (2014) give for the g′ band and the K_S band data point of Adams et al. (2013) for which no uncertainties were given. In Eq. (A.2) k is a scaling factor (corresponding to the squared radius ratio of the two stars ) and M_{HAT − P − 32B} and M_{HAT − P − 32A} are PHOENIX stellar models, interpolated to specific stellar parameters from the model grid provided by Husser et al. (2013) with trilinear interpolation. During the fit the stellar parameters of HAT-P-32A were allowed to vary within the uncertainties of the given literature values ([Fe/H] =−0.04 ± 0.08, log g (cgs) = 4.33 ± 0.01, Hartman et al. 2011, and T_eff = 6269 ± 64 K, Zhao et al. 2014) and an identical metallicity for both stars was enforced. For the comparison with the broadband points, we folded the PHOENIX model with the respective filter curves. In this step we used the exact instrument specific filter curve (downloaded from the facility web-pages) for each broadband data point. For the optical spectrum, we folded the PHOENIX model spectrum with a Gauss function to reduce the resolution to match the data.

We found that if we only fit the broadband measurements we arrive at similar values for the stellar parameters of HAT-P-32B as those derived by Zhao et al. (2014) and Ngo et al. (2015). When including the optical data, however, the fit clearly favored cooler temperatures. Owing to the inconsistencies in the literature regarding the results obtained for the H and K_S bands with Keck-II/NIRC2, we decided to include only our own data in the final optimization. The difference in the best-fit effective temperature due to this exclusion of broadband points is insignificant (6 K). Our results for the stellar parameters of HAT-P-32B are given in Table A.2. The uncertainties of these results were inflated with the red noise factor β, which was calculated as described in Sect. 3.2. The errors do not incorporate any intrinsic uncertainties of the PHOENIX stellar atmosphere models.

The best-fit model relative spectrum is plotted together with all data points (including those that were not regarded in the fit) in Figs. A.3 and A.4.

It stands out that the K′ band result of Ngo et al. (2015) diverges significantly from all K_S band values including their own. This offset is too large to be explained by the difference between the K′ and K_S band passes. Since the K_S and K′ observations by Ngo et al. (2015) were conducted at different dates (13 month apart), a more likely explanation for this inconsistency could be stellar variability of HAT-P-32B due to activity, i.e., star spots or flares. We measure a prominent H_α emission line at 656 nm in the optical spectrum (see Fig. A.3), which is indicative of such stellar activity.

All Tables

All Figures