© ESO 2019

1 Introduction

Ground-based and space surveys have uncovered an astonishing diversity in the physical properties of exoplanets (e.g., Fulton et al. 2017 and references therein). Nevertheless, there are still many unanswered questions relating to their true nature, including their internal composition and its possible dependence on the formation history. Exoplanet atmospheres provide a means to directly address these still open questions.

Atmospheric characterization measurements have been particularly successful for the class of transiting exoplanets. The spectra of transiting planets are obtained by studying how the total flux (star+planet) varies as a function of wavelength, when the planet passes in front of (transmission spectrum) or behind (emission spectrum) its host star along our line-of-sight. The probability that a planet transits is given by P_transit ~ R_⋆∕a, where a is the orbital separation between the star and the planet, and R_⋆ is the radius of the star. For close-in hot Jupiter exoplanets, this probability is on the order of ~10% or less. Unfortunately, exoplanets of this type are generally far away from us; the planets around nearby bright stars typically do not transit. For this reason, we need to devise a method capable of studying the atmospheres of nontransiting planets around bright stars in the solar neighborhood.

In the past, many studies have been carried out trying to characterize the atmosphere of nontransiting planets with the high-resolution spectroscopy (HRS) technique. Collier Cameron et al. (1999) first reported a probable detection of starlight from τ Boötis b. This was followed by other works (Charbonneau et al. 1999; Wiedemann et al. 2001; Leigh et al. 2003; Collier-Cameron et al. 2004; Lucas et al. 2009; Rodler et al. 2010), some very useful since they provided meaningful upper limits that confirmed low albedos for hot Jupiters. But unfortunately, searches for atmospheric compounds in these earlier works did not result in robust detections. A turning point came with Brogi et al. (2012) and Rodler et al. (2012). By using a cross-correlation technique similar to that developed by Snellen et al. (2010) in the case of the transiting planet HD 209458b, they succeeded in detecting a molecule (carbon monoxide) in the atmosphere of a nontransiting planet: the hot Jupiter τ Boötis b.

At highspectral resolution, molecular features are resolved into a dense forest of individual lines in a patternthat is unique for a given molecule (Birkby 2018), meaning that molecular species (H₂ O, CO, CO₂, CH₄, etc.) can be robustly identified by line matching with planetary models. In addition to providing information on atmospheric composition, the HRS technique can also provide an estimate of the true mass and the orbital inclination of nontransiting planets by breaking the M_P sini degeneracy intrinsic to the radial velocity (RV) measurements. In fact, HRS allows the orbital velocity of a planet to be detected in the form of a Doppler-shifted spectrum. In this way, planetary mass can be determined using Newton’s law of gravity, treating the star+planet system as a double-line spectroscopic binary.

Atmospheric characterization studies using HRS for non-transiting planets have been carried out using spectrographs mounted on the biggest telescopes on the ground, such as the Cryogenic InfraRed Echelle Spectrograph (CRIRES) at the ESO Very Large Telescope (VLT; Rodler et al. 2012; Brogi et al. 2013, 2014; Snellen et al. 2014; Schwarz et al. 2016; Birkby et al. 2017) or the Near InfraRed Spectrograph (NIRSpec) at Keck (Lockwood et al. 2014; Piskorz et al. 2016, 2017). In this paper, we show that smaller telescopes, if equipped with performing high-resolution spectrographs, can also be employed for molecular detections in the dayside spectra of nontransiting planets. Our group has worked with the Telescopio Nazionale Galileo (TNG) and we have recently demonstrated that despite its 3.6 m mirror, it can be successfully employed to characterize the atmosphere of transiting planets around bright host stars (Brogi et al. 2018). The TNG is relatively powerful because it is endowed with the high-resolution (R ~ 50 000) near-infrared (NIR) spectrograph GIANO (Oliva et al. 2006; Origlia et al. 2014). GIANO can cover a wide spectral range (0.95–2.45) μm in each individual exposure and provides a 25-fold increase in spectral range compared to CRIRES pre-upgrade at VLT. Furthermore, this large spectral range leads to signal-to-noise ratios (S/N) that allow GIANO to secure solid detections in the NIR.

In this work, we report spectroscopic observations of the dayside of the nontransiting hot giant planet HD 102195b with GIANO, to probe its atmospheric composition and to estimate the planet’s true mass and the inclination angle of the orbital plane.

The paper is organized as follows. In Sect. 2, we describe the system HD 102195. The observations and the standard calibration are discussed in Sect. 3, while in Sect. 4 we report how we extract the planetary signal. The results of the analysis are presented in Sect. 5. The summary and discussion are provided in Sect. 6.

2 The system HD 102195

HD 102195 is a mildly active, relatively bright (V = 8.06 ± 0.01 mag, J = 6.629 ± 0.024 mag, H = 6.268 ± 0.026 mag, K = 6.151 ± 0.018 mag) G8V dwarf (Ge et al. 2006; Melo et al. 2007) that may be part of the Local Association (Ge et al. 2006). We redetermined the stellar parameters using the Yonsei-Yale evolutionary tracks (Demarque et al. 2004) with constraints on the stellar luminosity L = 0.485 ± 0.014 L_⊙ from the Gaia DR2 parallax (Gaia Collaboration 2016, 2018), and the host star effective temperature (T_eff = 5291 ± 34 K) and metallicity ([Fe/H] = 0.05 ± 0.05 dex) as derived by Melo et al. (2007) from the analysis of HARPS spectra. We found M_⋆ = 0.890 ± 0.020 M_⊙ and R_⋆ = 0.837 ± 0.017 R_⊙, which are consistent with the values given in Ge et al. (2006) and Melo et al. (2007).

The discovery of a giant planet around HD 102195 was announced by Ge et al. (2006). It was discovered with the RV measurements gathered with the Exoplanet Tracker instrument, a dispersed fixed-delay interferometer mounted at the Kitt Peak National Observatory (KPNO) telescope. Additional RVs were also gathered with the HARPS spectrograph at the 3.6 m ESO telescope (La Silla) by Melo et al. (2007). Since HD 102195b is a nontransiting planet, only Doppler studies have been possible, from which we have an estimate of its minimum mass M_P sini = 0.488 ± 0.015 M_J (Ge et al. 2006), while its orbital inclination i is undetermined. In this work, we use the HRS technique to determine true planetary mass and system inclination.

To refine the orbital ephemeris, we modeled both the HET/High-Resolution-Spectrograph and ESO/HARPS RVs with a Keplerian orbit by employing a differential evolution Markov chain Monte Carlo Bayesian technique (Eastman et al. 2013) with the same implementation as in Bonomo et al. (2015). We found the improved orbital period P and inferior conjunction time to be P = 4.11390 ± 0.00072 d and T_c = 3897.068 ± 0.061 BJD_TDB, respectively.We did not include the RVs obtained at the KNPO telescopes because they do not yield any significant improvement in the orbital solution, their formal uncertainties being approximately one order of magnitude larger than those of the HRS and HARPS measurements (see Table 1 in Ge et al. 2006). The orbital eccentricity e = 0.056 ± 0.036 is consistent with zero (e.g., Lucy & Sweeney 1971; Bonomo et al. 2017).

3 Observation and standard calibration

We observed HD 102195 with the NIR high-resolution echelle spectrograph GIANO mounted at the Nasmyth-A focus of the TNG. GIANO operated in this configuration at the TNG in the March 2015–August 2016 period. GIANO was fiber-fed through two optical fibers (each with an angular diameter of 1′′) at a fixed distance of 3′′. GIANO covered four spectral bands in the NIR (Y, J, H, K) divided into 47 orders. The spectra were imaged on a 2048 × 2048 detector. Observations of science targets were performed with the nodding on fiber acquisition mode: target and sky were taken in corresponding pairs and alternatively acquired on nodding A and B, respectively, for an optimal subtraction of the detector noise and background. Each order was thus divided into two different tracks: orderA and orderB. Furthermore, the presence of an additional slit divided each position (orderA, orderB) into two tracks: A1, A2, B1, B2.

GIANO observations of HD 102195 were carried out during the nights of 9, 17, and 26 March 2016, collecting a total of 48, 50, and 46 spectra, respectively, each with an exposure time of 200 s. During these nights of observation, three ranges of planetary orbital phase were targeted: φ = 0.46−0.50, φ = 0.40−0.44 and φ = 0.59−0.63 (see Fig. 1). In these phase intervals, the planet is near the superior conjunction where it reaches the maximum rate of change of velocity. The radial component of the orbital motion of the planet changes by tens of kilometres per second, which greatly helps in disentangling the Doppler-shifted planet spectrum from the stationary signals in wavelength due to telluric contamination and stellar spectrum.

Fig. 1 Radial velocity curve for HD 102195b. The planet was observed during the nights of 9 (red), 17 (blue), and 26 (green) March 2016, covering three ranges of planetary orbital phases 0.46–0.50, 0.40–0.44, and 0.59–0.63. During these phase intervals, the planet orbital motion changes by tens of kilometres per second, while telluric and stellar contaminations are stationary signals in wavelength.

3.1 Extraction of 1D spectra

For these observations, we have written our own analysis pipeline in IDL, including the calibration of the data and the extraction of the 1D spectra. A master flat-field is produced by subtracting dark images, taken with the same exposure time, and averaging the difference images in time. The master flat-field is used to identify the position of the orders. The position of the apertures is then refined by fitting each spectral track with four Gaussians per order simultaneously (two sliced images and two nodding positions). The location of the orders is given by the average of the four Gaussian centroids. This operation is carried out for 32 equidistant positions along the dispersion direction. The sampled values are then fitted with a parabolic polynomial on the 2048 pixels of the dispersion direction. Since the spectral tracks are curved they are rectified via spline interpolation.

At this stage, it is possible to work with A and B images (i.e., images taken in the nodding positions A and B, respectively). Each image is then corrected for bad pixels by assigning zero flux to the damaged pixels. In order to remove the thermal background and sky emission lines, it is appropriate to work with difference (A −B) images. The1D spectra are then extracted via optimal extraction (Horne 1986). The B spectra need to be multiplied by − 1 to correct for the initial subtraction of the B images. The optimal extraction does not preserve the fluxes. To overcome this problem every spectrum of each aperture must be scaled so that its average is equal to that of the spectrum obtained with rectangular extraction.

Even if a correction before extracting the spectra has been performed, there could be residual bad pixels to be corrected. Order by order and image by image the remaining bad pixels are identified by applying a sigma clipping to each spectrum (after subtracting out the median-filtered version of the same spectrum). For the two nodding positions, numerical masks are built to account for the bad pixels present in each image and in each order.

3.2 Aligning the spectra

Each spectrum is then divided by its blaze function. Order by order, the blaze function is created fitting a four-order polynomial to the spectrum obtained from the master flat with the same extraction procedure described in Sect. 3.1.

The scientific analysis exploits the fact that the atmosphere transmission spectrum of the Earth and, to a first approximation, the stellar spectrum are stationary signals in wavelength (i.e., they fall always in the same pixels) while during each of the three nights of observations the planetary RV is significantly Doppler-shifted and changes by approximately tens of kilometres per second (i.e., it moves across many pixels during the observing sequence). For this reason, we have to correct for the fact that the wavelength solution of GIANO is not stable during the hours of observation and this causes small shifts in wavelength. We remove this effect by aligning all the observed spectra between them.

For each night, spectral alignment is then performed by computing the cross-correlation between a reference spectrum (taken as the average over time) and each spectrum order by order and for a lag vector of 401 elements between −2 pixels and +2 pixels.

The fractional lag corresponding to the maximum cross-correlation is adopted as the shift of each spectrum compared to the reference one. Figure 2 shows that there is a linear trend as a function of order between the dispersion solution for nodding position A (in black) and that for nodding position B (in red). This effect had already been noted in the GIANO data collected for the HD 189733b atmospheric characterization study (Brogi et al. 2018). Order by order, the trend is removed by re-aligning each spectrum to the reference one by spline interpolation based on the calculated shifts.

3.3 Wavelength solution

The absorption spectrum of the atmosphere of the Earth (telluric spectrum) in the NIR provides an excellent wavelength calibration source. For each order, we visually match the positions of a set of telluric lines in the time-averaged observed spectrum and in a high-resolution atmospheric transmission spectrum of the Earth generated via the ESO Sky Model Calculator¹. The telluric lines are chosen to be isolated and not saturated. The line centroids are then estimated by fitting each line with a Gaussian profile. The centroid positions are computed in pixels for the averaged spectrum, and in wavelength for the corresponding lines in the template spectrum. A first (pixel, wavelength) relation is obtained by fitting the centroids positions with a fourth-order polynomial. We then use this first solution to refine the centroid positions; we compute a super-sampling of each line via spline interpolation at 1/20 of a pixel and we record the fractional pixel and the corresponding wavelength at which the flux reaches a minimum. We re-fit this new (pixel, wavelength) relation with a fourth-order polynomial and we assume the fit to be the wavelength solution.

Unfortunately this type of calibration based on the telluric lines is not possible for all the orders. For those showing only a few or no telluric lines, the wavelength calibration steps are performed by working with the stellar lines instead of the telluric lines, and by using a Phoenix stellar spectrum (Allard et al. 2011), calculated for a star with properties closely matching those of HD 102195, as areference spectrum.

Since there are spectral orders where the atmosphere of the Earth is particularly opaque (i.e., they do not have enough flux or unsaturated spectral lines), and others in which there is a mismatch between modeled and observed spectral lines, we are not able to calibrate all the spectral regions, and we have to exclude orders 8, 9, 10, 23, 24, 30, 41, 43, 44, and 45 from our analysis. For the remaining orders, the precision of our calibration is estimated by computing the standard deviation of the residuals. We achieve a residual scatter per line well below 1 km s⁻¹, that is, approximately one third of a pixel (see Fig. 3).

Fig. 2 For each order, we plot the average shift (diamonds) over time (dots). The drift relative to the reference spectrum is expressed in pixels (1 pixel equals to 2.7 km s−1). In the left panels, for the three nights, we have taken as reference spectrum the average over the spectra taken in position A (in black) while in the right panels the average over the spectra taken in position B (in red). The average shifts of the spectra over time acquired with one nodding position are shifted with respect to those taken in the other position. These shifts present themselves in a linear fashion, as a function of spectral order. The amplitude of the shifts does not depend on the night considered.

4 Extraction of the planetary signal

In ground-based observations, the planetary signal is hidden under telluric contamination. Standard techniques to remove telluric lines typically use the spectrum of a telluric standard star observed before or after each science spectrum. However, because we want to observe how the planet changes its velocity during time and gain as greater a S/N as possible, it is not advisable to interrupt the series of observations to take telluric standards. Moreover, telluric standards are observed in regions of the sky that are differentfrom the studied target and with different airmasses. Since in the NIR the atmosphere varies on timescales of minutes, thecorrection with telluric standards would not reach the level of precision required by our analysis. There is an alternative approach, a key point of the HRS method, that allows us to disentangle the planet and the telluric spectrum and therefore to remove the telluric contamination (Birkby 2018). The planetary signal during the three nights of observation is significantly Doppler-shifted and changes by approximately tens of kilometres per second, while the atmosphere transmission spectrum of the Earth is a stationary signal in wavelength. The telluric spectrum can therefore be removed by modeling the flux of each spectral channel in time. For this analysis, we follow the approach of Brogi et al. (2018).

Fig. 3 Precision of the wavelength calibration for the orders of GIANO considered in the analysis of HD 102195b spectra.

Fig. 4 Example of our data reduction chain to remove telluric contamination: a: extracted spectra. b: data in logarithm space after throughput correction. c: data in logarithm space after airmass correction. d: GIANO residuals after telluric correction.

4.1 Removing the telluric lines

We remove the telluric contamination as follows (see Fig. 4). For each night we consider every order separately as a matrix with wavelength (pixels, 2048 channels) on the x-axis and time (frame number, n_f) on the y-axis. We prefer to treat A and B spectra separately so as to take into account possible variations between the two nodding positions. Therefore, we work with n_ord (orders positions remaining after wavelength calibration) × 2 (nodding position) matrices, each one of 2048 × n_f in size. We first apply the masks built in Sect. 3.1 and we update them including all the pixels that present a flux less than 0.1counts. We perform this operation because we work with flux in logarithmic space (since the depth of telluric lines depends on the exponential of airmass) and in this way we avoid infinite values. We also mask residual damaged pixels previously identified in Sect. 3.1. The masks are re-updated whenever additional bad pixels are identified. Once the natural logarithm of each matrix has been taken, to remove the instrumental throughput (due to pointing, seeing, and sky transparency) we subtract the median of the 300 brightest pixels from each spectrum.

Since the GIANO detector is divided into four sensors, each one sized 1024 × 1024, every order is therefore split into two sensors. To correct for the differences between the right and the left part of each order, pixels 0 through 1023 along the x-axis are separated from pixels 1024 through 2047. For the three nights of observation, we therefore work with n_ord (positions of the orders) × 2 (nodding position) × 2 (detector quadrants) matrices of 1024 × n_f in size. From every spectrum of these matrices a linear fit between the spectrum itself and a time-average spectrum is subtracted out. In this way, we eliminate possible discontinuities between spectra acquired with different nodding positions or with different detector quadrants. The linear fit is performed using the IDL routine ladfit.pro, which fits the paired data (xi, yi) with a linear model, y = A+ Bx, using a robust, outlier-resistant least-absolute-deviation method.

The most dominant feature in the spectra remains the telluric absorption, which has a fixed pattern in wavelength but has a depth which varies with time (due to changes in airmass during the night of observation). For this reason, we remove the telluric contamination by modeling, for every spectral channel, λ_i, how the airmass changes with time. For every order, we first merge the matrices corresponding to the right and left parts of the GIANO quadrant and then we subtract from each spectrum a second-order fit performed between the observed flux (in logarithmic space) and the airmass. Some residuals are still present at the position of some of the strongest absorption lines due to varying conditions in the atmosphere of the Earth during the night. This is expected because the amount of water vapor changes, generally not following the airmass, and it is not distributed homogeneously in the atmosphere, but rather confined to the troposphere of the Earth. We correct for this effect by measuring the flux changes in a few of the deepest H₂ O lines over time, and we use these to correct the rest of the matrix with a column-by-column linear regression (Brogi et al. 2014).

After telluric removal, the bad-pixel/low-count masks are applied again assigning zero flux to these pixels. To eliminate residual strong outliers, a sigma-clipping with a threshold of eight times the standard deviation of each matrix is applied. In the end, the matrices are returned to linear space.

4.2 Atmospheric models

We compute models for the atmosphere of HD 102195b using the line-by-line plane-parallel radiative transfer code of de Kok et al. (2014). The planet atmosphere is divided into 50 layers equally spaced in logarithm of pressure (log p), and the corresponding temperature is interpolated from a prescribed temperature–pressure profile given as input (see Fig. 5). Opacities from H₂ –H₂ collision-induced absorption (Borysow et al. 2001; Borysow 2002) are incorporated. CO, H₂ O, and CH₄ are added as trace species, with line position, intensity, and air-broadening parameters taken from HITEMP2010 (Rothman et al. 2010) for CO and H₂ O, and from HITRAN2012 for CH₄ (Rothman et al. 2013). Due to the overwhelming amount of transitions, the line list for water vapor is peeled to keep the strongest 10 000 lines every 0.042 μm in wavelength, resulting in a total of 380 000 H₂O lines. The line lists for CO and CH₄ do not require peeling, and contain 40 000 and 46 300 lines, respectively, in the range of these GIANO models (0.91–2.50 μm). Being extrapolated from low-temperature measurements, the HITRAN methane list is known to be incomplete, however it contains the majority of the strong lines shorter than 1.6 μm and longer than 1.7 μm, where most of the GIANO data is sensitive to Yurchenko et al. (2014). In order to keep the modeling computationally feasible, we calculatetwo sets of models:

• 1.

  In the first set, water, carbon monoxide, and methane are added as single trace species. The temperature–pressure profile is bounded by two isothermal zones, one at temperature T_base = 1450 K for p > 3 bar, and the other at T = T_up for pressures p < p₁, where log(p₁) = [−4, −3, −2] bar and T_up = [600, 850, 1650] K. Volume mixing ratios for the species range between 10⁻⁶ and 10⁻³ in steps of one decade. We thus compute 36 models per species.

• 2.

  In the second set, we mix the two species detected singularly (water and methane) and we further simplify the temperature–pressure profile by removing the thermal inversion case (T_up = 1650 K), which is only producing anti-correlation from the previous set of models (see Sect. 5.1). We thus produce a set of 54 models.

Fig. 5 Schematic representation of the parametrization used in our models for the planet temperature–pressure profile. The atmosphere is composed of two isothermal layers (at temperatures Tbase and Tup) separated bya layer with constant lapse rate d log(p)∕dt. While Tbase is fixed to 1450 K and the base pressure is fixed to 3 bars, the other quantities are varied in a small grid, as explained in Sect. 4.2.

4.3 Extraction of the signal by cross correlation

With the telluric lines removed as above, we are left with residuals (shown in the lower panel of Fig. 4). Since the individual planet lines have very low S/N (S∕N_line ≪ 1) the planet signal is deeply hidden in the noise. However, there are thousands of molecular lines in the GIANO spectral range, and therefore we can combine their signals to attempt a detection of the planet signature.

The spectral lines can be co-added by cross-correlating the residual GIANO data with a model spectrum appropriate for the planet atmosphere. In this way, observation by observation, we evaluate how well the model matches the residual spectra.

The cross-correlation is computed on a fixed grid of RV lags between −225 and 225 km s⁻¹, in steps of 2.7 km s⁻¹ (167 values), which matches the pixel scale for most of the GIANO orders. We must be careful to choose the step size because if it is taken as smaller than the velocity resolution of the pixel of the detector, we risk an oversampling of the data (Birkby 2018). The cross-correlation is calculated for every order and for every integration separately. Each spectrum is cross-correlated with the model, Doppler shifted for each RV value, and convolved with the instrument profile of GIANO (a Gaussian with FWHM 5.4 km s⁻¹). For every night the output of the cross correlation is a matrix with dimensions 167 × n_f × n_ord.

Each cross-correlation function (CCF) gives an estimate of how well the model matches the data for each RV step; a positive value indicates correlation, while negative values denote anti-correlation. This can arise if we cross-correlate emission lines with models in absorption or vice versa.

At this stage, the planetary signal is not yet visible; to maximize it, we co-add all orders, obtaining three matrices n_f × 167, one for eachnight.

The toppanel of Fig. 6 illustrates the co-added matrix of CCFs for the night of March 26, showing clear telluric residuals, which are at zero velocity in the telluric frame. We therefore apply the masking of Brogi et al. (2018) prior to cross correlation by assigning a zero value to those pixels corresponding to telluric lines deeper than 5–30%, where this threshold varies order by order. This is a procedure that is to be carried out iteratively by visually analyzing the cross-correlation of each order until telluric residuals disappear. Despite this masking, in some orders we continue to see some cross-correlation noise that does not resemble the signature of residual telluric lines (likely affected by e.g., modal noise). We therefore discarded these orders from the analysis. In this way we remain with the following groups of orders [0, 1, 2, 3,32, 33, 37, 38, 40, 41], [0, 2, 7, 17, 18, 19, 25, 27, 29, 31, 32, 33, 34, 38, 40, 41], and [1, 2, 3, 4, 12, 14, 16, 19, 27, 29, 32, 33, 40, 41], for the nights of 09, 17, and 26 March, respectively. Finally, we proceed to co-add the cross-correlation matrices of the three nights as follows. Since every CCF is taken at a different planetary phase and so at a different RV, we must shift all the CCFs in the planet rest frame. Doppler shifting the planet requires the knowledge of its RV (V_rv), the systemic velocity of the star-planet system with respect to the Earth, for which we adopt the Gaia DR2 value V_sys = 1.85 ± 0.15 km s⁻¹ (Gaia Collaboration 2016, 2018), and the velocity of the observer induced by the motion of the Earth around the Sun (V_obs; i.e., the barycentric velocity). Assuming a circular orbit and neglecting cross-terms in velocity, the total planetary velocity V_P (t) is given by: (1)

In the above equation, V_rv = K_P sin[2π φ(t)], where φ(t) is the planetary phase and K_P is the semi-amplitude of the planet RV. Since HD 102195b is a nontransiting planet, sin i is undetermined, and therefore we do not know K_P. We assume a range of possible K_P (0 ≤ K_P ≤ 170 km s⁻¹), which matches all the possible orbital inclinations of the system, plus a small interval of unphysical values (sin i > 1) as a sanity check (Brogi et al. 2014). The phase φ(t) at time t is obtained as the fractional part of (2)

where P is the orbital period and T_c is the time of inferior conjunction (see Sect. 2). Uncertainties in T_c and P translate into a 1σ uncertainty of Δφ = 0.15 on the orbital phases.

For eachvalue of K_P, we align the CCFs in the planet rest frame via linear interpolation (so that they are centered around V_P (t)), and we sumthem up in phase, thus combining the signal from all of the lines present in the spectra. In this way we obtain the total planetary cross-correlation signal as a function of the rest frame velocity V_rest and the planetary semi-amplitude K_P. If the planetary signal is detected, we measure a peak in the cross-correlation map at V_rest ~ 0. Significant deviations far from V_rest ~ 0 indicate that the signal does not have a planetary origin.

Fig. 6 Co-added matrix of CCFs for the night of March 26. In the upper panel at zero RV, telluric residual are visible. We remove them applying a mask to the data before the cross-correlation. The result of this operation is shown in the lower panel.

Fig. 7 Total planetary cross-correlation signal, obtained by combining the three nights, as a function of rest frame velocity Vrest and planetary semi-amplitude KP. Cross-correlation was performed using models containing only one molecule: CH4 (left panel) or H2O (right panel). Significance was calculated using the Welch t-test, as described in Sect. 5.3

5 Results

5.1 Molecule identification

Figure 7 shows the results of the cross-correlation obtained from models containing only one molecule, that is, water vapor or methane. H₂O and CH₄ are detected with a significance of 4.4 and 4.1σ, respectively. The significance level is estimated through the Welch t-test as explained in Sect. 5.3.

Carbon monoxide is not detected in these data at significance above 3σ. Although this might seem counter intuitive compared to past CRIRES studies in which this molecule produced the strongest detections, one ought to consider the difference in sensitivity between CRIRES and GIANO. The latter instrument, being mounted at a smaller telescope and having less resolving power, can only compete with CRIRES by significantly increasing the spectral range to include multiple NIR bands. However, this benefit only applies to the detection of molecular species broadly absorbing in the NIR, particularly water. Carbon monoxide, on the other hand, only absorbs in part of the Kband, and more weakly in part of the H band. At the S/N of these GIANO observations, there are not enough CO lines to build up a significant detection through cross correlation. We verified that injecting CO-only model spectra in the data and attempting to retrieve their signal indeed produces no firm detections, with peak values of significance below 3.0. For these reasons, we do not consider CO for the remainder of the analysis.

We choose the best-matching day-side model by requiring the highest significance of the cross-correlation signal obtained from the combination of the three nights. We can immediately rule out an inverted temperature profile for HD 102195b. Indeed, as the upper panels of Fig. 8 show, when we cross-correlate with grid models with a temperature inversion (i.e., T2 = 1650 K), they do not show significant cross-correlation and they peak at V_rest ≠ 0. Moreover, the emission lines of these models tend to negatively correlate with the absorption lines in the observed planetary spectrum (see e.g., Fig. 9, and the discussion in Birkby et al. 2017; Schwarz et al. 2015). For these reasons we conclude that inverted models in our grid do not correctly represent our data.

After this first screening, if we exclude from our analysis the models which peak far from V_rest ~ 0, and we consider the water vapor molecule alone (the panel at the bottom right of Fig. 8), we find three families of models that appear to behave better than others, that is, those with volume mixing ratio (VMR) ~ 10⁻³, (VMR) ~ 10⁻⁴, and (VMR) ~ 10⁻⁵. These are selected by simultaneously requiring proximity to V_rest ~ 0, high significance, and K_P similar to or less than the orbital velocity of the planet V_P (see Sect. 5.2). If we treat the models containing methane-only (the bottom-left panel of Fig. 8), the models with 10⁻⁵ < VMR(CH₄) < 10⁻³ seem to provide the best agreement.

From the single-molecule analysis, after testing the full set of models described in Sect. 4.2, we find that the cross-correlation signal comes from H₂O and CH₄. The VMRs appear to be between 10⁻³ and 10⁻⁵ for both molecules, perhaps with a (slight) preference towards 10⁻⁵ for CH₄ and towards 10⁻³ for H₂O. For both molecules, the preferred lapse rate brackets the range 150 K dex⁻¹ < dT/dlog₁₀(p) < 425 K dex⁻¹.

Finally, the total cross-correlation signal, obtained from models containing both methane and water vapor, appears to be maximized when we consider a mixed model with VMR(H₂O) ~ 10⁻³ and VMR(CH₄) ~ 10⁻⁵.

If we consider the nights taken individually and their combination, within the small group of models (VMR(H₂O) ~ 10⁻³, VMR(CH₄) ~ 10⁻⁵, 150 K dex⁻¹ < dT/dlog₁₀(p) < 425 K dex⁻¹) the one with dT/dlog₁₀(p) ~ 425 K dex⁻¹ is marginally preferred. To better illustrate the choice of our VMRs, we plot in Fig. 10 the different significance level of the mixed models, assuming a dT/dlog₁₀(p) ~ 425 K dex⁻¹. We do note that multiple models with different physical model parameters yield detections with similar statistical significance, and therefore the best model is preferred only marginally over the others.

Fig. 8 Significance, as a function of the rest-frame velocity Vrest and the planetary semi-amplitude KP, for the grid of models containing only water vapor (right panels) or methane (left panels). From the upper panels, we see that models with inverted T/P profiles (red bullets) do not describe the observations well; indeed they have a low significance and peak far away from Vrest ~ 0, where we would expect a signal of planetary origin. If however we restrict our analysis to the models which peak near V rest ~ 0, we find three families of models that provide better agreement than others both for water vapor (10−5 < VMR(H2O) < 10−3) and for methane (10−5 < VMR(CH4) < 10−3). The different bullets shown in the bottom panels correspond to the different lapse rates used.

5.2 Planetary mass and system inclination

The total cross-correlation signal obtained from the best model (VMR(H₂O) = 10⁻³, VMR(CH₄) = 10⁻⁵, dT/dlog₁₀(p) ~ 425 K dex⁻¹) by combining the three nights of observation is shown in the right panel of Fig. 11. The cross-correlation peak is detected at the orbital semi-amplitude of K_P ~ 128 ± 6 km s⁻¹. A phase shift of Δφ = 0.012 needs to be applied in order to match the signal at zero V_rest. This shift is well within the 1σ uncertainty range of Δφ = 0.15.

The upper panel of Fig. 12 shows the 1D cross-correlation function for the mixed model that best matches the data, assuming a planet orbital velocity corresponding to the cross-correlation peak. For completeness, the middle and lower panels show the 1D CCF computed for methane-only, and water-only models, respectively.

The measurement of K_P allows us to treat the star-planet system as a double-lined spectroscopic binary. We can therefore derive the true mass of HD 102195b from the velocity and mass ratios M_S∕M_P = K_P∕K_S. Given the stellar RV semi-amplitude K_S = 63.4 ± 2.0 m s⁻¹ (Ge et al. 2006), its mass M_S = 0.89 ± 0.02M_⊙, and our planetary semi-amplitude K_P = 128 ± 6 km s⁻¹, we obtain a first estimate of the planetary mass of M_P = 0.46 ± 0.03 M_J. The major source of uncertainty in the derivation of the planetary mass comes from the error on K_P which represents about 47% of the indetermination in the calculation of M_P.

The final step is to derive the orbital inclination. Indeed, using Kepler’s Third Law, the orbital velocity of the planet is found to be V_P = 128.0 ± 1.0 km s⁻¹, and under the assumption of a circular orbit sini = K_P∕V _P, therefore allowing us to determine sini = 1.00 ± 0.05. This is translated into 1σ, 2σ, and 3σ lower limits on the orbital inclination of 72.5°, 65.0°, and 59.2° respectively.The nondetection of the planet in previous photometric monitoring (Ge et al. 2006) places an upper limit on the system inclination of i_max = 84.79° ± 0.09°. We note that the resulting error on the assumed RV of the planet is 1.0 km s⁻¹, and being on the order of one third of the GIANO pixels, it does not significantly affect the co-addition of the planetary signal.

Fig. 9 Left panels: total cross-correlation values in the neighborhood of V rest ~ 0. Right panels: theoretical models they correspond to. We use the best models containing methane only (in green, first panels), and water vapour only (in blue, second panels), and two models with temperature inversion with methane only, VMR ~ 10−5, (in green, third panels), and water vapor only, VMR ~ 10−3, (in blue, fourth panels). Since the emission lines in the models with temperature inversion correlated negatively with the absorption lines in the observed planetary spectrum, unlike the absorption lines in the noninverted models, we can infer a noninverted temperature profile for HD 102195b.

Fig. 10 Significance as a function of rest frame velocity Vrest and planetary semi-amplitude KP, by assuming a dT/dlog10(p) ~ 425 K dex−1. From this mixed-models analysis, we found that the model with VMR(H2O) ~ 10−3 and VMR(CH24) ~ 10−5 seems to be marginally preferred over the others.

Fig. 11 Significance of the cross-correlation detection. Right panel: difference between the aligned CCF distributions “in-trail” (red-line) and “out-of-trail” (black-line) for the derived KP. The two distributions are normalized and the error bars correspond to the square-root of the bin occurrences. The in-trail distribution is centered to higher values than the out-of-trail distribution with a 5.3 σ confidence level. Left panel: results of the Welch t-test performed for the same range of vrest and KP as before.

Fig. 12 Total cross-correlation signal as a function of Vrest from all spectra taken during the three nights combined assuming a planet orbital velocity corresponding to the cross-correlation peak. We use our best water only model, methane only model, and mixed models, in the lower, middle, and upper panels, respectively.

5.3 The significance of the cross-correlation detection

To secure the spectroscopic detection of the exoplanet we have to determine the significance level of the signal. Following Brogi et al. (2012, 2013, 2014, 2018), we compare the “in-trail” values in each 2D aligned CCFs’ matrix to those “out-of-trail”. The right-hand side of Fig. 11 shows that the in-trail distribution, for the derived K_P (i.e., those within 4 pixels of the planetary RV curve) is shifted toward higher values with respect to the out-of-trail distribution (i.e., those beyond 10 pixels from the planet RV curve). In order to quantify this, a Welch t-test on the two samples is generated (Welch 1947), with the null hypothesis H₀ being that if they are drawn from the same parent distribution, they have the same mean. The corresponding t value is converted into a significance at which H₀ is rejected, based on the results of the test, that is 5.3σ. The left panel of Fig. 11 shows also the test performed for the same range of v_rest and K_P used before. Except for the peak of the cross-correlation map, no other region of the parameter space shows significant cross-correlation signal.

We may think that atmospheric models containing water vapor and methane would give a cross-correlation signal with a significance that corresponds to the quadrature sum of those obtained from models with H₂ O or CH₄ only; in this case we would obtain a significance of ~6σ. As we see above, this is not the case, as we obtain a lower combined signal at 5.3σ. However, this lower significance can be explained by noting that the two cross-correlation signals are not completely independent, as water and methane lines can overlap at the instrumental resolution of our observations, reducing the total amount of strong lines contributing to the cross-correlation function.

6 Summary and discussion

We studied the atmosphere of the nontransiting gas giant HD 102195b at high spectral resolution using the NIR spectrograph GIANO. Weobserved the system on three nights in March 2016, during which the planet was near the superior conjunction. This helped us todistinguish between the telluric contamination and the planetary signal. In this way, we detected absorption from water vapor and from methane in the planetary day-side hemisphere.

After correcting the telluric lines in order to detect the planetary signal, we cross-correlated the GIANO residual data with a wide grid of models for the planetary atmosphere. We find that the best model has minimum molecular abundances of VMR(H₂O) = 10⁻³ and VMR(CH₄) = 10⁻⁵ and a steep lapse rate of dT/dlog₁₀(p) ~ 425 K dex⁻¹. Our results are in general agreement with the 2D classification scheme suggested by Madhusudhan (2012). Indeed, for temperatures close to that of HD 102195b, they predict a difference of ~ 10⁻² between the two VMRs. In particular, a cursory comparison with chemical model predictions published in the latter-cited publication allows us to infer that a VMR difference of ~10⁻² is compatible with a possibly oxygen-rich atmosphere for HD 102195b, with a C∕O value qualitatively comparable to that of the Sun.

Moreover, our analysis indicates that models with a thermal inversion are disfavored. This result is in line with theoretical predictions. In fact, the host star being active, (Melo et al. 2007), according to Knutson et al. (2010), we should not have expected a temperature inversion, because the UV flux of active stars might destroy the molecular species responsible for this mechanism. Also, a noninverted atmosphere is expected given the relatively low temperature of the planet; indeed we believe that the condensation of TiO/VO plays an important role, potentially able to remove these species in the gas phase from the atmosphere.

The results of the cross-correlation analysis presented here also allowed us to derive the semi-amplitude of the planetary RV (K_P = 128 ± 6 km s⁻¹), which could be translated into a value for the true planetary mass (M_P = 0.46 ± 0.03 M_J), and into an estimate of the system inclination (72.5° < i < 84.79°), with the upper limit set by the nondetection of the planetary transit in previous photometric monitoring (Ge et al. 2006).

In this paper we not only demonstrate that GIANO is effective in studying the atmosphere of nontransiting planets, but thanks to its wide spectral coverage we also attain the first detection of methane in a hot Jupiter with the HRS technique. In past studies, methane had been sought in the atmospheres of exoplanets using both low- and high-spectral-resolution data, but without clear success². These difficulties could possibly be reconciled with the detection presented here if we consider that low-spectral-resolution space-based data obtained with HST WFC3 do not extend longward of 1.7 μm, not reaching the K-band and part of the H-band where the strongest methane lines are found. Instead, at high spectral resolution, CRIRES pre-upgrade and NIRSpec do have wide spectral coverage of (1.0–5.0) μm and (0.95–5.4) μm, respectively, but during one observation only a much smaller spectral region than that covered with GIANO is effectively used.

In addition to spectral coverage and resolution, there are two further aspects that could explain our success in detecting methane. The first is the differences in the fine detail of the data analysis that we have used and the second is the chemistry of the planet being studied. However, it is difficult at the moment, given the small number of observed comparison objects, to determine which element plays the most important role. A greater number of observed planets could help to clarify this aspect. However, what we can conclude is that our CH₄ detection cannot be explained by advances in the line-list quality. Indeed, we use a line list originally compiled for low-temperature applications that certainly does not provide a decisive improvement (in terms of e.g., line-list inaccuracies) with respect to those used in previous analyses that ended without detections. That being said, our line list is in need of improvement. Indeed, having been built for low-temperature applications, it is likely to miss spectral lines that are too weak to be relevant at such temperatures, but that are arguably important at the high temperatures of hot Jupiters. Recently, high-temperature methane line lists have become available (e.g., Yurchenko et al. 2014) and we are planning to expand our modeling framework to incorporate them in future work.

The presence of methane in the atmosphere of a hot Jupiter is not particularly unexpected if we look at theoretical prediction papers (e.g.,Madhusudhan 2012; Madhusudhan et al. 2016; Moses et al. 2011; Drummond et al. 2018; Tsai et al. 2018), although there is still a variety of possible predictions on the expected methane VMR. Our results appear consistent with recent studies conducted on atmospheres similar to that of HD 102195b by Drummond et al. (2018); Tsai et al. (2018).

As mentioned, the lack of coverage of some spectral regions could be a problem for some spectrographs in the search for methane. Unfortunately, even for GIANO, it can undermine the research of certain molecules, in particular carbon monoxide. Indeed, we argue that the nondetection of this molecule in the GIANO spectra is primarily due to the fact that CO absorbs only partially in the H and K bands and has a lower density of spectral lines.

The results presented here constitute a new demonstration of the efficacy of GIANO for the purpose of atmospheric characterization measurements, but they come with an important caveat. Overall, the useful spectral range is somewhat reduced with respectto that used in the Brogi et al. (2018) study, with fewer orders successfully calibrated in wavelength and a larger number of discarded orders in the cross-correlation analysis. These differences, as far as the analysis of pre-upgrade GIANO observations are concerned, can be understood once the following considerations are taken into account. First, the investigation carried out in this work is based on an independently developed analysis pipeline that is similar in philosophy to that used in Brogi et al. (2018), but differs in the fine structure details. Second, compared to the CalarAlto high-Resolution search for M dwarfs with Exoearths with Nearinfrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2014) and to the SpectroPolarimetre Infra-Rouge (SPIRou; Artigau et al. 2014), for example, GIANO in its fiber-fed configuration has a low instrumental stability, which prevents us from using offlinecalibration frames (e.g., lamp spectra taken before or after the observations) to obtain an accurate absolute wavelength solution across all orders. GIANO data rely instead on simultaneous calibration using telluric and (where possible) stellar lines. This means that the effectively usable spectral range can be a function of the night conditions (e.g., variable humidity affects the strength of telluric lines), the achieved S/N per spectrum (which impacts the quality of line centroid determination particularly in the spectral regions where modal noise is strongly present), and the spectral type of the exoplanet host star. As mentioned in Brogi et al. (2018), important improvements are expected from the recent GIANO upgrade to a slit-fed instrument (GIANO-B) at the Nasmyth B focus of the TNG, which can be coupled to HARPS-N in the GIARPS simultaneous observing mode (Claudi et al. 2017).

Finally, our analysis could have important consequences for the future characterization of exoplanetary atmospheres. In fact, HD 102195 being between one and three magnitudes fainter in K-band compared to nontransiting systems studied until now with 8 m telescopes, our study opens the doors to atmospheric characterization measurements with 4 m class telescopes for an important sample of exoplanets. We are now entering a new era for the characterization of exoplanet atmospheres with NIR, high-resolution spectrographs. CARMENES has started to give its first results (Nortmann et al. 2018; Alonso-Floriano et al. 2019), GIANO-B has recently come online, and SPIRou, the Near Infra Red Planet Searcher (NIRPS), and CRIRES+ will be available soon. These instruments are extraordinary both in terms of spectral range and throughput, and will be very well suited to performing studies of exoplanetary atmospheres such as the one presented here. The systematic measurements of molecular compounds at high spectral resolution will permit us to create an “atlas” of strongly irradiated exoplanetary atmospheres, and in the future, with the coming of the European Extremely Large Telescope (E-ELT), it will be possible to apply the same HRS technique to search for biomarkers in the atmospheres of Earth twins (Snellen et al. 2013; Rodler & López-Morales 2014).

Acknowledgements

We wish to thank an anonymous referee for many useful suggestions, which helped to improve an earlier version of the manuscript. We thank F. Borsa, I. Carleo, M. Damasso, R. de Kok, A. Maggio, I. Pagano for lending initial impetus to this analysis. G.G. acknowledges the financial support of the 2017 PhD fellowship programme of INAF. P.G. gratefully acknowledges support from the Italian Space Agency (ASI) under contract 2018-24-HH.0. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

All Figures

Fig. 1 Radial velocity curve for HD 102195b. The planet was observed during the nights of 9 (red), 17 (blue), and 26 (green) March 2016, covering three ranges of planetary orbital phases 0.46–0.50, 0.40–0.44, and 0.59–0.63. During these phase intervals, the planet orbital motion changes by tens of kilometres per second, while telluric and stellar contaminations are stationary signals in wavelength.

In the text

Fig. 2 For each order, we plot the average shift (diamonds) over time (dots). The drift relative to the reference spectrum is expressed in pixels (1 pixel equals to 2.7 km s−1). In the left panels, for the three nights, we have taken as reference spectrum the average over the spectra taken in position A (in black) while in the right panels the average over the spectra taken in position B (in red). The average shifts of the spectra over time acquired with one nodding position are shifted with respect to those taken in the other position. These shifts present themselves in a linear fashion, as a function of spectral order. The amplitude of the shifts does not depend on the night considered.

In the text

	Fig. 3 Precision of the wavelength calibration for the orders of GIANO considered in the analysis of HD 102195b spectra.
In the text

	Fig. 4 Example of our data reduction chain to remove telluric contamination: a: extracted spectra. b: data in logarithm space after throughput correction. c: data in logarithm space after airmass correction. d: GIANO residuals after telluric correction.
In the text

Fig. 5 Schematic representation of the parametrization used in our models for the planet temperature–pressure profile. The atmosphere is composed of two isothermal layers (at temperatures Tbase and Tup) separated bya layer with constant lapse rate d log(p)∕dt. While Tbase is fixed to 1450 K and the base pressure is fixed to 3 bars, the other quantities are varied in a small grid, as explained in Sect. 4.2.

In the text

	Fig. 6 Co-added matrix of CCFs for the night of March 26. In the upper panel at zero RV, telluric residual are visible. We remove them applying a mask to the data before the cross-correlation. The result of this operation is shown in the lower panel.
In the text

	Fig. 7 Total planetary cross-correlation signal, obtained by combining the three nights, as a function of rest frame velocity Vrest and planetary semi-amplitude KP. Cross-correlation was performed using models containing only one molecule: CH4 (left panel) or H2O (right panel). Significance was calculated using the Welch t-test, as described in Sect. 5.3
In the text

Fig. 8 Significance, as a function of the rest-frame velocity Vrest and the planetary semi-amplitude KP, for the grid of models containing only water vapor (right panels) or methane (left panels). From the upper panels, we see that models with inverted T/P profiles (red bullets) do not describe the observations well; indeed they have a low significance and peak far away from Vrest ~ 0, where we would expect a signal of planetary origin. If however we restrict our analysis to the models which peak near V rest ~ 0, we find three families of models that provide better agreement than others both for water vapor (10−5 < VMR(H2O) < 10−3) and for methane (10−5 < VMR(CH4) < 10−3). The different bullets shown in the bottom panels correspond to the different lapse rates used.

In the text

Fig. 9 Left panels: total cross-correlation values in the neighborhood of V rest ~ 0. Right panels: theoretical models they correspond to. We use the best models containing methane only (in green, first panels), and water vapour only (in blue, second panels), and two models with temperature inversion with methane only, VMR ~ 10−5, (in green, third panels), and water vapor only, VMR ~ 10−3, (in blue, fourth panels). Since the emission lines in the models with temperature inversion correlated negatively with the absorption lines in the observed planetary spectrum, unlike the absorption lines in the noninverted models, we can infer a noninverted temperature profile for HD 102195b.

In the text

	Fig. 10 Significance as a function of rest frame velocity Vrest and planetary semi-amplitude KP, by assuming a dT/dlog10(p) ~ 425 K dex−1. From this mixed-models analysis, we found that the model with VMR(H2O) ~ 10−3 and VMR(CH24) ~ 10−5 seems to be marginally preferred over the others.
In the text

Fig. 11 Significance of the cross-correlation detection. Right panel: difference between the aligned CCF distributions “in-trail” (red-line) and “out-of-trail” (black-line) for the derived KP. The two distributions are normalized and the error bars correspond to the square-root of the bin occurrences. The in-trail distribution is centered to higher values than the out-of-trail distribution with a 5.3 σ confidence level. Left panel: results of the Welch t-test performed for the same range of vrest and KP as before.

In the text

	Fig. 12 Total cross-correlation signal as a function of Vrest from all spectra taken during the three nights combined assuming a planet orbital velocity corresponding to the cross-correlation peak. We use our best water only model, methane only model, and mixed models, in the lower, middle, and upper panels, respectively.
In the text