© D. Shulyak et al. 2019

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1 Introduction

Studying exoplanetary atmospheres is a key to understanding physico-chemical processes, origin and evolution paths, and habitability conditions of these distant worlds. This is usually done with the help of spectroscopic and/or (spectro-)photometric observations from space and from the ground. However, obtaining radiation from an exoplanetary atmosphere at an accuracy sufficient for its detailed investigation is a difficult task due to the large distance to even nearby exoplanets, which results in a very weak photon flux reaching the top of the Earth’s atmosphere. Photometric observations are the more efficient than spectroscopic observations in terms of the time needed to achieve a required signal-to-noise ratio (S/N) for the exoplanetary signal. By using photometric observations obtained in a wide wavelength range (which may cover dozens of microns in the infrared domain) it is possible to constrain the temperature stratification and mixing ratios of most abundant molecules in the atmospheres of exoplanets. So far, the atmospheres of hot Jupiters (HJ) have been a subject of intense investigations (e.g., Brogi & Line 2019; Sing et al. 2016; Waldmann et al. 2015a; Schwarz et al. 2015; Line et al. 2013; Lee et al. 2012). These gas giants have masses between 0.5 and 13 M_J and short rotation periods P < 10 days (Wang et al. 2015) so that their atmospheres are hot due to a close distance to their parent stars (Winn et al. 2010). This results in strong irradiation from their atmospheres that is much more easy to detect than the atmospheres of similar exoplanets that are at larger distances from their stars. To date, there are about a dozen HJs that have been investigated using photometric observationsfrom different space missions (e.g., Sing et al. 2016).

Spectroscopic techniques, on the other hand, are capable of resolving individual spectral lines of atmospheric molecules, but require much longerobserving times and are usually limited to a narrow wavelength range compared to photometry. Nevertheless, modern high-resolution spectroscopy has been successfully used to detect the presence of many molecules in atmospheres of different types of exoplanets via cross-correlation techniques (e.g., Snellen et al. 2010; de Kok et al. 2014; Brogi et al. 2014, 2016; Kreidberg et al. 2015). Unfortunately, extracting profiles of individual molecular lines from available spectroscopic observations remains very challenging due to the high noise level, telluric removal problems, and purely instrumental effects. Recently, Brogi & Line (2019) showed that the cross-correlation technique can be used to asses the temperature structure in the atmospheres of HJs via the differential analysis of weak and strong features in the retrieved exoplanet spectra without the need to resolve profiles of individual lines (see, e.g., Schwarz et al. 2015). However, only very limited information can be retrieved and it is still not possible to obtain robust results, for example due to the quality of the observed spectrum. This is why most past and present studies utilize spectroscopy to detect the presence of molecules, but they do not attempt to study atmospheric temperature stratification and accurate mixing ratios.

The advantage of spectroscopy in providing accurate estimates of molecular mixing ratios is currently limited by the capabilities of the available instruments. In particular, the quality of the observed spectrum is often not good enough to attempt studying atmospheric temperature stratification and accurate mixing ratios. However, with the advent of future instrumentation it will be possible to study the accurate shapes of atmospheric lines on diverse exoplanetary atmospheres. This will open a way to constrain atmospheric chemistry, global circulation (winds), and signatures of trace gases produced by possible biological activity.

Motivated by recent progresses in exoplanetary science and advances in instrumentation, in this work we investigate the potential of applying very high-resolution spectroscopy to the study of HJ atmospheres. We focus our research on the CRyogenic high-resolution IR Echelle Spectrometer (CRIRES+) on the Very Large Telescope scheduled for the end of 2019 at the European Southern Observatory (ESO) (Dorn et al. 2014; Follert et al. 2014). CRIRES+ is an upgrade of the well- known old CRIRES; it has improved efficiency thanks to new detectors, polarimetric capability, and an order of magnitude larger wavelength coverage that can be achieved in single-setting observation thanks to a new cross-disperser. In particular, this improvement is essential for exoplanetary research because it opens a way to observe many molecular features and to study atmospheric physics in great detail. Similar to the previous version of the instrument, CRIRES+ will operate between 0.9 and 5.3 μm with a highest achievable resolving power of λ∕Δλ = R = 100 000¹.

Our main goal was to test the near-infrared wavelength domain at very high spectral resolution for the retrieval of atmospheric temperature profiles and mixing ratios of molecular species using CRIRES+ as a test case. The key difference between our research and similar projects is that our aim was to retrieve assumption-free temperature profiles and that we used different near-infrared spectral regions that contain numerous lines of atmospheric molecules expected to be observed with CRIRES+. We simulated observations at different wavelength regions and their combinations, spectral resolutions, and S/N to test our ability to derive accurate temperatures and chemical compositions of HJ atmospheres. We made predictions to identify observational strategies and potential targets to study exoplanet atmospheres with high-resolution spectroscopy. Finally, we note that our results can also be applied to any other present or future instrument with capabilities similar to CRIRES+.

2 Methods

2.1 Estimating the state of the atmosphere

In our retrieval approach we used an algorithm based on optimal estimation (OE; Rodgers 1976). The main idea of OE is to derive a maximum a posteriori information about parameters that define the state of the system under investigation. The OE approach finds an optimal solution that minimizes the difference between the model and observations under an additional constraint, which is our knowledge of a priori parameter values and their errors. Following Rodgers (2000) we solve the OE problem for the state vector of the system (1)

where x_i+1 and x_i are respectively the new and old solutions for the state vector at iteration i + 1 and i, x_a and S_a are the a priori vector and its co-variance matrix, y and S_y are the measurements and measurement errors, and y_i and K_i are respectively the model prediction for the state vector x_i and the corresponding Jacobian matrix calculated at x_i. An adjustable parameter (γ) is used to fine-tune the balance between the measurement and the a priori constraint (Rodgers 2000; Irwin et al. 2008). We searched for the solution that optimizes the total cost function ϕ that can be written in the form (2)

where the first term is the usual χ² computed as a weighted difference between the model and observations, while the second term accounts for the deviation of the state vector from its a priori assumption. After the solution has converged, the final errors on the state vector are (3)

It is seen that the obtained solution depends on how accurately we know our a priori. By setting large errors on the a priori, the OE eventually turns to a regular nonlinear χ² -minimization problem. In the case of solar system planets, the a priori values of physical parameters such as temperature stratification and surface pressures can be known from in situ measurements, if available. This helps to substantially narrow down the parameter range and the application of OE provides robust solutions with realistic error estimates (Hartogh et al. 2010; Irwin et al. 2008; Rengel et al. 2008). In the case of extrasolar planets, little to nothing is usually known about their physical and chemical structures. Nevertheless, even in these cases it is possible to build initial constraints upon some simplistic analytical models and to estimate realistic error bars a priori. Because of its relatively fast computational performance compared to other methods, such as the Markov chain Monte Carlo (MCMC) algorithm, the OE method was successfully used to study the atmospheres of exoplanets (see, e.g., Lee et al. 2012, 2014; Line et al. 2012, 2013; Barstow et al. 2013, 2014, 2017).

Like every retrieval technique, OE has its own advantages and weaknesses. An obvious limitation is the use of a priori information, which penalizes the solution toward the one that provides the best balance between the observed and a priori uncertainties. Because the inverse problem we are dealing with is ill posed, there can often be a varietyof solutions that provide the same good fit but very different values of fit parameters. This happens if, for example, some of the parameters are strongly correlated and/or the data quality is poor. In these cases using the a priori constraints helps to keep retrievals from finding the solution with the lowest χ² but not expected from the physical point of view (i.e., strong fluctuations in temperature profile, molecular mixing ratios that are too far from thepredictions of chemical models, among others.). Unfortunately, in the case of exoplanets the use of inappropriate a priori information can lead retrievals to fall into local minima rather than finding the true solutions. Inthe OE method this difficulty is easy to overcome by choosing large enough errors on (poorly known) a priori parameters whose influence on the final solution is thus substantially reduced. Usually, a family of retrievals with different sets of initial guesses and a priori errors are used to ensure that the found solution is robust (see, e.g., Lee et al. 2012).

Another limitation of the OE method is the assumption that the distribution of the a posteriori parameters are Gaussian. This is usually the case for data that is of high spectral resolution and S/N, which is satisfied in most of our retrievals presented below. When the data is of low quality the usual workaround is again to perform retrievals with multiple initial guesses to ensure that the global minimum is found. The relative comparison of different retrieval methods can be found in Line et al. (2013). In our work we do not favor any of the retrieval methods specifically. The choice of the retrieval approach normally depends on the goals of the particular investigation, quality of the available data, number of free parameters, available computing resources, among others. Here we chose the OE method because (1) it has long been used for the analysis of high-resolution spectroscopic data of planets in our solar system, which means that the expertise can be easily shared between research fields; (2) future instruments will eventually provide us with high-resolution data for the brightest exoplanets, and this data will be of better quality compared to what we have now, thus allowing the OE method to use all its advantages; and (3) we carried out assumption-free retrievals of temperature distribution, which means that the local value of temperature in each atmospheric layer is a free parameter in our approach; because the atmosphere is split into several tens of layers, this results in too many free parameters to be treated with purely Bayesian approaches.

2.2 Forward model

In order to use the OE method, an appropriate radiative transfer forward model is needed that computes the outgoing radiation from the visible planetary surface for a given set of free parameters. We did this by employing the Tau Retrieval for Exoplanets (τ-REx) software package (Waldmann et al. 2015b,a). The τ-REx forward model uses up-to-date molecular cross sections based on line lists provided by the EXOMOL² project (Tennyson & Yurchenko 2012) and HITEMP (Rothman et al. 2010). These cross sections are pre-computed on a grid of temperature and pressure pairs and are stored in binary opacity tables that are available for a number of spectral resolutions. The continuum opacity includes Rayleigh scattering on molecules and collisionaly induced absorption due to H₂ -H₂ and H₂-He either after Abel et al. (2011, 2012) or Borysow et al. (2001); Borysow (2002) and Borysow & Frommhold (1989), respectively.The original line opacity tables contain high-resolution cross sections for each molecular species and for a set of temperature and pressure pairs. Thus, we optimized the public version of τ-REx³ to compute spectra in wide wavelength intervals and with very high spectral resolution R = 100 000 expected for CRIRES+ via the more efficient usage of opacity tables inside the code. Finally, we adapted our retrieval code for multiprocessing which is essential for the analysis of altitude dependent quantities (e.g., T, mixing ratios, winds) when the number of free parameters can be very high.

In order to keep solutions stable and to ensure that the retrieved altitude dependent quantities (e.g., temperature) are smooth, it is usually assumed that these quantities are correlated with a characteristic correlation length l_corr expressed inpressure scale heights, for example. Then, the off-diagonal elements of the a priori co-variance matrix can be written as (4)

where p_i and p_j are the pressures at atmospheric layers i and j, respectively, and l_corr is the number of pressure scale heights within which the correlation between layers drops by a factor of e. We found that a value of l_corr = 1.5 provides a necessary amount of smoothing without losing the physical information that we retrieve, and thus we adopted this value in all our retrievals (the same value of l_corr was also used in Irwin et al. 2008).

If no correlation between atmospheric layers is used then the inverse of co-variance matrices can be analytically computed. To the contrary, when the correlation length is non-zero the S_a matrix contains off-diagonal elements whose values can differ by many orders of magnitude which leads rounding-off errors in the numerical matrix inversion. In order to reduce numerical problems we normalize the state vector co-variance matrix S_a by its diagonal elements and re-normalize obtained solution in Eq. (2) accordingly.

In our retrievals we make no assumptions about the altitude dependence of temperature profile. However, we assumed constant mixing ratios of molecular species. We do this to be consistent with the original work by Lee et al. (2012) (see next section). This also decreases the number of free parameters in our model and significantly reduces calculation time.

The final set of free parameters in our retrieval model includes 50 temperature values corresponding to different altitude levels, four molecular species (H₂ O, CO, CO₂, CH₄ ), and a continuum scaling factor for each wavelength region that we investigated. The scaling factors account for any possible normalization inaccuracies between observed and predicted spectra. This adjustment will be needed because after the reduction each CRIRES+ spectra will be normalized to an arbitrary chosen level (e.g., pseudo continuum). In this work we make no attempts to simulate these very complicated reduction steps that should be a subject of a separate investigation. Instead, our simulations represent an ideal case where we do not consider all instrument systematics and data reduction inaccuracies. Thus, we normalized the observed spectrum to some arbitrary value (mean flux level in our case) and did the same with the fluxes computed with our forward model. The additional scaling factor is then applied to our forward model to account for any remaining mismatch between the observed and modeled spectra. Clearly, these scaling factors should be very close to unity in the case of accurate retrieval because the observed and predicted spectra were computed with the same code. Still, these additional scaling factors will be needed in retrievals performed on real dataand we include them in our model.

In this work we concentrate only on emission spectroscopy. We do this because emission spectra are formed in a very wide range of atmospheric temperatures and pressures, while transmission spectroscopy (performed when a planet transits in front of the host star) senses mainly regions of high altitudes. In addition, the duration of transits of all potential targets are shorter than off-transit times when the day side of the planet is visible. This makes it easier to observe emission spectra with a required noise level.

3 Validation of the method

Before addressing the main goal of our work, we validated our method against some representative examples. To do this we chose a test case of a well-studied exoplanet, HD 189733b. This HJ transits a young main sequence star of spectral type K0 at a distance of about 19.3 pc (Bouchy et al. 2005). The atmospheric structure of HD 189733b was studied by different groups and with different methods in the past (see, e.g., Lee et al. 2012; Waldmann et al. 2015a; Brogi & Line 2019).

Figure 1 summarizes the retrieved temperature and mixing ratio of H₂ O, CO, CO₂ , and CH₄ in the atmosphere of HD 189733b using photometric observations obtained with different space missions (see figure legend). We used results of Lee et al. (2012) as our a priori and assumed 100 K uncertainties for our initial temperature profile. We note that here and throughout the paper our a priori is the same as an initial guess in our iterative retrieval approach, although in a more general case these two can differ. The analysis of corresponding averaging kernels for the temperature distribution (second plot in the bottom panel of Fig. 1) shows that the available observations probe temperature structures in a wide range of atmospheric depths between 10⁻⁶ and 1 bar. The averaging kernels describe a response of the retrieval to a small perturbation in temperature at each atmospheric depth (Rodgers 2000), and thus characterize the contribution of each depth to the final result. The zero values for averaging kernels correspond to the case when no information can be retrieved from corresponding depths and they are nonzero otherwise. It is seen that averaging kernels peak at appropriate atmospheric depths between 10⁻⁶ and 1 bar, thus indicating that our retrievals are robust at these altitudes; they approach zero outside that range.

Figure 1 shows that we retrieved parameters in agreement with original work by Lee et al. (2012). We also confirm the high content of water in the atmosphere of HD 198733 b, and the general shape of the temperature distribution. Finally, the concentrations of CO₂ and CH₄ are badly constrained and we could not measure their content with any initial guess assumed, again in agreement with Lee et al. (2012).

To check the robustness of our results we performed four retrievals with different initial temperature profiles to investigate the sensitivity of our solutions to the initial guess. The result are presented in Fig. 2 where we show retrieved temperature profiles assuming four different initial guesses. In all cases we constrain very similar temperatures which confirms the robustness of our retrieval approach. In the lower and upper atmosphere our solution remains unconstrained because the observations do not sense these altitudes, as can be seen from the right panel of Fig. 2 which shows corresponding averaging kernels for each of the retrieved temperature distribution.

Finally, we note that our temperature distribution and that published in Lee et al. (2012) are noticeably different from the temperature found by Waldmann et al. (2015a) (blue dashed line in Fig. 2). The largest deviation is found in the region of the temperature drop around 0.1 bar. This is likely because of the different opacity tables used in the old version of τ-REx (I. Waldmann, priv. comm.). We conclude that our retrieval method is accurate and robust, and that it can be used for atmosphericstudies.

Fig. 1 Top panel: observed day side fluxes of HD 189733b and our best fit model predictions (see legend inthe plot). Bottom panels: retrieved T-P profile along with corresponding averaging kernels and mixing ratios of molecular species. We used the T-P profile of Lee et al. (2012) a priori (shown as red crosses). Our best fit model is shown with a solid blue line and the shaded area represents 1σ error bars. The assumed a priori are shown as a green dashed line (which coincides with the red crosses in this particular case). The same color-coding is used for the four side plots of mixing ratios (on the right).

Open with DEXTER

Fig. 2 Left panel: retrieved temperature distribution assuming four different initial guesses: three homogeneous at 1000, 1500, and 2000 K (vertical dashed lines), and the one from Waldmann et al. (2015a) (blue dashed line). The retrieved best fit profiles are shown with solid lines along with corresponding error bars as shaded areas. The profile from Lee et al. (2012) is shown as a solid red line. Right panels: averaging kernels for the four retrieved temperature distributions (from left toright and from top to bottom, color-coded accordingly).

Open with DEXTER

Fig. 3 Retrieved temperature and mixing ratios from five spectral regions with R = 100 000 and S/N = 10. In each panel, the first plot compares the best fit predicted spectrum (solid red line) with the simulated observations (black). The second plot is the temperature distribution as a function of atmospheric pressure (solid blue line) and with error bars shown as shaded areas. The red crosses and green dashed line are the true solution and initial guess, respectively. The third plot shows averaging kernels derived for the temperature distribution. Here averaging kernels for different atmospheric depths are randomly color-coded for clarity. The next four plots are the values of the retrieved mixing ratios of four molecular species (color-coded as in the first plot). The mixing ratios were assumed to be constant with atmospheric depth, and we show their values on the y-axis.

Open with DEXTER

4 Results

Compared to its old version, the essential new capability of CRIRES+ is the increased wavelength coverage when observing with a single setting. In one shot, the instrument is able to cover a region corresponding to about one spectral band (i.e., ten times more than the old CRIRES). This gives us the possibility to study many spectral lines simultaneously with our retrieval method and to increase the amount of information that we can learn about exo-atmospheres. In what follows we study several cases of simulated observations at different spectral bands, their combinations, spectral resolving powers, and S/N values. We chose to use five different spectral regions that are free from strong atmospherictelluric absorption and that contain molecular bands of main molecules found in atmospheres of HJs: H₂ O, CO, CO₂ , and CH₄ . These regions are: 1.50–1.70, 2.10– 2.28, 2.28– 2.38, 3.80– 4.00, and 4.80– 5.00 μm. Each of these regions can be observed with single CRIRES+ setting; however, some gaps between spectral orders within each setting are expected. For each of these regions we computed the spectra of HD 189733b assuming spectral resolutions of R = 100 000 and R = 50 000 (corresponding to slit widths of 40′′ and 10′′ , respectively) and four S/N values of 5, 10, 25, and 50, respectively. The S/N values are given per spectral dispersion element (1 pixel), which would be obtained in real observations after integrating the spectrum profile along the spatial direction on the CCD. We note that the S/N of the planetary spectrum is obtained from the S/N of the stellar spectrum by taking into account a characteristic planet-to-star flux ratio S∕N(planet) = F_p∕F_s ⋅ S∕N(star). The values of S/N(star) can be directly computed using exposure time calculators, for example, as we discuss below.

In order to obtain accurate results, we performed retrievals with different initial guesses and their uncertainties. We used three starting guesses for molecular volume mixing ratios (VMRs) of ± 0.5 and − 1.0 dex from the true solution. For each of the initial guess we tried uncertainties of 0.5 and 1.0 dex, respectively. By choosing such large uncertainties on the VMRs we avoid the problem of falling into local minima, as was discussed in the previous section. We found that we always converge to nearly the same solution whether 0.5 or 1.0 dex error bars were assumed. However, this was not always the case, especially for data of low S/N, and we discuss these cases below. For the temperature retrievals we tried different uncertainties between 100 and 500 K and found that the value of 200 K allows us to retrieve a temperature structure that is smooth and very close to the true one. Choosing too high temperature errors causes strong nonphysical fluctuations in the finally retrieved temperature distribution and we do not discuss this solution further.

4.1 Single spectral region retrievals

In Fig. 3 we show the results of retrievals from five spectral regions and S/N = 10 (the results for other S/N values are provided in Figs. A.1–A.3). The data was simulated using temperature structure and mixing ratios of HD 189733 b as derived in Lee et al. (2012) (red dashed linesin Fig. 3). In all cases we start from an isothermal temperature distribution with T = 1500 K. In this way we can better study the ability of the EO method to recover the true temperature distribution at different altitudes. The initial values for mixing ratios were taken to be − 0.5 dex from the true ones with the uncertainty of 1 dex. It is seen that the spectra of HJs in the wavelength range covered by CRIRES+ is sensitive to a wide range of atmospheric pressures between 1 and 10⁻⁴ bar. At these altitudes we recover temperatures very close to the true solution. However, everything that is outside this altitude range remains unconstrained due to the lack of information provided by the observed spectra. The typical temperature errors in the pressure range between 1 and 10⁻⁴ bar is 120 K for S/N = 5, 100 K for S/N = 10, 90 K for S/N = 25, and 80 K for S/N = 50. All spectral windows except 4.80− 5.00 μm appear to be relevant to study temperature stratification in atmospheres of HJs, but only two of them (1.50–1.70) and (2.28–2.38) μm, allow us to constrain VMRs of at least some molecules.

As was mentioned in Lee et al. (2012), the degeneracy between fit parameters may lead to incorrectly derived mixing ratios. We find that this is also the case for high-resolution low S/N data. As an example calculation, our Fig. 4 shows retrievals of H₂ O and CO from the (2.28–2.38) μm region using three different initial guesses for the molecular VMRs. For S/N = 5 the retrievals converge to nearly the same solution if the initial values of VMRs are 0.5 and 1.0 dex below the true one. However, when the initial guess is 0.5 dex higher than the true value then the algorithm prefers to increase the local temperature in the atmosphere keeping the VMRs at their high values. The increase in local temperature makes lines of molecules weaker, thus compensating for the high initial VMR. The molecular VMRs and temperature are thus degenerated as all solutions result in the same final normalized cost function (ϕ = 0.73). As the S/N of the data increases, the solutions for temperature and VMRs both tend to converge to the true values, while for S∕N ≤ 10 the final values may deviate from the true ones (but are still within the 3σ error bars) depending on the initial guess.

Finally, degrading spectral resolution by a factor of two from R = 100 000 to R = 50 000 does not significantly affect the results of our temperature retrievals, while the accurate retrievals of molecular VMRs become problematic because the obtained uncertainties are large. An example is shown in Fig. 5 for the region 2.28–2.38 μm and S/N = 10. With the degrading of spectral resolution, the amount of information that we can learn from the observed spectra mostly depends on the number of spectral features resolved. Decreasing the resolving power weakens the depth of spectral lines, and thus directly affects retrieval of mixing ratios. At the same time, the profiles of strong molecular lines are still satisfactorily resolved, which provides enough information for the temperature retrieval. It is thus possible to rebin the obtained high-resolution spectrum to a smaller resolution and boost the S/N without greatly affecting our ability to retrieve temperature stratification. Alternatively, our Fig. 5 reflects the case of two observations with different slit widths, but after reaching the same S/N. This would require less telescope time for a wider slit width to reach the same S/N, which could be an affordable choice to reduce the integration time for a fixed S/N.

	Fig. 4 Retrievals of the H2O and CO mixing ratios from the (2.28–2.38) μm region assuming three different initial guesses (from top to bottom) and four S/N values (from left to right). The green dashed and red dotted lines are the initial guess and true solutions, respectively. The solid blue line and the shaded area are the retrieved value and its 1σ uncertainty, respectively.
Open with DEXTER

	Fig. 5 Comparison of retrieved temperature and mixing ratios of four molecules assuming S/N = 10 and two spectral resolutions of R = 100 000 and R = 50 000 for the top and bottom panels, respectively. The color-coding is the same as in Fig. 3.
Open with DEXTER

	Fig. 6 Retrieved temperature and mixing ratios of four molecular species from a combination of different spectral regions with a reference region 2.28–2.38 μm. The retrievals are shown for the case of S/N = 10. The color-coding is the same as in Fig. 3.
Open with DEXTER

4.2 Retrievals from multiple spectral regions

As a next step we studied the scenario in which the retrievals are done from a combination of different spectral regions. In particular, we investigated combinations of the (2.28–2.38) μm region with the others because, as noted in the previous section, this region is often used for the detection of molecular lines in atmospheres of HJs, and it also provides rich information about the altitude structure of their atmospheres. Figure 6 compares the results of retrievals considering S/N = 10 (the results for other S/N values are provided in Figs. A.4–A.6). As expected, we find that a combination of (1.50–1.70)+(2.28–2.38) μm regions is one of the best in terms of retrieved information. When analyzed separately both these regions already constrain accurate temperatures, and combining them helps to improve our retrievals even further. On the other hand, the combination of (2.28–2.38)+(4.80–5.00) μm gives us the possibility to probe even higher altitudes up to P = 10⁻⁶ bar, but the concentration of CO is then derived less accurately compared to the previous case with S/N = 10, and noticeably worse with S/N = 5 (Fig. A.4). In general, retrievals from two different spectral regions help to constrain molecular number densities better. For the same case of (1.50–1.70)+(2.28–2.38) μm we can now pin down the mixing ratio of CO and CO₂ by a factor of about two compared to the case when either of these regions is used (compare the VMR uncertainties in the top panel of Figs. 3 and 6). The typical temperature errors in the pressure range between 1 and 10⁻⁴ bar is 100 K for S/N = 5, 90 K for S/N = 10, 80 K for S/N = 25, and 70 K for S/N = 50. We conclude that retrievals from two different spectral regions should be preferred in order to constrain molecular number densities as accurately as possible compared to single-band retrievals, but such observations would require two times longer integration times for the same noise levels.

Finally, simultaneous retrievals from three and four spectral windows (with a maximum wavelength coverage of about 1 μm) did not significantly improve our results for the temperature and mixing ratios compared to the retrievals from two spectral ranges. Thus, the only benefit of observing at many different spectral regions is to measure molecular number densities of as many molecules as possible. However, enough lines of at least three molecules that we considered (H₂ O, CO, CO₂) can be captured already with two cross-disperser settings of CRIRES+. Thus, including additional spectral windows is not justified for the case of atmospheric retrievals.

5 Suggested observation strategies and potential targets

Our analysis shows that with a single setting, CRIRES+ covers a spectral range containing enough features for exoplanet atmospheric studies. However, not all the spectral regions that we considered in this study are equally good for the retrieval of temperature and molecular concentrations. In the case of a single observation setting, accurately determining the temperature and mixing ratios depends only on the S/N of the spectrum of an exoplanetary atmosphere. Our simulations showed that under the same observing conditions the retrievals of accurate mixing ratios will generally require higher S/N compared to the case when only temperature stratification is investigated. For instance, with S/N = 5 it is possible to derive accurate temperatures in the line forming region of an HJ atmosphere (P > 1 bar), while realistic values of the H₂O mixing ratio, for example, can only be obtained with two times higher S/N (see the results for the (2.28–2.38) μm region in Figs. 3 and A.1).

In order to retrieve accurate temperature stratification and mixing ratios of H₂ O, CO, and CO₂, the simultaneous observations at two spectral windows (1.50–1.70)+(2.28–2.38) μm and S/N = 10 should be used. As a next choice the (2.10–2.28)+(2.28–2.38) μm or (2.28–2.38)+(4.80–5.00) μm regions can also be used for robust retrievals. We note that the M band contains many more telluric lines and the performance of the old CRIRES in that region was noticeably worse compared to the performancein K band.

When only S/N = 5 could be reached, the choice of a single observation setting is not obvious, and there are two options. On the one hand, using (1.50–1.70) μm region could provide VMRs of H₂O, CO, and CO₂, but accurate temperatures could only be retrieved up to 10⁻³ bar. On the other hand, the (2.28–2.38) μm region would still provide accurate measurements of H₂O and temperatures up to 10⁻⁴ bar, but less accurate CO and no constraints for CO₂ (see Fig. A.1). Furthermore, the planet-to-star flux contrast in the (1.50–1.70) μm wavelengths is expected to be lower in HJs compared to the (2.28–2.38) μm region, thus requiring longer integration times for the same S/N.

Observing in two spectral regions with S/N = 5 in each of them, we still favor the combination of (1.50–1.70)+ (2.28–2.38) μm. However, for the same integration time, we suggest observing only the (2.28–2.38) μm region obtaining higher S/N = , which allows more precise measurements of H₂O and CO.

In order to provide estimates for the exposure times needed to study atmospheres of HJs we used predictions of our models and the ESO Exposure Time Calculator (ETC) for CRIRES⁴ . Because there is no ETC for CRIRES+ available yet, we notethat the listed exposure times are upper limits due to the expected improved performance of the new instrument compared to the old one. From the list of known HJs we selected several of the best targets that can be observed with high-resolution ground-based spectroscopy similar to the capabilities of CRIRES+. We additionally restricted the list of objects by including only those that require less than 100 h of integration time. We list these targets in Table 1. The integration times were computed using stellar and planetary parameters listed in the current version of the Extrasolar Planets Encyclopaedia⁵ and The Exoplanet Orbit Database⁶. The stellar magnitudes in the V and K bands were extracted from the SIMBAD database⁷. The planetary equilibrium temperatures were computed assuming atmospheric albedo a = 0.03 expected for HJs (Sudarsky et al. 2000). We calculated planet-to-star flux ratios as an average between 1 and 5 μm assuming a blackbodyapproximation for stellar and planetary flux. This values will be larger in the red part of the spectrum (typically by a factor of two), and thus the corresponding exposure times that we list in the table will be shorter for the same S/N at the red end of CRIRES+ spectra. However, the performance of old CRIRES detectors in the M band was noticeably lower compared to the K band so that the estimation of realistic exposure times for CRIRES+ at different spectral bands is currently not possible. This is why we provide exposure times for the mean flux contrast and use one of the most efficient instrument settings at the echelle order No. 25 (λ_eff = 2.267 μm). We note that the S/N for the planetary spectrum is defined by the planet-to-star flux ratio S∕N(planet) = F_p∕F_s ⋅ S∕N(star).

Due to the demand for very high-quality exoplanetary spectra the number of potential targets is very small. Among them, it will be possible to observe only two HJs from the ESO Paranal observatory and with affordable integration times. These exoplanets are 51 Peg b and τ Boo b. For 51 Peg b, a little more than one full night (e.g., 12 h) will be needed to reach S/N = 5 using two instrument settings or S/N = 7 in one setting. This will be sufficient for accurate atmospheric inversions. For our preferable configuration of S/N = 10 and two settings a 48 h integration time will be needed. In case of τ Boo b two full nights will be needed to reach S/N = 5 in two settings, or eight nights for S/N = 10 and two settings. The choice of the corresponding instrument settings will depend on the performance of the CRIRES+ at different infrared bands. Our simulations favors the (1.50–1.70)+(2.28–2.38) μm region because the VMR of three molecules can be accurately constrained, as can the temperature distribution, until P < 10⁻⁴ bar. If one wants to look at higher altitudes then the (2.28–2.38)+(4.80–5.00) μm region should be used, but the uncertainties on CO will be slightly larger. In the case of a single setting we suggest the (2.28–2.38) μm region; however, the (1.50–1.70) μm region can also be an option if one wants to constrain the VMRs of as many molecules as possible. Should the sensitivity of CRIRES+ be considerably improved, then the integration times listed in Table 1 would be smaller, thus favoring a detailed analysis of the atmosphere of 51 Peg b with two setting observations.

6 Discussion

In our analysis of high-resolution simulated observations we used the OE method for the retrieval of atmospheric parameters, such as temperature and molecular mixing ratios. We did not attempt to utilize other methods (e.g., MCMC) primarily because the number of free parameters in our model is high. In this particular case the OE method has the advantage of being relatively fast but still robust. It should be noted that retrieving a detailed structure of planetary atmospheres relies primarily on the data quality and the presence of lines that probe different atmospheric altitudes. In cases when the temperature distribution is expected to be very smooth, or when the quality of the spectrum is poor so that the altitude structure could not be satisfactorily resolved, it is possible to parameterize altitude-dependent quantities with simple functions containing only few free parameters. For example, we could have parameterized the temperature distribution using a commonly used profile after,Parmentier & Guillot (2014), among others, which requires only four free parameters instead of retrieving a temperature value at each atmospheric layer. This would allow us to use a standard Bayesian parameter search. However, in this work we decided to carry out assumption-free retrievals. Hence, our application of the OE method should not be considered as the only suggested approach to the atmospheric retrievals at high spectral resolution, but rather as one possibility among others, with the ability to retrieve accurate parameters with a proper application.

Near-future instruments will deliver data of much better quality and wider spectral range coverage compared to what is available now. That means that at high spectral resolution the shapes of spectroscopic lines originating from the planetary atmospherescould be studied directly. This was not possible to achieve before, and the high-resolution spectroscopy of exoplanets was pushed forward mainly by the detection of molecular species via the cross-correlation technique. This was proven to be very efficient in disentangling the stellar and planetary fluxes from the noisy observations that are additionally contaminated by telluric absorption. Recently, Brogi & Line (2019) made another step forward and suggested the attractive approach of using cross-correlation technique for the retrieval of atmospheric properties. This is done by the mapping of the cross-correlation coefficients to log-likelihood. After this mapping was performed, the standard Bayesian methods could be used to derive the best fit parameters and their errors. The choice of the function that maps cross-correlation coefficients to the log-likelihood space is arbitrary, but should satisfy some basic criteria to make retrievals unambiguous (see Brogi & Line 2019, for more details). The new mapping proposed by Brogi & Line (2019) has certain benefits. For instance, it can be applied to already existing data when the planetary lines are drawn in the noise-dominated spectra.

Mapping of the cross-correlation coefficient and the OE method both rely on the presence of hundreds or even thousands of atmospheric lines, although the former method works with the residual spectra, while our implementation of the OE method works with spectra in arbitrary units and the corresponding flux scaling is a free parameter(s) of the model that could be optimized simultaneously with other (atmospheric) parameters. In our particular case the application of the new cross-correlation mapping was not strongly justified because the number of free parameters in our model is too high, making Bayesian post-processing very time consuming. In addition, we note that with the sufficiently high S/N values of the planetary spectra that are expected to be obtained with future instruments, the profiles of individual planetary lines could be satisfactorily resolved and directly analyzed. This will likely make the use of any cross-correlation mapping unnecessary favoring other approaches. We note that cross-correlation itself will still be used as a powerful tool for disentangling stellar and planetary spectra.

In this work we analyzed five different spectral regions between 1 and 5 μm. In these regions the telluric absorption is expected to be relatively weak compared to the rest of the spectrum. In addition, these regions contain strong bands of molecular species that we investigated (H₂ O, CO, CO₂, CH₄). Another region where it is possible to study planetary atmospheres that we did not consider in this study is the one around 3.5 μm, which is dominated by numerous lines of CH₄ and H₂ O. We carried out several additional retrievals using the 3.50–3.75 μm spectral region, but retrieved much less information compared to the other four regions. We conclude that although this region contains numerous molecular lines that can be easily detected from the usual cross-correlation analysis as shown in de Kok et al. (2014), among others, all these lines are strongly blended with each other and are not favorable for retrieving accurate molecular number densities.

In this work we did not account for the instrumental effects such as flat-fielding, order merging, or wavelength calibration. Some of these effects may influence the shape of the observed spectrum (which is a common problem for every echelle spectrograph) and therefore have direct impact on the retrieved parameters. Unfortunately, at this moment we cannot quantify these effects for CRIRES+ because many parameters of the instrument are still under investigation and the actual values will only be known once the instrument has been tested against real observations. Obtaining the best quality planetary spectrum requires a very accurate removal of instrument effects corresponding to the ratio between stellar and planetary continuum fluxes, which is about ≈ 10⁻³ for known HJs. Obviously, this accuracy will be very hard to achieve. On the other hand, the global shape of individual molecular bands will likely be preserved. In this case only arbitrary flux normalization to each spectral region under investigation will be needed (as we did in this study) in order to carry out accurate retrievals. Such a normalization ideally keeps the ratio of the planetary (pseudo)continuum-to-line depth unchanged even if the information regarding the true stellar continuum was lost during the data processing stage. For instance, in the wavelength domain covered by CRIRES+, the stellar continuum is very smooth. Therefore, instead of using simple continuum scaling, we also attempted to fit low-order polynomials to the simulated observations in each spectral region and did the same to the predicted spectra at each iteration in our retrievals. By doing this exercise we found that the choice of normalization function does not affect the final result very much. The temperature and mixing ratios varied from retrieval to retrieval, but always within the error bars. This exercise is very simplistic and cannot reflect all the details of the data reduction process planned for CRIRES+ and should not be considered as an approximation for reality. Nevertheless, it tells us that information about the true stellar continuum is not likely to be a major obstacle toward atmospheric retrievals as long as the relative shape of spectral lines is preserved. We note that in our approach the pseudo-continuum could be fit using low-order polynomials to be found during the retrieval process. This leaves the S/N of the planetary spectrum as the main limitation for atmosphericretrievals. Thus, accurate retrievals can still be performed provided that the order merging within a single CRIRES+ setting is accurately done by the pipeline of the instrument. We plan to improve our simulations by considering realistic instrumental effects in our next work.

Even if all the instrumental effects are properly taken into account, the telluric absorption represents the next obvious problem. Telluric lines are present in all near-infrared bands covered by CRIRES+. Even in the five carefully chosen regions that we investigated in this paper, many of the telluric lines can be strong depending on local weather conditions. Normally, moderately strong telluric lines can be successfully removed (e.g., Yan & Henning 2018; Brogi et al. 2016), while spectral regions containing the strongest telluric lines are excluded from the analysis. However, in cases when telluric lines need to be removed, the accuracy of this procedure should again be very high, about ≈ 10⁻³ or better. More importantly, there might be numerous very weak telluric lines, with intensities of just a fraction of a percent relative to the stellar continuum, that are unaccounted for in existing telluric absorption models and that could potentially distort the shape of planetary lines. We note that this is of minor importance if only a detection of exoplanetary lines using cross-correlation technique is needed. In that case, the amplitude of the cross-correlation function depends on matching the position of exoplanet lines against a chosen template spectrum. Thus, even if the depth of exoplanetary lines is affected by unremoved telluric absorption, it will not hide the detection of the molecular species. However, it could still be an issue when making retrievals from those distorted planetary lines. The overall impact of these very weak telluric lines should be tested with real observations. One obvious advantage of CRIRES+ is the ability to observe many molecular lines simultaneously, which will surely help to reduce the impact of telluric absorption in atmospheric retrievals.

In our simulated observations we used mixing ratios of four molecules (H₂ O, CO, CO₂, CH₄) as derived in Lee et al. (2012). Our retrievals showed that we can derive number densities of H₂ O, CO, and CO₂ by using different spectral regions. At the same time, we could not constrain a mixing ratio of CH₄ because the assumed mixing ratio is very low (1.9 × 10⁻⁷) and the lines of CH₄ are too weak compared to the lines of other three molecules. This is in agreement with the predictions of equilibrium chemical calculations (e.g., Gandhi & Madhusudhan 2017; Malik et al. 2017; Mollière et al. 2015). We note that CH₄ has a very rich spectrum in the near-infrared with many bands present in the wavelength range of CRIRES+. In atmospheres of HJs with equilibrium temperatures below 1000 K or with C/O > 1, the lines of CH₄ should become strong enough for the purpose of abundance analysis. As a test case we simulated an additional set of observations with the increased methane mixing ratio of CH₄ = 5 × 10⁻⁴ and performed retrievals from the (3.80–4.00) μm region and from the combined (2.28–2.38)+(3.80–4.00) μm region. In both cases we were able to obtain definite constraints on the methane mixing ratio, with slightly more accurate values when derived from the combined (2.28–2.38)+(3.80–4.00) μm region, just as expected. Thus, using CRIRES+ it will also be possible to constrain methane concentration in planets with cooler atmospheric temperatures than those considered here.

7 Summary

High-resolution spectroscopy in the near-infrared is a powerful tool for the study of absorption and emission lines to estimate the state of extra-solar atmospheres. Until now, it was mainly used for the detection of different molecular species via transit spectroscopy. In this work we made another step forward and investigated the potential of high-resolution spectroscopy to study molecular number densities and temperature stratification in exoplanet atmospheres by applying the optimal estimation method to the simulated observations at different near-infrared bands. We used a modified τ-REx forward model to simulate emission (i.e., out-of-transit) spectra for different molecular mixing ratios, spectral resolutions, and S/N values. The results of our investigation are summarized below.

• The spectra of HJs in the wavelength range of CRIRES+ is sensitive to a very wide range of atmospheric pressures between 1 and 10⁻⁶ bar.

• Avoiding regions with strong telluric contamination, the best region for deriving mixing ratios of H₂O, CO and CO₂ is around 1.6 μm, and the best region for studying temperature stratification and mixing ratios of H₂O and CO is around 2.3 μm. The mixing ratio of CH₄ is impossible to constrain accurately in any wavelength range and at the S/N that we considered, due to the very low number density that we assumed in our synthetic observations. It can be retrieved though, for example in cooler exoplanets that may have a CH₄ mixing ratio of ≈10⁻⁴–10⁻³.

• Retrieving mixing ratios simultaneously from two separate spectral regions helps to obtain accurate results for H₂O, CO, and CO₂. In this regard the combination of the 1.6 and 2.3 μm regions looks promising at any S/N, and the combination of 2.3 μm with any other region except 3.9μm when S∕N ≥ 10. However, the retrieval of accurate number densities always requires higher S/N compared to the S/N required for the temperature retrievals. The latter can be accurately retrieved even with the lowest S∕N = 5 that we assumed in our simulations.

• In the case of single setting retrievals, observing at the 2.3 μm region provides accurate VMRs of H₂O and CO only if S∕N ≥ 10. When S∕N ≥ 25 the 1.6 μm region provides very accurate VMRs of H₂O, CO, and CO₂. However, to achieve such high S/N values would exceed typical allocated telescope observational times. The best strategy for obtaining reliable mixing ratios of as many molecules as possible is thus to observe at two different infrared bands, at 2.3 μm and at any of the others except the one at 3.9 μm.

• Degrading spectral resolution results in decreasing the sensitivity of our retrievals of molecular concentrations, but the temperature distribution can still be accurately derived. Rebinning the spectra to a lower resolution could help to increase S/N and to obtain estimates on the molecular number densities and atmospheric temperature structure in cases when the S/N of the original data is not sufficient (e.g., S/N < 5).

• With CRIRES+ it will bepossible to carry out detailed atmospheric retrievals in two HJs known so far. Planet-hunting missions such as TESS⁸ , PLATO⁹ , and CHEOPS¹⁰ are expected to significantly increase the number of exoplanets accessible for current and future ground-based spectroscopic studies.

Our study suggests that even though the spectra of exoplanets will likely be obtained with very high noise levels, the atmospheric retrievals can still benefit from numerous molecular lines observed thanks to the wide wavelength coverage of CRIRES+. Moreover, additional important constraints can be provided by utilizing low-resolution data (Brogi et al. 2017). This data is already available for a number of HJs and could be used to break the degeneracy between the retrieved parameters. Thus, in our future work we plan to make another step forward and create a suite for the simultaneous retrieval of atmospheric parameters from low- and high-resolution data.

Acknowledgement

We wish to thank L. Fletcher for providing us with the data of HD 189733 b. We acknowledge the support of the DFG priority program SPP-1992 “Exploring the Diversity of ExtrasolarPlanets” (DFG PR 3624 602/41 and RE 1664/16-1). This work made use of the GWDG IT infrastructure. The GWDG is a shared corporate facility of the Georg-August-University Göttingen and the nonprofit Max-Planck Society (Max-Planck Gesellschaft zur Förderung der Wissenschaften e. V.). We also acknowledge the use of electronic databases SIMBAD and NASA’s ADS.

Appendix A Retrievals of temperature and mixing ratios from high-resolution spectroscopy

Fig. A.1 Retrieved temperature and mixing ratios from five spectral regions at spectral resolution R = 100 000 and S/N = 5. In each panel the first plot is the temperature distribution as a function of atmospheric pressure (solid blue line) and with error bars shown as shaded area. The red crosses and green dashed line are the true solution and initial guess, respectively. The second plot shows averaging kernels derived for the temperature distribution. Here averaging kernels for different atmospheric depths are randomly color-coded for better representation. The next four plots are the values of the retrieved mixing ratios of four molecular species (color-coded as in the first plot). The mixing ratios were assumed to be constant with atmospheric depth, and we show their values on the y-axis. The last plot compares the best fit predicted spectrum (solid red line) with the observed one (black circles).

Open with DEXTER

	Fig. A.2 Same as Fig. A.1, but for S/N = 25.
Open with DEXTER

	Fig. A.3 Same as Fig. A.1, but for S/N = 50.
Open with DEXTER

	Fig. A.4 Retrieved temperature and mixing ratios of four molecular species from a combination of different spectral regions with a reference region 2.28–2.38 μm. The retrievals are shown for the case of S/N = 5 and spectral resolution R = 100 000. The color-coding is the same as in Fig. 3.
Open with DEXTER

	Fig. A.5 Same as Fig. A.4, but for S/N = 25.
Open with DEXTER

	Fig. A.6 Same as Fig. A.4, but for S/N = 50.
Open with DEXTER

References

All Tables

All Figures

	Fig. 1 Top panel: observed day side fluxes of HD 189733b and our best fit model predictions (see legend inthe plot). Bottom panels: retrieved T-P profile along with corresponding averaging kernels and mixing ratios of molecular species. We used the T-P profile of Lee et al. (2012) a priori (shown as red crosses). Our best fit model is shown with a solid blue line and the shaded area represents 1σ error bars. The assumed a priori are shown as a green dashed line (which coincides with the red crosses in this particular case). The same color-coding is used for the four side plots of mixing ratios (on the right).
Open with DEXTER
In the text

	Fig. 2 Left panel: retrieved temperature distribution assuming four different initial guesses: three homogeneous at 1000, 1500, and 2000 K (vertical dashed lines), and the one from Waldmann et al. (2015a) (blue dashed line). The retrieved best fit profiles are shown with solid lines along with corresponding error bars as shaded areas. The profile from Lee et al. (2012) is shown as a solid red line. Right panels: averaging kernels for the four retrieved temperature distributions (from left toright and from top to bottom, color-coded accordingly).
Open with DEXTER
In the text

	Fig. 3 Retrieved temperature and mixing ratios from five spectral regions with R = 100 000 and S/N = 10. In each panel, the first plot compares the best fit predicted spectrum (solid red line) with the simulated observations (black). The second plot is the temperature distribution as a function of atmospheric pressure (solid blue line) and with error bars shown as shaded areas. The red crosses and green dashed line are the true solution and initial guess, respectively. The third plot shows averaging kernels derived for the temperature distribution. Here averaging kernels for different atmospheric depths are randomly color-coded for clarity. The next four plots are the values of the retrieved mixing ratios of four molecular species (color-coded as in the first plot). The mixing ratios were assumed to be constant with atmospheric depth, and we show their values on the y-axis.
Open with DEXTER
In the text

	Fig. 4 Retrievals of the H2O and CO mixing ratios from the (2.28–2.38) μm region assuming three different initial guesses (from top to bottom) and four S/N values (from left to right). The green dashed and red dotted lines are the initial guess and true solutions, respectively. The solid blue line and the shaded area are the retrieved value and its 1σ uncertainty, respectively.
Open with DEXTER
In the text

	Fig. 5 Comparison of retrieved temperature and mixing ratios of four molecules assuming S/N = 10 and two spectral resolutions of R = 100 000 and R = 50 000 for the top and bottom panels, respectively. The color-coding is the same as in Fig. 3.
Open with DEXTER
In the text

	Fig. 6 Retrieved temperature and mixing ratios of four molecular species from a combination of different spectral regions with a reference region 2.28–2.38 μm. The retrievals are shown for the case of S/N = 10. The color-coding is the same as in Fig. 3.
Open with DEXTER
In the text

	Fig. A.1 Retrieved temperature and mixing ratios from five spectral regions at spectral resolution R = 100 000 and S/N = 5. In each panel the first plot is the temperature distribution as a function of atmospheric pressure (solid blue line) and with error bars shown as shaded area. The red crosses and green dashed line are the true solution and initial guess, respectively. The second plot shows averaging kernels derived for the temperature distribution. Here averaging kernels for different atmospheric depths are randomly color-coded for better representation. The next four plots are the values of the retrieved mixing ratios of four molecular species (color-coded as in the first plot). The mixing ratios were assumed to be constant with atmospheric depth, and we show their values on the y-axis. The last plot compares the best fit predicted spectrum (solid red line) with the observed one (black circles).
Open with DEXTER
In the text

	Fig. A.2 Same as Fig. A.1, but for S/N = 25.
Open with DEXTER
In the text

	Fig. A.3 Same as Fig. A.1, but for S/N = 50.
Open with DEXTER
In the text

	Fig. A.4 Retrieved temperature and mixing ratios of four molecular species from a combination of different spectral regions with a reference region 2.28–2.38 μm. The retrievals are shown for the case of S/N = 5 and spectral resolution R = 100 000. The color-coding is the same as in Fig. 3.
Open with DEXTER
In the text

	Fig. A.5 Same as Fig. A.4, but for S/N = 25.
Open with DEXTER
In the text

	Fig. A.6 Same as Fig. A.4, but for S/N = 50.
Open with DEXTER
In the text