A Principal Component Analysis-based method to analyse high-resolution spectroscopic data

TL;DR

This paper introduces an automatic pipeline utilizing PCA and CCF for analyzing high-resolution spectroscopic data of exoplanet atmospheres, successfully detecting key molecules and aligning with existing results.

Contribution

A novel PCA-based analysis pipeline for high-resolution spectroscopic data, improving automation and validation in exoplanet atmospheric studies.

Findings

Successful detection of CO and H₂O in exoplanet atmospheres

Pipeline results agree with previous studies

Enhanced automation in spectral data analysis

Abstract

High-Resolution Spectroscopy (HRS) has been used to study the composition and dynamics of exoplanetary atmospheres. In particular, the spectrometer CRIRES installed on the ESO-VLT has been used to record high-resolution spectra in the Near-IR of gaseous exoplanets. Here we present a new automatic pipeline to analyze CRIRES data-sets. Said pipeline is based on a novel use of Principal Component Analysis (PCA) and Cross-Correlation Function (CCF). The exoplanetary atmosphere is modeled with the $\tau$-REx code using opacities at high temperature from the ExoMol project. In this work, we tested our analysis tools on the detection of CO and H$_2$O in the atmospheres of the hot-Jupiters HD209458b and HD189733b. The results of our pipeline are in agreement with previous results in the literature and other techniques.

Tables

Table 1. Table 1: Relevant parameters of the studied targets

Parameter	HD189733	HD209458
Stellar Parameters
$R_{\star }$ (R⊙)	$0.756\pm 0.018$ 11footnotemark:	(${1.155}_{-0.016}^{+0.014}$) 11footnotemark:
$T_{eff}$ (K)	$5040\pm 50$ 11footnotemark:	$6065\pm 50$ 11footnotemark:
$M_{\star }$ (M⊙)	$0.806\pm 0.048$ 11footnotemark:	$1.119\pm 0.033$ 11footnotemark:
$log\left(g_{\star }\right)$ (csg)	$4.587\pm 0.015$ 11footnotemark:	$4.361\pm 0.008$ 11footnotemark:
$v_{sys}$ (kms-1)	$-2.361\pm 0.003$ 22footnotemark:	$-14.7652\pm 0.0016$ 55footnotemark:
Planet Parameters
$T_{eq}$ (K)	$\left({1201}_{-12}^{+13}\right)$ 11footnotemark:	$1449\pm 12$ 11footnotemark:
$a$ (AU)	$0.03120\left(27\right)$ 33footnotemark:	$\left({0.04707}_{-0.00047}^{+0.00046}\right)$ 11footnotemark:
$R_p$ (RJup)	$\left({1.178}_{-0.023}^{+0.016}\right)$ 33footnotemark:	$\left({1.359}_{-0.019}^{+0.016}\right)$ 11footnotemark:
$M_p$ (MJup)	$\left({1.144}_{-0.056}^{+0.057}\right)$ 11footnotemark:	$0.685\pm 0.015$ 11footnotemark:
$P$ (days)	$2.21857567\left(15\right)$ 44footnotemark:	$3.52474859\left(38\right)$ 66footnotemark:
$T_0$ (BJDUTC)	$2454279.436714\left(15\right)$ 44footnotemark:	$2452826.629283\left(87\right)$ 66footnotemark:
$I$ (deg)	$85.710\pm 0.024$ 44footnotemark:	$86.71\pm 0.05$ 11footnotemark:

Table 2. Table 2: This work and previous results

Parameter	HD189733	HD209458
Previous results	Brogi et al. (2016)	Snellen et al. (2010)
S/NCO	-	-
$K_{p,CO}$ (kms-1)	${205}_{-51}^{+38}$	-
$v_{p,CO}$ (kms-1)	-	$140\pm 10$
$v_{rest,CO}$(kms-1)	$-{1.6}_{-1.8}^{+2.0}$	$\sim$ 2
S/N${}_{H{}_2{}^{}O}$	5.5	-
$K_{p,H{}_2{}^{}O}$ (kms-1)	${183}_{-59}^{+38}$	-
$v_{rest,H{}_2{}^{}O}$(kms-1)	$-{1.58}_{-1.50}^{+1.65}$	-
Results	This Work	This Work
S/NCO	5.24	5.7
$K_{p,CO}$ (kms-1)	$190\pm 16$	${148}_{-15}^{+16}$
$v_{rest,CO}$(kms-1)	$-{3.0}_{-1.3}^{+1.0}$	$-{3.0}_{-1.1}^{+1.3}$
W T-TestCO ($\sigma$)	-	21.62
S/N${}_{H{}_2{}^{}O}$	3.69	3.95
$K_{p,H{}_2{}^{}O}$ (kms-1)	${167}_{-21}^{+32}$	${140}_{-16}^{+25}$
$v_{rest,H{}_2{}^{}O}$(kms-1)	$-{4.0}_{-1.8}^{+2.0}$	$-{4.0}_{-1.6}^{+1.4}$
W T-Test${}_{H{}_2{}^{}O}$ ($\sigma$)	5.21	6.56

Equations

$$CCF(d)=\frac{{\sum }_i\left(\left(x(i)-\overline{x}\right)\cdot \left(y(i-d)-\overline{y}\right)\right)}{\sqrt{{\sum }_i{\left(x(i)-\overline{x}\right)}^2}\cdot \sqrt{{\sum }_i{\left(y(i-d)-\overline{y}\right)}^2}}$$

$$V_p=K_{p}sin\left[2\pi \varphi (t)\right]+v_{sys}+v_{bary}(t).$$

$$K_p=v_{orb}sin(i)$$

$$v_{orb}=\frac{2\pi a}{P_{orb}}$$

$$\varphi (t)=\frac{t-T_0}{P_{orb}}$$

$${\sigma }_{value}=S{F}^{-1}(p\text{-}value/2)$$

$$SF=1-CDF.$$

Full text

Here we present a new automatic pipeline to analyse CRIRES data-sets. Said pipeline is based on a novel use of Principal Component Analysis (PCA) and Cross-Correlation Function (CCF). The exoplanetary atmosphere is modelled with the $\mathcal{T}$-REx code using opacities at high temperature from the ExoMol project. In this work we tested our analysis tools on the detection of CO and H2O in the atmospheres of the hot-Jupiters HD209458b and HD189733b. The results of our pipeline are in agreement with previous results in the literature and other techniques.

methods: data analysis — planets and satellites: atmospheres — planets and satellites: individual (HD189733b and HD209458b) — techniques: spectroscopic

††journal: ApJ

\turnoffedit\turnoffediting\published

June 25, 2019

1 Introduction

More than 4000 confirmed exoplanets are currently listed in the catalogues, together with basic planetary, stellar and orbital parameters as they become known. Transit and direct imaging spectroscopy from space and ground facilities are enabling the study of the physical and chemical properties of some of these exoplanets. From space, one can observe exoplanet spectra in the UV, VIS and IR at low spectral resolution, without the hurdle of telluric contamination. Molecules, ions, atoms or absorbers able to imprint strong modulations in the recorded spectra, can be detected by using space-borne facilities, (e.g. Charbonneau et al. (2002); Linsky et al. (2010); Tinetti et al. (2007); Grillmair et al. (2008); Sing et al. (2016); Fraine et al. (2014); Damiano et al. (2017); Tsiaras et al. (2016a, b, 2018)). By contrast, observations from the ground at high-resolution (R$>$25,000) have enabled the detection of molecules or atoms whose weak absorptions are hard to detect at low spectral resolution. This is particularly true for alkali metals and CO which have been found in the atmospheres of most hot-Jupiters analysed (Redfield et al., 2008; Birkby et al., 2013, 2017; Birkby, 2018; Brogi et al., 2014, 2016; de Kok et al., 2013; Snellen et al., 2010).

High-resolution spectroscopy (HRS) allows to resolve molecular bands into individual lines. Using radial velocity measurements and techniques such as Cross-Correlation Function (CCF) we may separate three physically different sources: telluric absorption, stellar absorption and the planetary spectrum, which are normally entangled. The aim – but also the biggest challenge – is to recognise the planetary signal among the telluric and the stellar signals, which can be orders of magnitude stronger. The standard method used in the literature to analyse HRS data is to apply a number of corrections which involve the correction of the airmass, the subtraction of a modelled stellar spectrum from the data and the use of ad-hoc masks to eliminate residual strong features (Birkby et al., 2013, 2017; Birkby, 2018; Brogi et al., 2014, 2016; Snellen et al., 2010).

In this paper we present and assess an alternative automatic procedure to analyse HRS data , which requires no manual intervention that could interfere with the objectivity and repeatability of the analysis. Our analysis method is based on a novel use of Principal Component Analysis (PCA) and Cross-Correlation Function (CCF). The exoplanetary atmosphere has been simulated using $\mathcal{T}$-REx (Waldmann et al., 2015b, a) and line lists have been adopted from the ExoMol project (Tennyson et al., 2016).

We applied our analysis method to two datasets recorded with VLT/CRIRES freely available on the ESO archive. The exoplanets observed are HD209458b and HD189733b (see Tab. 1), the most studied planets up to date, and therefore good examples for testing new and/or different data analysis techniques. HD209458b (Mazeh et al., 2000) was the first planet analysed with high-resolution spectroscopy: Snellen et al. (2010) reported a detection of CO in its atmosphere. CO is absent in the Earth’s atmosphere but also in the stellar spectrum due to the relatively hot temperature of HD209458. The CO signal in the exoplanetary atmosphere should not be contaminated by the star and Earth’s atmosphere. By contrast the star hosting HD189733b is a K-type (Bouchy et al., 2005) showing CO absorption features in its spectrum: additional caution is therefore needed to remove the potential stellar contamination. Brogi et al. (2016) have reported the detection of H2O and CO in the atmosphere of HD189733b.

In Section 2 we describe our analysis method, in Section 3 we show the results and in Section 4 discussion and conclusions are presented.

2 Data Analysis

We selected datasets relative to HD189733b and HD209458b which are publicly available on the ESO archive. These are part of 289.C-5030(A) and 383.C-0045(A) programs (PI Snellen, I.) (Fig. 1 and 2 panels (a)). The observations have been recorded by using VLT/CRIRES at the highest resolution available (R$=100,000$) through the $0^{\prime\prime}.2$ slit. Both datasets cover a narrow wavelength range, i.e. 2287.54-2345.34\nm and 2291.79-2349.25\nm respectively, with three gaps ($\sim 200$ pixels per gap) due to the physical separation of CRIRES’ detectors. Both datasets have been recorded with the nodding method ABBA for a better background subtraction (Snellen et al., 2010; Brogi et al., 2016). The steps of the analysis process are represented in Fig. 3 and they are described in following sections.

2.1 Data reduction and calibration

We adopted the pipeline provided by ESO (Crire kit Version-2.3.3) to process the raw data. The CRIRES’ reduction pipeline has been embedded into our code thanks to ESO’s EsoRex which is a command-line driven utility that can launch pipeline reduction routines (they are referred as recipes). These are individual scripts that perform specific actions to the input data. The reduction process performs the following steps:

• •

dark subtraction;

• •

correction for detector non-linearity;

• •

flat-fielding;

• •

combination of nodding exposures;

• •

spectrum extraction;

• •

wavelength calibration.

The master reduction files (e.g. dark and flat) are provided with the raw data, while the specific non-linearity correction files need to be downloaded from the archive111https://www.eso.org/sci/facilities/paranal/instruments/crires/doc/VLT-MAN-ESO-14200-4032_v91.pdf. The 1D spectrum is extracted from the reduced images via an optimal extraction (Horne, 1986). By using the ABBA nodding method, we obtained 45 spectra for HD189733b and 51 for HD209458b dataset.

To subtract and correct the telluric absorbtion, the calibration from the ESO pipeline is not accurate enough, we followed instead the procedure described in the literature (Birkby et al., 2013, 2017; Brogi et al., 2013, 2014, 2016; de Kok et al., 2013; Snellen et al., 2010), which involves a further calibration using the SKYCALC tool 222https://www.eso.org/observing/etc/bin/gen/form?INS.MODE=swspectr+INS.NAME=SKYCALC. This simulates the telluric absorption spectrum for a specific night.

The first step is to normalise each spectrum of each detector by dividing it by its median. This step is necessary to avoid differences of baseline across spectra. After the normalisation, working on one detector at a time, we consider the mean spectrum. Here, the strongest lines (all the lines reaching a minimum $<0.8$) have been identified as homogeneously distributed as possible to cover the whole x-axis range. These same lines are also been identified within the telluric template. A Gaussian fit is then performed for each of these lines and the centroid is taken. The extracted spectrum centroids indicate the pixel number position. In the telluric template, instead, they indicate wavelength positions of the lines. We performed a fourth order polynomial fit to establish the relationship between pixels and wavelengths (Snellen et al., 2010). All the single spectra are then interpolated via a third order spline to the derived wavelength grid to have the same grid for all the spectra.

We analysed each detector separately as a two-dimensional matrix, where the x-axis contains wavelengths and y-axis time: every row of this matrix is a spectrum, every column is a temporal-series at a given wavelength (see Fig. 1 and 2 panels (a)). We have therefore four different matrices. Finally, the pipeline removes all the cosmic rays or spikes that could occur at the edges of the spectra due to the spline interpolation to the wavelength grid. The pipeline takes one column at a time of each 2D matrix, it calculates the median of the column and all the values outside $3\sigma$ from the median are set to the median value.

2.2 Decomposition Analysis (PCA)

The next steps involve the correction for telluric absorption, the subtraction of stellar signal and subtraction of correlated noise.

More recently, the algorithm SYSREM developed by Tamuz et al. (2005) has been adopted to perform a similar task (Birkby et al., 2017; Nugroho et al., 2017). SYSREM allows to extract components iteratively one by one, however, the orthogonality of the extracted components is not guaranteed (Tamuz et al., 2005).

As PCA is highly sensitive to data scaling, we subtracted each column of the data matrices by its mean (Fig. 1 and 2 panels (b)). On a typical spectroscopic dataset, the number of spectra are less than the wavelength bins, resulting in matrices that have more columns than rows. Here we adopt the eigenvalue decomposition (EVD) of the covariance matrix (Jolliffe, 2002). The dimension of the covariance matrix and the number of principal components (eigenvectors) are equal to the number of rows of the input matrix. Two cases are then considered:

• •

time domain matrix (TDM); we use the individual spectra as rows and the wavelength bins as columns.

• •

wavelength domain matrix (WDM); we transpose the matrix to have the spectra as columns and wavelength bins as rows;

In the WDM/TDM case the principal components (eigenvectors) contain the information of the correlations in the wavelength/time domain. We consider, for example, the first detector of the HD189733b dataset: Fig. 4 shows the first five components of the TDM case (left) and the first five of the WDM (right). The TDM components contain the time-domain information and the first one, in particular, is linked to the variations of the airmass: these are linearly correlated as we can appreciate from Fig. 5. The WDM components show the correlation in the wavelength domain and they appear to be correlated with the telluric transmission spectrum. A good example is the strong feature around $200$ (Fig. 4 x-axis unit, $\sim$2290 nm) that persists in all the components.

The TDM case has been chosen as best method for the following reasons:

• •

the WDM component space cannot be fully described since there are more variables (1024 spectral bins) than observations (45 spectra for HD189733b and 51 for HD209458b dataset). The eigenvalues are null after the 44th or 50th component, depending on the dataset;

• •

the application of a telluric mask is required if the WDM case is chosen to remove most prominent telluric features that persist after PCA has been applied.

Following this choice, we calculated 50 TDM components for the HD209458b dataset and 44 components for the HD189733b dataset. They are equal to the number of recorded spectra minus one, due to the normalisation performed before the PCA decomposition. From the eigenvalues we calculated the explained variance ratio (EVR) as follows: EVR${}_{j}=\lambda_{j}/\sum\lambda_{i}$, where $\lambda_{i}$ are the eigenvalues. The EVR estimates the information carried by each principal component in percentage. The EVRs of each principal component for every detector of both datasets are shown in Fig. 6 and 7. The first component has always the largest variance ($\sim$ 80%) as the telluric signal is the most significant.

After the application of PCA, each column of the output matrix was divided by its standard deviation to restore the S/N of the processed data (de Kok et al., 2013; Birkby et al., 2013; Ridden-Harper et al., 2016; Nugroho et al., 2017) (Fig. 1 and 2 panels (c)).

2.3 Cross-Correlation Function (CCF)

The cross-correlation function measures the similarity of two signals. It is also often called sliding dot product since it returns a single value from the product of two signals when one slides over the other. Considering two series x and y, the normalised cross-correlation $CCF$ at the delay $d$, for discrete series, is defined as follows (Bracewell, 1965)

[TABLE]

where $\overline{x}$ is the mean of the array x, $\overline{y}$ is the mean of the array y and $i=0,1,2...N-1$. The idea of using such function is to find possible correlations between the data and an atmospheric model. The cross-correlation aims at matching similarities between the two signals.

The exoplanet atmospheric models have been simulated using $\mathcal{T}$-REx  (Waldmann et al., 2015b, a). The CO and H2O line lists at the planetary temperature were provided by ExoMol (Tennyson & Yurchenko, 2012; Tennyson et al., 2016).

Every row of the data matrix (every single spectrum), after the application of PCA, is cross-correlated with the simulated exoplanet atmospheric spectrum. This spectrum is interpolated to the same wavelength grid of the data, and it is then shifted from $-100$ to $100$ km/s with $1.0$ km/s as step. The step is chosen based on the precision obtained during the calibration step ($\sim$ 1.0 km/s) and on the velocity resolution of the instrument (1.5 km/s).

The CCF transforms the matrices (one for each of the four CRIRES’ detectors) from the wavelength domain to the velocity domain. The CCF matrices are then added together to obtain one single matrix (we will refer to it as CCF matrix).

At this stage the exoplanetary signal is not visible (see Fig. 8 top left panel). We then injected a synthetic signal with the orbital parameters of the planet (Fig. 8 bottom left panel) to predict the position of the signal (Fig. 8 bottom left panel) and to calculate the area of the S/N matrix interested by the planetary signal (see Sec. 2.5 and Fig. 9, 10). This step is performed after the calibration (before the PCA decomposition) and the effects of this process are not visible until the CCF is performed, because the signal intensity is at least three order of magnitude weaker than the telluric and stellar signals. The injected model cross-correlates with itself resulting, for example, in the signal shown in Fig. 8 bottom panels. Finally, the injection process allows us to monitor the signal during the PCA decomposition. That helps to determine when the component subtraction starts to erase part of the signal.

2.4 Signal extraction

v

At this stage the planetary signal was barely visible or completely invisible, therefore we co-added the single CCFs (rows of the CCF matrix) in-transit to obtain the integrated signal from the planet. As the data were aligned to the telluric spectrum reference system, the planetary spectrum moved across time, we then re-aligned the single CCFs to the reference system of the planet by computing the following correction:

[TABLE]

where $K_{p}$ is the radial velocity amplitude of the planet (Eq. 3) and $\phi(t)$ is the orbital phase (Eq. 5)

[TABLE]

[TABLE]

[TABLE]

All parameters are listed in Table 1.

Once all CCFs were aligned to the planetary rest frame, we co-added in time only the in-transit CCFs. These were selected by computing the transit time (Seager & Mallén-Ornelas, 2003; Kipping, 2010). When all the in-transit cross-correlations are summed together, the 2D cross-correlation matrix is reduced to a 1D signal, which is connected to the theoretical orbital velocity of the planet. To explore different orbital velocities we proceeded as follows:

we let $K_{p}$ varying from [math] to $250$ km/s with $1$ km/s step; 2. 2.

for each $K_{p}$ we applied the correction in Eq. 2 to every single CCF in the CCF matrix; 3. 3.

we summed only the in-transit cross-correlations.

In this way we were able to explore all possible orbital velocities including those corresponding to the host star. Following the previous steps, we obtained a matrix with $K_{p}$ on y-axis and velocity rest frame along x-axis ($v_{rest}$). From this matrix two different outputs were extracted: the S/N map and the T-test statistic.

2.5 S/N matrix

We considered the last matrix obtained, i.e. $K_{p}$ on y-axis and $v_{rest}$ on x-axis. We calculated the standard deviation of this matrix excluding those points potentially correlated to the planetary signal (|v_{rest}|\ <\ 15\km/s) and we divided the entire matrix by this value. We refer to the obtained matrix as the S/N matrix (Fig. 11 and 12 left panels).

To assign an uncertainty to the $K_{p}$ value we followed the same procedure as reported in Brogi et al. (2016): i.e. we took the maximum value of the matrix and, fixing the relative $v_{rest}$, we calculated the $K_{p}$ interval where the S/N dropped by a unit around the $K_{p}$ peak. The same approach was used to determine the uncertainty for the $v_{rest}$.

The S/N map is not only useful to represent visually the results but also to inspect whether spurious signals or telluric residuals are present. These signals may have high S/N value but are located at different $K_{p}$ and/or $v_{rest}$ from those expected for the planetary signal.

We calculated the S/N matrix for each excluded principal component. Two loops need to be performed to explore the entire principal component space: the first loop subtracts higher variance components onwards and aims to remove the most correlated signal (e.g. telluric absorption and stellar signal). The second loop subtracts lower variance components backwards and aims to remove uncorrelated noise from the data. Finally, the principal components were selected to maximise the peak of the S/N matrix in correspondence of the expected planetary $K_{p}$ and $v_{rest}$.

2.6 Welch’ T-test statistics

The Welch’s T-test is used to test the hypothesis that two populations have equal means. This test compares the population of points on the CCF map connected to the planetary signal with those that are not.

From the CCF matrix we defined, as done in the literature (Brogi et al., 2016; Nugroho et al., 2017):

• •

in-trail, those values inside a squared box centred on the CCF’ peak with a radius of $\pm 15$ km/s;

• •

out-trail, those values outside the in-trail box

We extracted two families of values from the CCF matrix and these were compared through Welch’s T-test (Fig. 11 and 12 right panels). The test, calculated using scipy.stats.ttest_ind in python, provides a p-value (two-tailed) which was converted into $\sigma$-value (significance interval) through the inversion of the survival function (SF)

[TABLE]

where the $SF^{-1}$ is the inverse of the survival function that is calculated from the cumulative density function (CDF) as follows:

[TABLE]

3 Results

HD209458b

for the cross-correlation process, we assumed an isothermal $T-p$ profile at $T=1400$K, with the pressure varying from $10^{-5}$ to $10^{4}$ Pa. We did not include clouds nor line-broadening due to the rotation of the planet. We used $10^{-3}$ as Volume Mixing Ratio (VMR) for both molecules, this value is compatible with chemical models’ predictions for Hot-Jupiters atmospheres (Venot et al., 2012). The same value was also used by Snellen et al. (2010) for the CO.

The signal obtained for CO peaks at S/N=5.7 (Fig. 11 top left panel and Tab 2). The signal is compatible with the planetary orbital parameters ($K_{p}=148_{-15}^{+16}$ km/s, $v_{rest}=-3.0_{-1.1}^{+1.3}$ km/s). This result has been obtained by considering components from the 7$th$ to the 28$th$ . The statistical significance of the result is also confirmed by the Welch’s T-Test (Fig. 11 top right panel). Using a box of radius 15 km/s the null hypothesis is rejected with a confidence greater than 7$\sigma$, the shift of the in-trail population is noticeable with respect to the out-trail values that are, instead, distributed as a Gaussian centred to zero.

The signal of the water vapour is more difficult to detect since the Earth’s atmosphere also contains water. To extract the planetary signal a robust telluric correction is required, and therefore several components need to be subtracted using PCA. A signal at the compatible planetary parameters is observable in the S/N map in Fig. 11 bottom left panel. The maximum peaks at S/N=3.95, $K_{p}=140_{-16}^{+25}$ km/s and $v_{rest}=-4.0_{-1.6}^{+1.4}$ km/s and it is obtained considering components from the 33$th$ to the 43$th$ . To demonstrate that the H2O planetary signal survives after 33 components have been subtracted, Fig. 9 shows the in-transit co-added cross-correlation relative to the range of components aforementioned. Both the injected and non-injected signal survive to the PCA correction (note that the injected signal does not include any atmospheric dynamics, so it is not blue-shifted as the planetary signal). Moreover, the co-added cross-correlation value is lower with respect to the CO case meaning that the concentration of water is lower than CO or that PCA has erased part of the signal. Finally, the Welch’s T-Test is performed on the in-trail and out-trail populations (Fig. 11 bottom right panel and Tab 2). In this case the shift of the in-trail population is not as strong as in the CO case but the null hypothesis is rejected with a confidence greater than 6$\sigma$.

HD189733b

the planetary transmission spectrum was modelled with an isothermal $T-p$ profiles at $T=1000$K. The pressure varies from $10^{-5}$ to $10^{4}$ Pa and we did not include clouds or any lines broadening due to the rotation of the planet. We used $10^{-3}$ as Volume Mixing Ratio (VMR), this value is compatible with chemical models’ predictions for Hot-Jupiters (Venot et al., 2012).

The CO detection is highly difficult since the star, being a K-type star (T$\sim 4900$ K), contains CO in the outer regions. In Brogi et al. (2016) a master stellar spectrum has been simulated and subtracted to the data, but the stellar contamination continued to be persistent also in the result. In this work, PCA was not as effective as in the HD209458b case because the star spectrum moves 1-2 pixels on the detector preventing an optimal correction. The result (see Fig. 12 top left panel and Tab. 2) is compatible with the one claimed by Brogi et al. (2016) (S/N=5.1, $K_{p}=194_{-41}^{+19}$ km/s, $v_{rest}=-1.7_{1.2}^{1.1}$ km/s), however the error on the $K_{p}$, being smaller than the one reported in literature, does not include the theoretical value of the orbital velocity of the planet ($K_{p}=152.564$ kms*-1*). The signal determined at lower $K_{p}$ ($\sim 85$ km/s) is due to stellar contamination, that results from the Rossiter-McLaughlin effect combined with the change of reference frame from the Earth to the barycentric one (Brogi et al., 2016).

Concerning water vapour, the same discussion done for HD209458b can be applied here. The planetary water signal needs to be disentangled from the telluric absorption. The result obtained is compatible with both the literature and the theoretical parameters, e.g. see Fig. 10 where the planetary signal is compared with the injected one. The injected signals do not account for $v_{rest}\neq 0$ km/s, we can appreciate the data being blue-shifted. Here the Welch’s T-Test confirms that the null hypothesis can be rejected with a confidence greater than 5$\sigma$ (Fig. 12 bottom right panel and Tab 2).

We performed an additional test by cross-correlating the telluric model used in the calibration process with the data, to check if any telluric signal still persists. Using the components reported in the results we did not notice any significant correlation with the telluric signal at the position of the planetary parameters.

4 Discussion and conclusions

We presented here and tested a new automatic method based on the iterative use of PCA and CCF to re-analyse two CRIRES datasets observed with high resolution spectroscopy technique.

CO and H2O have been detected in the HD209458b dataset, and H2O in the HD189733b dataset. The detection of CO in the HD209458b atmosphere is supported by an S/N peak of 5.7 at $K_{p}$ and $v_{rest}$ compatible with the planetary orbital parameters. Contrary to CO, H2O is present in the Earth’s atmosphere and therefore an accurate telluric correction is required. The lower S/N peak may be due to a lower concentration of H2O with respect to CO in the atmosphere of HD209458b, or part of the signal might have been removed by PCA. In both detections a blueshift has been observed and this could be explained with high altitude winds. The results presented here are in agreement with the results published by Snellen et al. (2010).

Concerning HD189733b, using our method, we have been able to detect H2O. Even in this dataset a blueshift of the signal has been observed and also in this case it could be associated with high altitude winds. The detected CO signal, is compatible with the literature (Brogi et al., 2016) but it is not in agreement with the theoretical radial velocity of the planet, and could be due to stellar contamination (K-type star ).

We note that the EVR is different for each detector (Fig. 7 and 6) and this means that the planetary signal is contained in different components in each detector. An optimal approach should adapt the number of components per detector based on their variance.

Future work will consider the use of the algorithm presented here to analyse high-resolution observations taken by other instruments. These include CRIRES+ (Follert et al., 2014), GIANO-B, a high-dispersion spectrograph at TNG (Oliva et al., 2012) which covers $0.9$ to $2.5\ \mu$m with a resolution of (R=50,000); IRCS-SUBARU (Kobayashi et al., 2000), which uses a lower resolution (R=20,000) but covers a broader range (from $1$ to $5\ \mu$m) and CARMENES at Calar Alto Observatory (Quirrenbach et al., 2014) with a spectral resolution up to 80,000 in the near-IR ($0.9-1.7\ \mu$m).

Acknowledgments

M.D and G.T. are supported by the European Research Council (ERC) grant 617119 ExoLights. M.D. and G.M. are supported by the project ”A Way to Other Worlds” funded by Italian Minister for University and Research. This work was also supported by STFC (ST/P000282/1). The authors wish to thank Dr. Antonino Petralia, Dr. Ingo Waldmann, Prof. Jonathan Tennyson, and Dr. Sergey Yurchenko for their input.