© ESO 2018

1. Introduction

In the past few years, spectroscopy at resolving powers above 20 000 has become a valuable tool to unravel the composition, structure, and dynamics of exoplanet atmospheres. The key aspect of this technique is its ability to resolve molecular bands into the individual lines. This enables line-by-line comparison of model spectra to observations, typically obtained by cross-correlation. Hundreds or thousands of molecular lines are co-added into a single cross-correlation function (CCF), enhancing their faint signal by a factor of approximately $\sqrt{N_{\mathrm{lines}}}$. High spectral resolution also allows us to resolve the change in the orbital radial velocity of the planet. For exoplanets orbiting their parent star in just a few days, these changes are in the order of tens of km s⁻¹ across a few hours. Thanks to its unique and highly variable Doppler signature, the spectrum of the planet can be disentangled from the main contaminants (i.e. the spectrum of the Earth’s atmosphere and the stellar spectrum), which are instead stationary or quasi-stationary. Thus, detections at high spectral resolution are particularly robust against contamination from spurious sources. Furthermore, they are self-calibrated, eliminating the necessity of a reference star of similar brightness and spectral type in the same field of view.

After more than a decade of inconclusive attempts at using high-resolution spectroscopy for investigating exoplanet atmospheres (Collier Cameron et al. 1999; Charbonneau et al. 1999; Wiedemann et al. 2001; Leigh et al. 2003; Lucas et al. 2009; Barnes et al. 2010; Rodler et al. 2010), Snellen et al. (2010) devised a successful observational strategy and detected CO absorption lines in the transmission spectrum of exoplanet HD 209458 b via near-infrared (NIR) spectroscopy with the Cryogenic InfraRed Echelle Spectrograph (CRIRES, Kaeufl et al. 2004) at the ESO Very Large Telescope (VLT). Two years later, Brogi et al. (2012) and Rodler et al. (2012) demonstrated that the same measurements could be extended to non-transiting planets and used to estimate their orbital inclination and hence true mass. While CRIRES began to routinely detect CO and H₂O in the atmospheres of transiting (de Kok et al. 2013; Birkby et al. 2013) and non-transiting (Brogi et al. 2013, 2014; Birkby et al. 2017) planets, Lockwood et al. (2014) showed that similar measurements were possible even at a lower resolving power of 25 000, obtained with NIRSPEC at Keck. They designed an alternative analysis where the planet’s Doppler shift no longer needs to change during an observing night (Piskorz et al. 2016, 2017), potentially enabling the study of longer-period bodies. Besides detecting molecular species, determining their relative abundances, and constraining the overall thermal structure of the atmosphere, high-resolution spectroscopy is also capable of measuring planet rotation and even global winds (Schwarz et al. 2016; Brogi et al. 2016) thanks to the broadening and asymmetry in the cross-correlation function.

So far the major downside of the method has been the enormous request in terms of signal-to-noise. Only the biggest telescopes on the ground, namely VLT and Keck, have succeeded in providing an adequate collective area for securing solid detections in the near-infrared, and only for the brightest known exoplanet systems (stellar K magnitudes below 7.5). This has limited the sample to a handful of hot Jupiters orbiting very bright stars with spectral types ranging from late F to early K. This limitation is, however, tied to the technology of near-infrared spectrographs that were available in the past, and can be overcome by the recent development of instruments with bigger throughput and/or larger spectral range. The high resolution spectrograph GIANO (Oliva et al. 2006; Origlia et al. 2014), mounted at the Nasmyth-A focus of the 3.6-m Telescopio Nazionale Galileo (TNG), belongs in the latter category. With its nearly contiguous coverage of the entire Y, J, H, and K bands, it provides a 21-fold increase in spectral range compared to CRIRES at the VLT.

In this work we have successfully tested the capabilities of GIANO for characterising exoplanet atmospheres. We analysed spectra of HD 189 733 b (Bouchy et al. 2005) taken in and around transit and extract the signal of the transmission spectrum of the planet. HD 189 733 b is considered one of the beststudied exoplanets to date, due to the brightness of its parent star (K = 5.54). Therefore it serves as the ideal testbed for the scalability of high-resolution spectroscopy to smaller telescope facilities.

2. The GIANO spectra and their calibration

Two transits of HD 189 733 b were observed on July 11 and July 30, 2015, with GIANO at the TNG. The data consist in a sequence of nodded observations following an ABAB pattern, each exposed for 60 s. At each nodding position, a fibre of 1″ in diameter (on-sky) feeds the spectrograph. A second fibre, located at the other nodding position (3″ away) is instead looking at the sky, providing an accurate reference for subtracting the thermal background and telluric emission lines. Each time the telescope nods, the fibres swap, i.e. the object fibre becomes the sky fibre and vice versa.

The fibres were re-imaged onto a 2× slicer, resulting in two spectral tracks per nodding position, or four tracks per order. The entire NIR spectrum (0.95–2.45 μm) was cross-dispersed at a resolving power of R = 50 000 and orders 32–80 are imaged on a 2k × 2k HAWAII-2 detector. Contiguous wavelength coverage is ensured for orders higher than 45 (bluer than ~1.7 μm), whereas for lower orders (redder wavelengths) the coverage is incomplete, down to 75% for order 32. In this work, we have renumbered the orders so that order 0 is the reddest and order 46 the bluest. We skipped the extraction of the two bluest orders in the cross-dispersed spectrum for lack of photons.

Both observations consist of spectra taken before, during, and after the transit of the planet (Fig. 1). We collected a total of 90 and 110 spectra in the night of July 11 and 30, respectively. We measured a signal-to-noise ratio (S/N) of 50–60 per spectrum per pixel averaged across the entire spectral range. Peak values of S/N = 70−80 were reached in the H band and correspond to particularly transparent regions of the Earth’s atmosphere.

Fig. 1. Orbital phase of exoplanet HD 189 733 b and corresponding airmass during the two nights of observations analysed in this work. Plus and star symbols denote out-of-transit and in-transit data, respectively.

2.1. Extracting the one-dimensional spectra

Each science spectrum is processed with our implementation of the IDL algorithm fixpix¹ with 10 iterations to correct for bad pixels. Difference images are created by subtracting each B- image (i.e. the spectrum taken in nodding position B) from the A-image, after dividing each image through the master flat field. This subtraction removes the thermal background and sky emission lines.

Differently from common practice, we did not merge the A and B spectra into a combined AB spectrum. Conversely, we extracted the spectra separately at the A and B positions to increase the time resolution of our observations. This choice also overcomes the necessity of measuring the shift between the dispersion solution in A and B spectra at this stage. This step would be redundant since we re-align the spectra to the telluric reference frame in Sect. 2.2.

We used the master flat-field to identify the apertures corresponding to each order, each nodding position, and each sliced spectral trace. On the (A–B) flat-fielded science images, we refined the position of the apertures by fitting the spectral trails with four Gaussian profiles per order simultaneously (two sliced images and two nodding positions) at 32 equally spaced positions along the dispersion direction. These sampled values were then fitted with a second-order polynomial to the full 2048 pixels of the dispersion direction and for each pixel the spectral traces are rectified in the spatial direction via spline interpolation. Lastly, the extraction of the one-dimensional spectra is performed by optimal extraction (Horne 1986) on the rectified images. B spectra need to be multiplied by −1 to correct for the initial subtraction of the B images. The count levels of each spectrum are also rescaled to that obtained with a rectangular extraction. This is useful for S/N estimates, see Sect. 5.

2.2. Aligning the spectra to a common reference frame

The high spectral resolution of GIANO allows us to resolve the orbital motion of the planet HD 189 733 b during transit. Its radial velocity indeed changes by approximately 32 km s⁻¹ from ingress to egress. In contrast, the absorption spectrum of the Earth’s atmosphere (telluric spectrum) is stationary in wavelength. We intend to exploit this intrinsic difference in the Doppler signature of the planetary and telluric signals to disentangle the former from the latter. As part of this strategy, we align all the observed spectra to the telluric reference frame.

We measured the position of a set of telluric lines in every order and every spectrum in the sequence. The left-hand panels in Fig. 2 show that there is a clear trend in wavelength between the dispersion solution for the A positions (in red) and that of the B positions (in blue). This trend seems consistent between the two observing nights (top and bottom row). The right-hand panels show the temporal variations in the wavelength solution once the overall trend is removed. Not only does the wavelength recorded on a certain pixel of the detector vary as a function of time, but also the amplitude of this variation differs between the two observing nights. Whereas on July 11 the dispersion solution of GIANO is stable to 1/20 of a pixel (1 pixel equals to 2.7 km s⁻¹), on July 30 it varies by at least 0.2 pixels. We further discuss the stability of the night of July 30 in Sect. 4.2. For both nights, we proceeded to re-align each spectrum in the temporal sequence by spline interpolating it based on the measured shifts.

Fig. 2. Stability of the GIANO spectrograph. The dispersion solution of the spectra taken at A and B position differs as a function of the spectral order considered (left panel, averaged over the observing night). In addition, we measure temporal variations in the global dispersion solution (right panel). Their amplitude differs between the first and the second night of observations (top-right and bottom-right panels, respectively). These measurements suggest that although A and B spectra move coherently (i.e. their relative shift stays the same), their absolute position is variable on a timescale of minutes. Furthermore, the overall stability of GIANO differs between observing nights.

2.3. Calibrating the spectra in wavelength

Once the spectra are all in the telluric reference frame, we calibrated their wavelengths. We note that the set of arc frames taken after the observing night are not suitable for calibration at the level of precision required by our analysis. Since we showed in Sect. 2.2 the instability of the dispersion solution in time, any calibrations non-simultaneous to the science observations could result in a biased wavelength solution. In the NIR telluric absorption lines provide instead an excellent simultaneous calibration source, except for some orders of the J band where the Earth’s atmosphere is particularly transparent. However, the star HD 189 733 is a K1-2V dwarf and it shows significant absorption lines at near-infrared wavelengths. These can also be used for calibration.

We built a template for wavelength calibration by multiplying a Phoenix stellar spectrum (Allard et al. 2011) matching the properties of HD 189 733 by the atmospheric transmission spectrum generated via the ESO Sky Model Calculator². The stellar spectrum is Doppler-shifted based on the barycentric and systemic velocities at the time of the observations. For each order, we visually paired absorption lines in the template and in the averaged observed spectrum. We then computed the wavelengths of the line centroids in the template spectrum, and the pixel value of the corresponding lines in the averaged spectrum. The centroids were determined by computing a super-sampled version of each line via spline interpolation. Super-sampling is achieved via spline interpolation at the 1/20 of a pixel. We recorded the fractional pixel at which the flux in the super-sampled observed line reaches a minimum, and the corresponding wavelength in the super-sampled template spectrum. We exclude from the analysis lines that appear saturated, close to saturation, or blended. The (pixel, wavelength) relation obtained from the above is fitted with a fourth-order polynomial, and the fit is assumed as the wavelength solution. For most of the orders we achieve a residual scatter per line well below 1 km s⁻¹. Figure 3 shows the residuals as a function of wavelength, and its +1σ (0.75 km s⁻¹) and −1σ (1.25 km s⁻¹) dispersion on the global sample. Some of the orders corresponding to spectral windows where the Earth’s atmosphere is particularly opaque do not have enough flux or unsaturated spectral lines to calibrate in wavelength. Furthermore, for a few orders in the J band, we notice an evident mismatch between modelled and observed spectral lines, which also prevents us from obtaining a reliable calibration. Due to these challenges, we exclude six orders from the night of July 11 (orders 8, 9, 10, 23, 24, 25, corresponding to the ranges 983.06–983.30 nm, 1352.7–1408.7 nm, and 1804.8–1937.8 nm) and six orders from the night of July 30 (orders 8, 23, 24, 43, 45, 46, corresponding to 970.8–983.4 nm, 1009.7–1032.3 nm, 1352.8–1408.7 nm, and 1895.1–1937.8 nm).

Fig. 3. Residuals of the wavelength calibration (in km s−1) vs. wavelength for the full spectral coverage of GIANO. The dashed lines denote the ±1σ in the distribution of the residuals. For most of the orders, a sub-km s−1 per line wavelength solution is achievable.

3. Extracting the planet signal

At this stage of the analysis the planet signal is outshone by orders of magnitude by the stellar and telluric spectra. However, the planet’s orbital velocity has a non-zero radial component during transit, which amounts to −16 km s⁻¹ in ingress and +16 km s⁻¹ in egress. Consequently, while telluric and stellar lines are stationary or quasi-stationary (the stellar barycentric velocity changes by 0.2 km s⁻¹ during transit) in wavelength, the planet spectrum experiences a detectable change in Doppler shift during the 110 min of transit. The analysis described in this section aims to make use of this peculiar Doppler signature to remove the telluric and stellar spectra while preserving the planet signal as much as possible.

3.1. Removing the Earth’s absorption spectrum

We have treated each order separately, as a bidimensional matrix in which wavelength (λ) is on the x-axis and time (t) on the y-axis. Due to the possible differences in the properties of the detector (quantum efficiency and bad pixels) between the two nodding positions, A and B spectra (i.e. odd and even spectra in the temporal sequence) are also treated separately. Lastly, we report a clear discontinuity between the right and left quadrant of the detector. Hence, we also separated x-pixels 0–1023 from 1024–2047. We therefore process 41 × 2 × 2 = 164 matrices of spectra for each observing night. We note here that we already discarded the six orders with inaccurate wavelength calibration (Sect. 2.3).

The analysis begins by taking the natural logarithm of each matrix. This choice reflects our expectation that at first order the depth of telluric lines depends on the exponential of the geometric airmass (a(t), see below). Before taking the logarithm, all the pixels below a linear count level of 0.1 are masked to avoid divergent or infinite values. Residual bad columns of data (i.e. certain wavelengths falling on damaged portions of the detector), previously identified by visual inspection, are also masked.

With the data in logarithmic space, the median of the 300 brightest pixels is subtracted from each observed spectrum to remove variations in throughput due to pointing, seeing, and sky transparency. A temporal average of each matrix is then taken to construct the mean observed spectrum for the night. A robust linear fit between each observed spectrum and the mean spectrum is computed and subtracted out, removing most of the spectral features. However, an additional correction is required to model the change with time of the depth of telluric lines. For each λ_i in the data, this is done by subtracting the fit of the observed flux log [F(λ_i, t)] with the geometric airmass a(t). If all the species in the Earth’s atmosphere were vertically mixed and had constant abundances, a linear fit with airmass would suffice. However, the water vapour content is likely to change during the night, and in addition water is confined to the Earth’s troposphere. Consequently, the amount of water vapour along the line of sight will change in a non-trivial fashion. In these data, we fitted the functional dependence by airmass as$$\log \left[F({\lambda }_i,t)\right]=c_{i0}+c_{i1}a\left(t\right)+c_{i2}{a}^2\left(t\right),$$(1)

which we find far superior to a linear fit in modelling the depth of telluric lines. The fit is computed and subtracted out for each of the spectral channels in the data (each of the columns in the (λ, t) matrix). Although we still detected residual telluric signal in our data (see Sect. 3.2), we verify that it is not possible to improve the correction by just increasing the degree of the polynomial in a(t). We suspect that this happens because even on a photometric night the water vapour content changes irregularly on timescales of minutes, and these variations cannot be easily modelled with a smooth function of time.

With telluric lines removed as above, the previous mask for bad columns and/or low-count pixels is applied again, setting the corresponding pixels to zero. An additional sigma-clipping is applied with threshold of eight times the standard deviation of each matrix to ensure that no residual strong outlier is present. Lastly, the exponential of the data is taken to restore linear flux values.

Since we removed only the time-averaged quantities from the data, or we modelled each spectral channel as function of time separately, we preserved most of the planet signal in the process. This is due to the fact that absorption lines in the planet spectrum will shift wavelength during the night, and therefore will fall on a certain spectral channel only for a limited amount of time.

3.2. Enhancing the planet signal via cross-correlation

Even if preserved, the planet spectrum is still below the noise threshold at this stage. The deepest absorption lines in the planet transmission spectrum barely reach a depth of 1% compared to the stellar continuum, while our typical S/N is 50–60 per spectrum per spectral channel. However, there are thousands of molecular lines in the spectral range covered by GIANO. Their signal can be co-added by cross-correlating the residual data after Sect. 3.1 with template spectra for the planet atmosphere. These are computed via plane-parallel, line-by-line radiative transfer calculations following the prescriptions of Brogi et al. (2016) for what concerns T –p profile and H₂– H₂ collision-induced absorption. In this work CO is not included as trace gas. Given that this molecule mainly absorbs at selected wavelength intervals (longer than 2290 nm, and much more weakly in the range 1560–1600 nm), that GIANO data is affected by strong modal noise (i.e. noise related to fibre modes) in the K band, and that the stellar Rossiter-McLaughlin effect would likely dominate the CO signal without an appropriate correction of the stellar spectrum (Brogi et al. 2016), we opt to explore H₂O and CH₄ as trace gases in this analysis, with opacities taken from HITEMP2010 (Rothman et al. 2010) and HITRAN2012 (Rothman et al. 2013) respectively. Not only are these two species strong tracers of chemistry in hot-Jupiter atmospheres (Moses et al. 2011, 2013; Madhusudhan 2012), but they also show a dense forest of spectral lines across most of the NIR. Hence, they are well suited to exploit the large spectral range of GIANO.

In our calculations we first set the volume mixing ratio (VMR) of methane to 10⁻¹² and explore VMRs of 10⁻⁵, 10⁻⁴, and 10⁻³ for water. We then set VMR(H₂O) = 8 × 10⁻⁴, i.e. the chemical equilibrium value for T = 1500 K, solar C/O, and solar metallicity (Madhusudhan 2012), and add CH₄ at VMRs between 10⁻⁷ and 10⁻⁴, in increments of one order of magnitude. We kept the T –p profile in the lower part of the atmosphere (p > 0.1 bar) fixed and constrained to results from past theoretical and observational studies (Madhusudhan & Seager 2009). However, we have explored two different temperatures (500 K and 1500 K) for the upper part of the atmosphere. Although these would produce completely different emission spectra (i.e. absorption and emission lines, respectively), they have a much smaller impact on the transmission spectrum. The main effect is on the overall depth of absorption lines, with the hotter upper atmosphere producing models with up to 4× the contrast ratio of the colder models. The influence on the relative strength of absorption lines is instead minimal. In conclusion, for this study we test a set of 14 models (seven relative abundances and two T –p profiles).

Each model has the average transit depth subtracted before cross-correlation. In this way both the models and the data have their continuum set to zero. The cross correlation is computed on a fixed grid of radial velocity lags between −225 and +225 km s⁻¹, in steps of 2.7 km s⁻¹, which approximately matches the pixel scale for most of the GIANO orders. This results in a lag vector of 167 radial velocity values. For each spectrum, each order, and each radial-velocity lag, the Doppler-shifted model spectrum is computed by spline interpolation and cross-correlated with the data. The output of the cross-correlation is a set of 41 matrices (one per order), each with dimensions (167 × n_spec), where n_spec is the number of spectra on each night (90 and 110 on July 11 and 30, respectively).

Even at this stage the planet signal is not expected to be detectable yet. We therefore proceeded to co-add the 47 cross-correlation matrices with equal weights, even though the signal in each order will depend on (i) the transparency of the Earth’s atmosphere, (ii) the density and depth of absorption lines in the spectrum of HD 189 733 b, and (iii) the position of exoplanetary lines compared to strong telluric lines. Factors (ii) and (iii) could be computed from our set of model spectra used for cross-correlation, but this would make the weighting modeldependent. In order to avoid this bias, we therefore assumed equal weighting.

The remainder of this Section describes the analysis tuned specifically for the night of July 11, 2015. This is due to the fact that on the night of July 30, 2015, no signal is detected (see Sect. 4 and more specifically Sect. 4.2). However, the qualitative behaviour of the cross-correlation signal and the consequent strategy, discussed below, apply to both nights. As illustrated in the left panel of Fig. 4, the co-added cross-correlation signal shows a prominent feature spanning the entire sequence and smoothly shifting in time from anti-correlation (lighter colours in the plot) to correlation (darker colours). We identified this feature with spurious residual telluric absorption because it is centred at zero radial-velocity lag in the telluric frame, it appears regardless of the planet orbital phase (hence it is not correlated with the transit of the planet), and it is enhanced when cross-correlating the data with a model for the Earth’s transmission spectrum (not shown here), containing the cold H₂O spectrum rather than the more complex and rich hot H₂O spectrum of the planet HD 189 733 b.

Fig. 4. Cross-correlation between the sequence of spectra of HD 189733 observed with GIANO on July 11, 2015, and a model spectrum for the planet HD 189 733 b containing H2O at VMR = 10−4, shown as function of radial velocity and planet orbital phase. Ingress and egress of HD 189 733 b are marked with horizontal dashed lines. In the left panel, we co-added all the GIANO orders after removing telluric lines similarly to past CRIRES observations (Brogi et al. 2014, 2016). In the right panel, we applied an additional mask at the position of telluric H2O lines. The latter step removes the majority of time-correlated residual telluric noise, leaving a faint trail of positive correlation along the expected planet radial velocity (marked by slanted dashed lines). Details can be found in Sect. 3.1.

In past CRIRES data where we de-trended telluric lines by airmass (Brogi et al. 2012, 2013, 2014), residuals were clearly visible in the processed spectra prior to cross-correlation (i.e. at the end of the analysis in Sect. 3.1). Therefore, the strongest measured residuals at the position of telluric lines could be used to de-trend the rest of the spectral channels, generally achieving a near photon noise correction. In these GIANO data, however, no residual is clearly visible prior to cross-correlation. We therefore took a more straightforward approach where we mask the data prior to cross-correlation by assigning zero value to those spectral channels corresponding to telluric lines deeper than 5–30%, with the actual threshold varying based on the order but not significantly influencing the outcome. We repeated the cross correlation, and visually inspect the matrices for each order. For 24 of the 41 orders with good wavelength calibration, any residual telluric signal is suppressed below the cross-correlation noise. For the remaining 17 orders, we still identified some cross-correlation noise by visual inspection, although in broader patterns spanning tens of km s⁻¹ in radial velocity. These noise structures do not resemble the signature of residual telluric lines. We therefore also discarded these additional orders from the analysis, ending up utilising approximately 50% of the available spectral range of GIANO (orders 2, 7, 12–17, 19, 27–30, 33–35, 38–44, and 46). We note that some of the discarded orders are affected by either strong telluric absorption (orders 8–10, 22–25) or modal noise (orders 0–6). They are both likely causes of this extra cross-correlation noise.

With the 24 good orders selected as above, we achieved a co-added matrix of CCFs relatively free from telluric contamination (Fig. 4, right panel). A slanted trail of darker cross-correlation values (higher correlation in the adopted colour scheme) is barely visible at the expected planet radial velocity (along the slanted guidelines) and only in-transit. This is the CCF of the planet transmission spectrum, which we then proceeded to co-add in the rest-frame of HD 189 733 b to obtain the total signal from the planet.

Co-adding requires computing the planet radial velocity ν _P in the telluric reference system. Assuming a circular orbit, this is given by:$${\nu }_{\mathrm{P}}(t,K_{\mathrm{P}})={\nu }_{\mathrm{bary}}\left(t\right)+{\nu }_{\mathrm{sys}}+K_{\mathrm{P}} \sin \left[2\pi \phi \left(t\right)\right],$$(2)

where ν_bary is the velocity of the barycentre of the solar system with respect to the observer, ν_sys is the systemic velocity of HD 189 733 (−2.3 km s⁻¹, Triaud et al. 2009), K_P is the planet orbital radial velocity semi-amplitude, and φ the planet orbital phase. The latter is zero at the centre of the transit, 0.5 at mid-secondary eclipse, and 0.25–0.75 when the planet is in quadrature. It is computed from the time t of the observations, the time of mid-transit T₀ and the planet orbital period P (Agol et al. 2010) as the fractional part of$$\phi \left(t\right)=\frac{t-T_0}{P}\mathrm{.}$$(3)

For each value of K _P we shifted each cross correlation function at time t via linear interpolation so that it is centred around v_P. We then co-added the shifted cross-correlations in time to obtain the total cross-correlation signal from the planet as a function of rest-frame velocity ν_rest and planet radial velocity semi-amplitude K_P (Fig. 5, left panels). If the planet transmission spectrum is detected, we measure a signal at ${\nu }_{\mathrm{rest}}={0.00}_{-0.16}^{+0.19} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$ and $K_{\mathrm{P}}={152.5}_{-1.8}^{+1.3} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$. The latter being computed from the semi-major axis, orbital period, and orbital inclination in the literature (Triaud et al. 2009; Agol et al. 2010), and by error propagation. Since the uncertainty in ν_rest varies with orbital phase via Eq. (2), we computed it for every value of φ and adopted the mean value of the series. This simulates the shifting and co-adding of the CCF illustrated above.

Fig. 5. Detection of H2O in the transmission spectrum of HD 189 733 b observed on July 11, 2015 (top panels) and the non detection on the night of July 30, 2015 (bottom panels). The total cross-correlation signal is shown as function of rest-frame velocity νrest and planet radial velocity semi-amplitude KP. The measured S/N and corresponding significances are shown in the left and right panels, respectively. These are computed as explained in Sect. 3.3. The cross-correlation shown here is obtained with a model spectrum containing H2O at VMR = 10−4 and CH4 at VMR = 10−12.

We note that although K _P is generally well known for transiting exoplanets, we still explore a full range of values. This allows us to verify that no other spurious signals produce a significant detection near the planet’s position. In Fig. 5 (top panels), for instance, residual uncorrected signal is still present at very low K _P. This is typically due to either stellar residuals (if at zero rest-frame velocity) or telluric residuals (at non-zero rest-frame velocity, as in this case). Exploring a sufficiently large parameter space offers a strong diagnostic power on all sources of noise and we strongly recommend its implementation every time the signal from the planet can be separated in Doppler space from stellar and telluric lines.

3.3. Determining the significance of the signal

Following Brogi et al. (2012, 2013, 2014), we estimated the significance of the detection with two different methods. Firstly, we divided the peak value of the cross-correlation function at each K_P by the standard deviation of the noise away from the peak. In this way we computed the S/N of the detection, which is a good proxy for its significance in the absence of significant correlated noise. The latter hypothesis is verified by studying the distribution of the cross-correlation values more than 10 km s⁻¹ away from the planet radial velocity (hereafter denoted as “out-of-trail” values). In Fig. 6, we compare it to a Gaussian curve with the same mean and standard deviation as the sample. We did not measure any deviations from a Gaussian distribution down to approximately 4σ, the limit being imposed by the number of available cross-correlation values. In the same figure, we show in red the distribution of the cross-correlation values within 4 km s⁻¹ from the planet radial velocity (“in-trail” values). It is already evident even from visual inspection that while for the first night (top panel) the in-trail distribution has a higher mean that the out-of-trail distribution (consistent with a higher degree of correlation), in the second night (bottom panel) the two distributions broadly overlap.

Fig. 6. Distribution of cross-cross correlation values away from the radial velocity of HD 189 733 b (out-of-trail values, in blue), compared to a Gaussian distribution with the same sample mean and variance (dashed line, light blue). The top panel is for the first night of observations, the bottom panel for the second. In both datasets, the sample distribution shows no deviation down to approximately four times the standard deviation. For comparison, the corresponding distribution of the cross-correlation values around the planet radial velocity (in-trail) is shown in red. In the data from the first night, the two distributions appear shifted as one would expect if the planet transmission spectrum is detected. Conversely, no evident shift is detected in the data from the second night. The statistical significance of these shifts is quantified through a generalised t-test as explained in Sect. 3.3.

We quantified the above statement statistically by performing a generalised t-test on the data. The null hypothesis H₀ is that out-of-trail and in-trail values have the same mean. The test is performed for the same range of ν_rest and K_P as before. For each value of the two parameters, different sets of cross-correlation values populate the in-trail or out-of-trail sample, and we translated the corresponding t value into a significance at which H₀ is rejected (Fig. 5, right panels), which is by definition the significance of our detection.

4. Results

Despite the exclusion of a fraction of the data due to telluric residuals, we clearly detect the absorption of water vapour in the transmission spectrum of HD 189 733 b taken on July 11, 2015. When cross correlating with a model containing H₂O at VMR = 10⁻⁴ and CH₄ at VMR = 10⁻¹², we measure a S/N = 4.9 and a significance of 5.5σ (Fig. 5, top panels). Varying the VMR of H₂O, as well as including non-negligible VMRs of CH₄ in our models result in a marginal decrease in S/N for the entire range of VMRs tested. This suggests that CH₄ is not detected in the planet’s transmission spectrum, in line with previous measurements in the K band (Brogi et al. 2016). We do not detect the transmission spectrum of HD 189 733 b on the second available transit (July 30, see Fig. 5, bottom panels). In Sect. 4.2 we further discuss this non-detection.

By averaging the maxima in the CCF and t-test matrices (Fig. 5, top panels), we derived a planet maximum radial velocity of $K_{\mathrm{P}}={154}_{-40}^{+36} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$. The 1σ uncertainty corresponds to a drop of one in either S/N or significance. The measured K_P is fully compatible with the literature value reported in Sect. 3.2 ($K_{\mathrm{P}}={152.5}_{-1.8}^{+1.3} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$), and also with the value measured in Brogi et al. (2016), that is $K_{\mathrm{P}}={194}_{-41}^{+19} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$. We note that the large uncertainty in K_P is due to the relatively small change in planet redial velocity during transit (a few tens of km s⁻¹), insufficient to tightly constrain the orbit at this level of S/N. Past dayside observations, particularly those co-adding spectra spanning a wide range orbital phases (Brogi et al. 2012, 2014) provided instead much tighter constraints on K_P, with error bars of a few km s⁻¹. We measured a net blue-shift of the planet CCF ($-{1.6}_{-2.7}^{+3.2} \mathrm{km} {\mathrm{s}}^{\mathrm{-}\mathrm{1}}$). Although very marginal, this is in line with previous measurements by Brogi et al. (2016) and Louden & Wheatley (2015) regarding the strength of day-to-night side winds at the planet’s terminator. We note that high-altitude winds flowing from the dayside to the night-side hemisphere of hot Jupiters are predicted by global circulation models, and expected to produce a radial-velocity anomaly of 1–2 km s⁻¹ in the peak position of the CCF (Miller-Ricci Kempton & Rauscher 2012; Showman et al. 2013; Kempton et al. 2014). This is in line with the measurement in this paper, although the signal detected with GIANO is both significantly broader and at lower significance than the combined CO + H₂O detection of the transmission spectrum of HD 189 733 b made in the K band with CRIRES (7.6σ, see Brogi et al. 2016). This contributes to the bigger uncertainties in estimating the net blue-shift from the peak position of the co-added planet signal.

In Fig. 7 we compare the width of our detection on the first night of observations to a simulated profile. The latter is the convolution of the autocorrelation function of the best-fitting model, a Gaussian with FWHM = 6.0 km s⁻¹ approximating the instrument profile of GIANO, and a rotational profile obtained with the rigid-rotation model of Brogi et al. (2016) for an equatorial rotation of 3.4 km s⁻¹ (i.e. their best-fitting value). The agreement between the measured and the theoretical widths of the signal is excellent. This further suggests that no other spurious source of correlated noise contribute to the data. We note that if we increase the level of telluric noise by progressively re-including those discarded orders where the masking of telluric lines did not work, and/or by reducing the threshold of the mask (see details in Sect. 3.1), we observe a progressive broadening of the planet signal especially in the red wing, which is coincident with the position of telluric lines (red-shifted) relative to the planet signal on July 11 (see Fig. 8). The correspondence between the expected and the measured width of the signal is hence an additional indicator for an adequate correction of residual telluric signal.

Fig. 7. Shape of the measured planet signal from the first night of observations (black line) compared to the theoretical analogue (cyan line) computed by convolving the autocorrelation function of the bestfitting model, a Gaussian instrumental profile with FWHM = 6.0 km s−1, and a rotational profile (rigid-body rotation, from Brogi et al. 2016) with equatorial velocity of 3.4 km s−1 and day-to-night wind speed of −1.6 km s−1. The agreement between the two profiles is remarkable.

Fig. 8. Radial velocity of HD 189 733 b as a function of planet orbital phase, relative to the observer. This is the sum of the systemic velocity of the system, and the velocity of the barycentre of the solar system compared to Earth. Symbols and colours are the same as in Fig. 1. It shows that on the second observing night (light blue) the planet has zero radial velocity close to mid-transit, i.e. when the strength of the transmission spectrum of HD 189 733 b is maximum. We comment on the differences between the two observing nights in Sect. 4.2.

The marginal detection of a global blue shift on these data suggests that the lower spectral resolution of GIANO compared to CRIRES is unlikely to provide meaningful constraints on the planet rotational rate and/or equatorial winds, which were already hard to constrain with CRIRES at a significance of 7.6σ. We therefore refrain from applying a more sophisticated model for the planet rotation in this study.

4.1. Residual telluric absorption and alternative analysis techniques

In Sect. 3.2, we showed that telluric lines need to be masked even after they are de-trended by airmass to achieve a cross-correlation signal free from telluric contamination. We explained that this is likely due to rapid modulations in the flux of telluric lines due to changes in the water vapour content and possibly instrumental profile, which we cannot model with low-order polynomials. However, these variations are largely common-mode between all the spectral channels affected by telluric lines. This opens up to the possibility of using de-trending algorithms such as principal component analysis (PCA) to model and remove these correlated (in time and in wavelength) signals.

A PCA algorithm was successfully applied by de Kok et al. (2013) to CRIRES data, and more recently Piskorz et al. (2016, 2017) have also utilised it – although along the spectral direction rather than the temporal direction – to clean their data from telluric contamination. We have developed our own version of PCA for de-trending GIANO data. As for the standard pipeline, we process each order separately and then co-add their CCFs at a later stage. So far our results have been sub-par compared to the masking technique, mainly due to the fact that the algorithm struggles in isolating telluric residuals from other time-dependent noise components:

• for an unweighted PCA (i.e. where each spectral channel is equally weighted), the most significant eigenvectors do not contain telluric residuals, but rather trends in spectral channels with very low flux (for instance at the centre of a saturated telluric line) and broad-band variations of the throughput (for instance due to the instrument modal noise);

• even when deep absorption lines are masked and the signal-to-noise of each channel is used as weighting, the eigenvectors still model broad-band sources of correlated noise. The main reason for this is that broad-band variations affect a larger portion of the data than narrow telluric lines;

• if more eigenvectors are linearly combined and removed from the data, telluric residuals are eventually suppressed, but so is the planet signal, resulting in a marginal detection around S/N = 3.

4.2. The non-detection from the second transit

Although the weather and the seeing were good on both nights, several indicators point to an inferior quality of the data taken on July 30. Firstly, the bottom-right panel of Fig. 2 shows that the shifts in the dispersion solution as recorded on the GIANO detector are one order of magnitude bigger than on the first night, pointing to a lower stability of the spectrograph. Counterintuitively, we suspect that this is due to a particularly good seeing. This enhances the non-uniform illumination affecting the fibres of GIANO, and indeed we also measure an enhanced amount of modal noise in the K-band orders (0–7 in our numbering). Secondly, on July 30 the observations were interrupted twice for technical problems with the instrument and its interface. It resulted in a loss of 20% of the planetary transit (see gaps in Figs. 1 and 8) and it potentially impacted the overall stability of the instrument. Finally, Fig. 8 shows the radial velocity of HD 189 733 b in the observer’s reference frame. This also equals the relative velocity shift between telluric and planet spectra. Ideally one would want to have this shift be as large as possible to minimise the contamination between the two spectra. During the first night of observations (dark blue) the relative radial velocity between HD 189 733 b and the Earth was largely non-zero, except for part of the egress. However, during the second observing night (light blue) the planet radial velocity is zero exactly at mid-transit, which is when the strength of the planet’s transmission spectrum is maximum. As a consequence, the masking applied to telluric lines prior to cross-correlation will also potentially affect planet lines, as they are superimposed at mid-transit. The spectra on the second night were also taken at higher airmass on average (Fig. 1), requiring a more aggressive masking. If this choice succeeded in minimising correlated noise at low K _P (Fig. 5, bottom panels), it also eroded a larger fraction of the planet spectrum, decreasing our sensitivity.

We substantiate the claim of an inferior data quality on July 30 by injecting the best-fitting planet model in the spectra from both nights, at three times the nominal level, and at the planet radial velocity obtained from Eq. (2) with the literature values listed in Sect. 3.2. After running the injected spectra through the full data-analysis pipeline, we recover a signal-to-noise ratio of 9.4 and 4.3 on the first and second night, respectively. Given that the real planet signal is detected at S/N = 4.9 on the first night, scaling arguments point to an expected S/N of 2.2 on the second night, well below detectability. Based on the evidence presented in this Section, we conclude that the absence of signal on July 30 is in line with the expectations.

4.3. Comparison with previous detections

It is useful to assess whether the detection presented here is compatible with the previous detection of water at 4.8σ with CRIRES (Brogi et al. 2016). For GIANO, we measure 2200e⁻ per 60 s integration per resolution element (1 pixel = 2.7 km s⁻¹) in the K band. With CRIRES, we measure 6000e⁻ per 10 s integration, per resolution element (1 pixel = 1.5 km s⁻¹). Scaling the latter signal based on telescope size and pixel width while assuming equal throughput, we would expect from GIANO approximately 12 500e⁻ per 60 s integration per pixel. Observed GIANO spectra are thus deficient by a factor of (12.5/2.2) = 5.7, which we attribute to a much lower throughput. Based on flux alone, GIANO should only produce 40% of the S/N of CRIRES for the same exposure time, insufficient to detect the transmission spectrum of HD 189 733 b. However, we should also consider the increased spectral range of GIANO. We estimated this additional factor by adding random noise to the best fitting model for the planet transmission spectrum, cutting it accordingly to the wavelength range of GIANO and CRIRES, and cross-correlating it with the noiseless model. For each of the GIANO orders, the noise is rescaled to properly account for the wavelength-dependent transparency of the Earth’s atmosphere. We repeat the simulations 100 times, and we compare the distributions of the retrieved S/N for CRIRES and GIANO. On average, GIANO delivers 4.5× the S/N of CRIRES based on the spectral range alone. However, GIANO has half of the spectral resolution, producing a loss of 30% in S/N for unresolved planetary lines. Putting all the scaling factors together, H₂O should be detectable with GIANO in one transit at a significance of 4.8σ × 4.5 × 0.7 × 0.4 = 6σ. This is consistent with our reported detection at 5.5σ. Not only does this show that a detection at high significance is possible with just 110 min of TNG time, but also that the overall strength of the planet’s transmission spectrum has not changed significantly since the previous CRIRES observations from July 2012. Future monitoring of the planet via the same technique will allow us to put constraints on the level of variability of the NIR transmission spectrum due for instance to aerosols at the planet terminator.

5. Discussion

We have tested whether or not a 4-m telescope with a performing high-resolution spectrograph can successfully study the atmospheres of exoplanets at high spectral resolution. We have shown that GIANO in its initial fiber-fed configuration, despite the much smaller throughput, is capable of rivalling CRIRES at the VLT in detecting molecular absorbers with broad-band opacity, such as water vapour. We have demonstrated that the analysis devised for CRIRES still works at half of the spectral resolution, but residual telluric residuals need to be treated more aggressively (e.g. masked) in order to avoid contamination of the planet signal. Smart scheduling of the observations in order to maximise the modulus of the combined barycentric plus systemic velocity of the target is strongly advised to further reduce time-correlated noise coming from uncorrected telluric absorption.

GIANO has been recently updated to a slit-fed instrument and coupled to HARPS-N to deliver quasi-contiguous wavelength coverage from 0.383 to 2.45 µm (Claudi et al. 2017). In this configuration, the instrument will allow us to access in a single observation much smaller atmospheric pressures through sodium absorption (Wyttenbach et al. 2015; Louden & Wheatley 2015). Furthermore, it will be possible to search for the signature of the weaker absorption bands of water vapour (Allart et al. 2017; Esteves et al. 2017), which are very sensitive to the presence of clouds or hazes. On the instrumental side, this upgrade enhanced the throughput of GIANO by at least a factor of four. Furthermore, it eliminated the modal noise affecting the orders in the K band (approximately orders 0–7 in our analysis, which we largely discarded). When the simulations of Sect. 4.3 are updated to account for these better performances, and neglecting decreased performances in correcting for telluric lines, we conclude that GIANO could potentially deliver 2.5–3× the S/N of CRIRES at the VLT in the same observing time, when looking for broad-band absorbers such as H₂O or CH₄. This advantage could be used either to open up high-resolution spectroscopy to fainter systems or to start investigating the atmospheres of smaller and cooler exoplanets.

The planned upgrade for CRIRES (CRIRES+, see Follert et al. 2014) will transform the instrument into a cross-dispersed spectrograph, bringing a 10× increase in spectral range. Together with state-of-the-art detectors, a refurbished AO system, and improved calibration sources, CRIRES+ is likely to overcome any gaps in sensitivity with GIARPS. Whereas in the past highdispersion observations have been limited to either the southern or the northern hemisphere, with both spectrographs online we will be able to access the entire sky.

Adding to the fleet of near-infrared, high-resolution spectrographs, IGRINS (Park et al. 2014) and CARMENES (Quirrenbach et al. 2014) have recently come online, and SPIROU (Artigau et al. 2014) will be available soon. These instruments are revolutionary both in terms of spectral range and throughput compared to the old CRIRES. If the NIR sky at their observatories is of comparable quality to La Palma or Cerro Paranal, and if their initial specifications will be confirmed (especially in terms of throughput), these instruments will be well suited to perform studies of exoplanet atmospheres such as the one presented here.

Importantly for future comparative characterisation of exoplanets, spectra taken at high- and low-resolution contain highly complementary information (Pino et al. 2018). Whereas at high resolution temperature and opacity in the planet’s atmosphere are encoded in the relative ratio between hundreds of thousands of spectral lines, at low resolution the same information translates into broad-band spectral variations. Brogi et al. (2017) recently demonstrated that constraints on abundances, temperature–pressure (T – p) profile, and consequently metallicity of exoplanets can be greatly enhanced when combining the confidence intervals from space (HST/WFC3 and Spitzer/IRAC data) and ground observations (CRIRES). To enable this synergy, especially in the era of TESS and JWST, it is necessary to extend high-resolution spectroscopic observations of exoplanets to as many ground-based telescopes as possible. Based on the results of this paper, and considering the increased ease of scheduling, a 3–4-m telescope equipped with a state-of-the art high-resolution spectrograph could not only assist, but even outperform the biggest ground-based telescopes for detecting key molecular constituents such as H₂O and CH₄.

The above observations make use of both spectral and temporal resolution to extract exoplanet signals from the overwhelming stellar and telluric contaminants. With a total integration time of ten hours, core-to-wings line ratios as small as 10⁻⁵ (1σ) have been detected through their combined cross-correlation signal. Contrast ratios exceeding a billion could potentially be reached if the spatial separation is also used, i.e. by coupling a high-resolution spectrograph to a coronographic or extreme-AO system (Snellen et al. 2015). Early demonstrations of the principle were recently presented (Snellen et al. 2014; Schwarz et al. 2016), and the prospect of targeting biomarkers in terrestrial exoplanets orbiting in the habitable zones of M-dwarfs seems feasible when the next generation of extremely large telescopes come online (Snellen et al. 2013; Rodler & López-Morales 2014).

Acknowledgments

We thank the referee M. Kuerster for his insightful comments, which contributed to improving the quality of the manuscript. We thank R. Claudi and S. Benatti for their insightful discussion about GIARPS. We thank F. Borsa, A. Maggio, N. Nikolov, and I. Pagano for their initial contribution to the design of these high resolution observations. M.B. acknowledges partial support by NASA, through Hubble Fellowship grant HST-HF2-51336 awarded by the Space Telescope Science Institute. P.G. gratefully acknowledges support from INAF through the Progetti Premiali funding scheme of the Italian Ministry of Education, University, and Research. G.G. acknowledges the financial support of the 2017 PhD fellowship programme of INAF. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007–2013) under Grant Agreement No. 313014 (ETAEARTH).

References

All Figures

	Fig. 1. Orbital phase of exoplanet HD 189 733 b and corresponding airmass during the two nights of observations analysed in this work. Plus and star symbols denote out-of-transit and in-transit data, respectively.
In the text

Fig. 2. Stability of the GIANO spectrograph. The dispersion solution of the spectra taken at A and B position differs as a function of the spectral order considered (left panel, averaged over the observing night). In addition, we measure temporal variations in the global dispersion solution (right panel). Their amplitude differs between the first and the second night of observations (top-right and bottom-right panels, respectively). These measurements suggest that although A and B spectra move coherently (i.e. their relative shift stays the same), their absolute position is variable on a timescale of minutes. Furthermore, the overall stability of GIANO differs between observing nights.

In the text

	Fig. 3. Residuals of the wavelength calibration (in km s−1) vs. wavelength for the full spectral coverage of GIANO. The dashed lines denote the ±1σ in the distribution of the residuals. For most of the orders, a sub-km s−1 per line wavelength solution is achievable.
In the text

Fig. 4. Cross-correlation between the sequence of spectra of HD 189733 observed with GIANO on July 11, 2015, and a model spectrum for the planet HD 189 733 b containing H2O at VMR = 10−4, shown as function of radial velocity and planet orbital phase. Ingress and egress of HD 189 733 b are marked with horizontal dashed lines. In the left panel, we co-added all the GIANO orders after removing telluric lines similarly to past CRIRES observations (Brogi et al. 2014, 2016). In the right panel, we applied an additional mask at the position of telluric H2O lines. The latter step removes the majority of time-correlated residual telluric noise, leaving a faint trail of positive correlation along the expected planet radial velocity (marked by slanted dashed lines). Details can be found in Sect. 3.1.

In the text

Fig. 5. Detection of H2O in the transmission spectrum of HD 189 733 b observed on July 11, 2015 (top panels) and the non detection on the night of July 30, 2015 (bottom panels). The total cross-correlation signal is shown as function of rest-frame velocity νrest and planet radial velocity semi-amplitude KP. The measured S/N and corresponding significances are shown in the left and right panels, respectively. These are computed as explained in Sect. 3.3. The cross-correlation shown here is obtained with a model spectrum containing H2O at VMR = 10−4 and CH4 at VMR = 10−12.

In the text

Fig. 6. Distribution of cross-cross correlation values away from the radial velocity of HD 189 733 b (out-of-trail values, in blue), compared to a Gaussian distribution with the same sample mean and variance (dashed line, light blue). The top panel is for the first night of observations, the bottom panel for the second. In both datasets, the sample distribution shows no deviation down to approximately four times the standard deviation. For comparison, the corresponding distribution of the cross-correlation values around the planet radial velocity (in-trail) is shown in red. In the data from the first night, the two distributions appear shifted as one would expect if the planet transmission spectrum is detected. Conversely, no evident shift is detected in the data from the second night. The statistical significance of these shifts is quantified through a generalised t-test as explained in Sect. 3.3.

In the text

Fig. 7. Shape of the measured planet signal from the first night of observations (black line) compared to the theoretical analogue (cyan line) computed by convolving the autocorrelation function of the bestfitting model, a Gaussian instrumental profile with FWHM = 6.0 km s−1, and a rotational profile (rigid-body rotation, from Brogi et al. 2016) with equatorial velocity of 3.4 km s−1 and day-to-night wind speed of −1.6 km s−1. The agreement between the two profiles is remarkable.

In the text

Fig. 8. Radial velocity of HD 189 733 b as a function of planet orbital phase, relative to the observer. This is the sum of the systemic velocity of the system, and the velocity of the barycentre of the solar system compared to Earth. Symbols and colours are the same as in Fig. 1. It shows that on the second observing night (light blue) the planet has zero radial velocity close to mid-transit, i.e. when the strength of the transmission spectrum of HD 189 733 b is maximum. We comment on the differences between the two observing nights in Sect. 4.2.

In the text