A Near-infrared Chemical Inventory of the Atmosphere of 55 Cancri e

After the initial reduction process described above, the resultant spectra were largely dominated by both stellar and telluric absorption lines, seen as vertical lines in, for example, the second left and second right panels of Figure 1. However, due to the essentially stationary nature of these lines when compared with absorption due to the exoplanet's atmosphere (which varies in radial velocity from ∼−60 km s⁻¹ to ∼+60 km s⁻¹ over the course of the transit), these features can be removed using the sysrem algorithm, first described by Tamuz et al. (2005). sysrem allows for linear, systematic effects that appear in many data sets of the same sample to be removed through an algorithm that, in the case of equal uncertainties, reduces to a principal component analysis (PCA) algorithm.

The correction was determined using both in- and out-of-transit frames, and each order of the data was treated separately. For all observations, the average air mass at each exposure was taken as an initial approximation of the systematic effects to be removed. For consistency, we chose to apply the same number of iterations of the algorithm to each order; to determine this number, we calculated the rms of the residuals after subsequent iterations of sysrem. We found that on average, six iterations of the algorithm was the point at which the rms of the residuals had plateaued and additional applications of the algorithm did not yield a significant improvement. The third left and third right panels of Figure 1 show the results of applying six iterations of the sysrem algorithm to CARMENES and SPIRou spectra, respectively. The results of applying six iterations to all nights and orders of each data set are shown in the figures in Appendix A.

The algorithm performs poorly around strong or closely spaced lines; in particular, several orders of our observations (i.e., those at the edges of the Y, J, H, and K photometric bands) suffer from significant telluric absorption. These strong lines make the blaze response correction difficult, and as a result, nonlinear effects are introduced into the spectra that cannot be removed using sysrem. To reduce contamination from these poorly constrained frames, we followed, for example, Snellen et al. (2010), Esteves et al. (2017), and Deibert et al. (2019) in weighting each frame and each pixel by its standard deviation (shown in the fourth panels of Figure 1, as well as the figures presented in Appendix A). This step reduces the contribution from the noisier portions of the data. Additionally, we have chosen to fully exclude several particularly contaminated orders (those located within the gaps between the Y, J, H, and K bands) from further analysis. These are described in greater detail in Appendix A.

3.1.1. Telluric Feature Removal with Molecfit

Our analysis involves independently testing for the presence of various molecules (HCN, NH₃, C₂H₂, CO, CO₂, and H₂O) in the atmosphere of 55 Cnc e through cross-correlation with so-called "parametric models" similar to those used in Jindal et al. (2020). We explore a range of VMRs and mean molecular weights (μ's) of each of these molecules, described in further detail below. We note that these models are not self-consistent and are instead meant to provide us with an initial insight into the exoplanet's atmosphere and the limitations of our method. Furthermore, various models (i.e., those with low mean molecular weights; see the top left panels of the figures in the following section) are nonphysical. Such models allow us to test our method, as they are easily recoverable by the model injection/recovery process (see Section 4.3). We describe the calculations of all models below.

4.1.1. Atmospheric Model Calculation

For each night of observations, our data were cross-correlated with models at Doppler shifts spanning −250 km s⁻¹ to 250 km s⁻¹, with 1 km s⁻¹ steps between each velocity. This is shown in the top left panels of Figure 2. Next, we phase-folded the correlation signal from the in-transit frames by shifting each correlation to the reference frame of 55 Cnc e. The in-transit frames were determined using a model light curve generated with the occultquad package from Mandel & Agol (2002). The parameters used for this model light curve, including limb-darkening parameters from Claret (2004), are summarized in Table 2. In order to account for the fact that the planetary radial velocity is not known with perfect precision, as well as to better understand the behavior of the cross-correlation function at velocities outside of the planetary rest frame (e.g., to investigate spurious correlation peaks due to residual telluric lines), we created a correlation map over a wide range of systemic velocities (V_sys) and planetary orbital velocities (${K}_{p}\sin 2\pi \phi (t)$, where ϕ(t) is the orbital phase). We created the map for planetary RV semiamplitude (K_p ) values ranging from 1 to 300 km s⁻¹, with a 0.5 km s⁻¹ step between each value. An example is shown in the bottom left panels of Figure 2.

For each model, these correlation grids were summed over orders containing significant molecular absorption, and then summed over all observations.

To assess the robustness of our pipeline and to place constraints on the atmospheric composition of 55 Cnc e, we performed injection and recovery tests for each model used in our analysis. The goal of these tests was to determine whether or not our data-reduction process, including the sysrem algorithm, removes any signal that may be present, as well as to place sensitivity limits on our results.

We carried out these injection and recovery tests by multiplying the in-transit frames of each observation (prior to any data reduction past initial processing at the telescope) by an atmospheric model shifted to the frame of the exoplanet to create a synthetic "model + data" data set. Note that the model was only added to the in-transit frames of our data. This was done for each model described in Section 4.1, each order of every observation, and each night of observation. The synthetic "model + data" spectra were then processed through our data-reduction pipeline described in Section 3 and analyzed using the Doppler cross-correlation method described in Section 4.2. An example of the process is shown in the middle two panels of Figure 2.

As a result, we were able to assess the sensitivity limits of our analysis, by determining which synthetic spectra were recovered by our pipeline. For example, the top left panel in Figure 3 shows a case where a correlation between a model atmosphere and the corresponding synthetic data set resulted in a detection (magenta line), whereas a correlation between that same model atmosphere and the true data set (black line) did not result in a detection. Likewise, the top right panel in Figure 3 shows a case where neither a correlation between a model atmosphere and the corresponding synthetic data set (magenta line) nor the correlation between that same model atmosphere and the true data set (black line) resulted in a detection. The model atmosphere in question (HCN with a VMR of 0.1% and a mean molecular weight of 10 amu) was beyond the sensitivity limits of our detection pipeline.

We note that various models with nonphysical, low mean molecular weights (e.g., the top left panel in Figure 8) allow us to confirm that the model injection/recovery process is working as expected. Such models are easily recovered by the injection/recovery process, as seen in the top left panels of the figures in the following section.

This model injection/recovery process thus allowed us to place constraints on both the atmospheric makeup of the exoplanet and the detection limits of our observations.

In order to assess the significance of our results, we followed the methods described in Esteves et al. (2017) and Deibert et al. (2019). We derived 1σ and 3σ confidence levels for each feature in each night of the data by randomly selecting a set of out-of-transit frames corresponding to the number of in-transit frames for each night, assigning an in-transit phase to each of these spectra, and following through with the cross-correlation and phase-folding process as described in Section 4.2. We repeated this process 10,000 times, sorted the data, and then selected the 1σ and 3σ confidence levels based on the outcome. These levels were then compared to the correlation strengths of each model in order to assess their significance. An example is shown in the rightmost panel of Figure 2.

Figure 3 shows the results of our search for HCN in the atmosphere of 55 Cnc e. Our aim was to further investigate the results of Tsiaras et al. (2016), who used HST/WFC3 observations to report the detection of an atmosphere and suggested that HCN is the most likely molecule able to account for the observed absorption features.

We analyzed a range of atmospheric models varying in VMR from 0.1% down to 5 × 10⁻⁵% and in mean molecular weight from 2 to 20 amu. Our model spectra were generated using HITRAN2016 (Gordon et al. 2017). We note that the results depend on the accuracy of the line list used, so future studies with updated line lists may yield different outcomes.

As seen in Figure 3, we can rule out the presence of HCN in the atmosphere of 55 Cnc e at a mean molecular weight of 2 amu with a VMR of 0.001%; if the mean molecular weight is increased to 5 amu, the lowest VMR that we can rule out is 0.02%.

Zoom In Zoom Out Reset image size

Figure 4 shows the results of our search for NH₃ in the atmosphere of 55 Cnc e.

Zoom In Zoom Out Reset image size

We analyzed a set of models ranging in VMR from 5% to 10⁻⁸% and mean molecular weight ranging from 2 amu to 20 amu.

As can be seen in Figure 4, we are able to rule NH₃ out of the atmosphere at a mean molecular weight of 2 amu with a VMR as low as 0.0025%; if the mean molecular weight is increased to 5 amu, we can still rule NH₃ out with a VMR as low as 0.08%.

Figure 5 shows the results of our search for C₂H₂ in the atmosphere of 55 Cnc e.

Zoom In Zoom Out Reset image size

We analyzed a set of models ranging in VMR from 20% to 10⁻⁸% and mean molecular weight ranging from 2 amu to 20 amu.

As can be seen in Figure 5, we are able to rule C₂H₂ out of the atmosphere of 55 Cnc e at a mean molecular weight of 2 amu with a VMR as low as 0.08%, and at a mean molecular weight of 5 amu with a VMR as low as 1.0%.

The results of our search for CO in the atmosphere of 55 Cnc e are summarized in Figure 6.

Zoom In Zoom Out Reset image size

We analyzed a range of atmospheric models varying in VMR from 10% down to 10⁻⁶% and in mean molecular weight from 2 to 30 amu.

As shown in Figure 6, we are unable to detect atmospheric CO in our data. However, we are able to tentatively rule out the possibility of CO being present in the atmosphere of 55 Cnc e at a mean molecular weight of 2 amu and a VMR as low as 1.0% at the 3σ level. Note, however, that there are numerous additional features with peaks >1σ that are likely due to noise in the data (see, e.g., the discussion in Esteves et al. 2017, who noticed a similar phenomenon for optical transit data of 55 Cnc e). For this reason, we caution that our limits on the presence of CO are only tentative and warrant further investigation.

The results of our search for CO₂ in the atmosphere of 55 Cnc e are summarized in Figure 7.

Zoom In Zoom Out Reset image size

As was the case with CO (see Section 5.4), we analyzed a range of atmospheric models varying in VMR from 10% down to 10⁻⁶% and in mean molecular weight from 2 to 30 amu.

As can be seen in Figure 7, we are neither able to detect atmospheric CO₂ in our data nor place any significant constraints on its presence. We note that while a peak is present in the injected data at a systemic velocity of 0 km s⁻¹ and the orbital velocity of 55 Cnc e, as would be expected for a detection, it does not surpass 3σ and therefore is not sufficiently significant to make any conclusions about the presence or lack of atmospheric CO₂.

The results of our search for H₂O in the atmosphere of 55 Cnc e are summarized in Figure 8.

Zoom In Zoom Out Reset image size

We analyzed a range of atmospheric models varying in VMR from 20% down to 10⁻⁸% and in mean molecular weight from 2 to 20 amu.

As was the case with CO₂ (see Section 5.5), none of the injected models seen in Figure 8 surpass the 3σ confidence level, and thus we are unable to place any limits on the presence of H₂O in the atmosphere of 55 Cnc e with these data alone. This is likely due to the imperfect removal of telluric water lines in the data (see Section 3.1 and the additional discussion in Section 6.2 and Appendix D for further details).

Our model injection/recovery tests for HCN (see Figure 3) significantly narrow the parameter space of potential atmospheres that are consistent with the results reported in Tsiaras et al. (2016). These authors sampled VMRs from 1 to 10⁻⁸ and mean molecular weights from 2 to 10 amu. Their analysis indicates that the mean molecular weight of the atmosphere peaks at roughly 4 amu, with higher values being unlikely. Furthermore, strong absorption seen near 1.4 and 1.6 μm in their analysis indicates that the mean molecular weight is relatively low. They also reported that (of the molecules considered in their fit) HCN is the most likely absorber able to explain their observed absorption features, and that a scenario with a high VMR for HCN is favored, though values as low as 0.001% are acceptable.

While our results do rule out the most likely scenarios considered by Tsiaras et al. (2016), they are not in complete disagreement. A mean molecular weight of 2 amu at VMRs consistent with the Tsiaras et al. (2016) analysis is ruled out by our observations; however, it is possible that the mean molecular weight is slightly greater (∼4 or 5 amu) and that the VMR is between 0.02% and 0.001%. If the detection in Tsiaras et al. (2016) is real, a higher mean molecular weight would be unlikely, since that would significantly mute the features in the transmission spectrum.

While additional observations will be necessary to determine whether a scenario with low VMR and high mean molecular weight can be ruled out completely, our analysis has provided significantly improved constraints on the presence of HCN and vastly decreased the parameter space of possible atmospheres. Our analysis has also provided constraints on the presence of NH₃ and C₂H₂ and suggests that if any of these three molecules are indeed present in the atmosphere of 55 Cnc e, they are likely to have high mean molecular weights, low VMRs, or both. We note that previous observations (e.g., Ehrenreich et al. 2012; Demory et al. 2016b) do support an atmosphere with high mean molecular weight.

Finally, we caution that these constraints rely on the accuracy of the line lists used to generate our models, and future updates to these line lists could yield different results. The issue has been discussed in detail in a number of previous studies. In the case of TiO, for example, Hoeijmakers et al. (2015) demonstrated that inaccurate line lists can hamper high-resolution retrievals. Likewise, Webb et al. (2020) observed absorption due to H₂O in high-resolution VLT/CRIRES spectra of HD 179949 b, but found a weak dependence on the line list used for the cross-correlation. Updates to the available line lists for each of the molecules used in this analysis could similarly impact our detection capabilities, and we therefore note that the particulars of line lists should be kept in mind when interpreting our results.

We note that previous studies (Esteves et al. 2017; Jindal et al. 2020) have been able to place significant constraints on the presence of H₂O in the atmosphere of 55 Cnc e using high-resolution spectra at optical wavelengths. In particular, recent work by Jindal et al. (2020) resulted in a 3σ lower limit of 15 amu on the mean molecular weight of 55 Cnc e's atmosphere, assuming it is water-rich (i.e., has a VMR > 0.1%) and that it does not have thick clouds. Our observations do not improve on these constraints.

Given the fact that the NIR wavelength ranges of CARMENES and SPIRou contain more absorption lines due to H₂O than the optical wavelength ranges of these previous works (e.g., compare our Figure 28 with Figure 8 of Esteves et al. 2017 and Figure 2 of Jindal et al. 2020), it could be expected that our data (which cover the same number of transits as Jindal et al. 2020) should give the same or even better constraints than those analyses. To investigate the discrepancy, we carried out a number of tests designed to probe the limits of our technique in the context of NIR H₂O (and CO₂) absorption. The full details of these tests are presented in Appendix D.

We conclude that the discrepancies between this work and Esteves et al. (2017) and Jindal et al. (2020) are largely due to our inability to fully remove telluric water features from our data. While the overall wavelength coverage of the observations used in this work is indeed broader than that of Jindal et al. (2020), and therefore contains a greater number of absorption lines due to H₂O, the telluric contamination at these redder wavelengths is far more severe than at the wavelengths considered in Jindal et al. (2020). As the individual absorption lines due to the exoplanet's atmosphere are weak, our method relies heavily on our ability to combine the signals from many absorption lines, which in turn relies on our ability to remove telluric contamination adequately in the regions where these lines are present. sysrem is less efficient in regions of strong telluric contamination; this can be seen in a comparison between the residuals after applying sysrem in Jindal et al. (2020) and those in this work (see the plots in Appendix A). While additional iterations of sysrem may aid in removing some residual telluric contamination, we show in Appendix D that additional iterations may also begin to remove the model itself and thus will not improve our detection capabilities. We also note that the signal-to-noise raitos (S/Ns) of our observations (see Table 1) are in some cases significantly lower than those of Jindal et al. (2020), which would limit our sensitivities further.

As mentioned in Section 3.1.1, we also tried correcting for telluric absorption using Molecfit. We applied Molectit and Calctrans to every spectrum in each night of our CARMENES and SPIRou data. We then injected an atmospheric H₂O model with a VMR of 20% and a mean molecular weight of 2 amu (i.e., the model shown in Figure 28, the strongest water model of our sample) into the data and repeated this process. Then, we carried out the Doppler cross-correlation process as described in Section 4.2 in order to determine whether or not we were more readily able to recover an injected atmospheric H₂O model with sysrem or with Molecfit.

We found that Molecfit does not improve our ability to detect an injected atmospheric H₂O model. Instead, the data corrected with Molecfit yielded a weaker detection of the injected model than those corrected with sysrem. This is evident in Figures 9 and 10, where we compare the two for CARMENES and SPIRou observations, respectively. We tested two separate methods for removing stellar lines after applying Molecfit to our data: first, we subtracted a mean template of the stellar features remaining; and second, we applied two iterations of sysrem to the data after correcting them with Molecfit. In the case of CARMENES observations, we are unable to recover the injected model (the middle two panels of Figure 9); for SPIRou, on the other hand, we can recover the injected model at a lower significance than our original analysis (the second panel of Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a check for our CARMENES observations, we also tested a case where we injected the model at 100 times its nominal strength into a subset of our data and followed the same reduction process using Molecfit. In this case, we were able to recover the injected model at the expected location.

We note that previous studies that have compared Molecfit and sysrem in the NIR have had similar results. In particular, Sánchez-López et al. (2019) made a detection of water vapor in the atmosphere of HD 209458 b using CARMENES data that had been corrected with sysrem. After correcting their data with Molecfit, however, they were unable to recover the planetary signal. The authors noted that recent works have found the scatter in residuals for CRIRES data corrected with Molecfit to be 3%–7% (Ulmer-Moll et al. 2019) and concluded that the fit uncertainties over the wavelength coverage of CARMENES are likely too large to detect the H₂O features of the planet. Our analysis yields similar results: although we do not strongly detect the atmospheric H₂O model with sysrem, there is a correlation peak at the expected location that approaches 3σ. When running Molecfit on our CARMENES data, however, no such peak is present. In the case of SPIRou, a peak is visible at the expected location, but this signal is weaker than the case where sysrem was used instead (Figure 10).

As a final test of what might be impacting our ability to recover injected H₂O models, we followed Alonso-Floriano et al. (2019) and Sánchez-López et al. (2019) in separately analyzing individual bands of our observations. We separated the CARMENES data into the Y, J, and H bands and separated the SPIRou data into the Y, J, H, and K bands.

Figures 11 and 12 show these results for CARMENES and SPIRou, respectively. In both cases, the contribution to the signal comes almost exclusively from the J band (and, to a lesser extent, the H band). This suggests that although our observations span a wide wavelength range, only a small fraction of that wavelength range is contributing to our detection capabilities. The other water features present (see Figure 28) were likely not readily recovered by our model injection/recovery tests because they are comparatively weaker and because sysrem was not able to remove telluric features completely in these regions (see the figures in Appendix A, particularly the bottom two panels of each).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also ran the individual band analysis on our strongest HCN model as a comparison. The results are shown in Figures 13 and 14 for CARMENES and SPIRou, respectively. In this case, we see that a signal is recovered in all bands except the K band. For CARMENES the signal is strongest in the Y and H bands, whereas for SPIRou the signal is noticeably strongest in the J and H bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Several previous studies have pointed to a nitrogen-dominated atmosphere for 55 Cnc e. Hammond & Pierrehumbert (2017) used 3D calculations to show that the offset observed in the phase curve (Demory et al. 2016b; Angelo & Hu 2017), as well as the large day–night temperature contrast, may be explained by an N-dominated atmosphere. The analysis by Angelo & Hu (2017) also supports this scenario. To this end, Miguel (2019) explored the observable features in the spectra of an N-dominated atmosphere for 55 Cnc e. Through analytical arguments, equilibrium chemistry calculations, and adopting Titan's elemental abundances as a potential composition, Miguel (2019) showed that although N₂ is expected to be the most abundant molecule in the atmosphere (followed by H₂ and CO), the transmission spectra should show strong features of NH₃ and HCN. However, a decrease in the N/O ratio would tend to weaken these NH₃ and HCN features.

The transmission spectra calculated in Miguel (2019, their Figure 3, right panel) range between 3 and 20 μm and are thus outside the range of our own observations. Figures 6 and 7 of Miguel (2019) show the mixing fraction of the most abundant observable molecules as a function of temperature at a pressure of ∼1.4 bar (as calculated by Angelo & Hu 2017) and the mixing fractions as a function of pressure at different N/O ratios, respectively. For all temperatures considered by Miguel (2019) at a pressure of ∼1.4 bar, the VMR of NH₃ is less than ∼10⁻⁴%, while at all pressures considered and for all N/O ratios considered, the VMR of NH₃ is less than ∼10⁻⁵%. Our analysis rules out most low-mean-molecular-weight, high-VMR scenarios, which is consistent with expected VMRs from Miguel (2019) and with previous work that has combined mass and radius measurements with interior modeling to show that the atmosphere should have a high mean molecular weight (Demory et al. 2011; Winn et al. 2011; Bourrier et al. 2018).

We note that the limits we have placed throughout this section assume that any signals present in the atmosphere are not obscured by clouds or hazes. However, the presence of clouds or hazes could act to obscure atmospheric signals, limiting our ability to make detections. Therefore, our results may also be consistent with a cloudy atmosphere. The limits we have placed throughout this analysis are only relevant in the case of a cloud-free atmosphere, as the models were injected into the data under the assumption that they were not obscured at any altitude by clouds or hazes. We note as well that various other investigations into the composition of 55 Cancri e's atmosphere (Tsiaras et al. 2016; Esteves et al. 2017; Jindal et al. 2020) made similar assumptions, meaning that those limits also pertain to the case in which 55 Cnc e does not harbor clouds.

While an in-depth exploration of the effects of clouds or hazes on the atmosphere could prove insightful, such an analysis is beyond the scope of this work. However, we note that an analysis by Mahapatra et al. (2017) found that despite 55 Cnc e's high equilibrium temperature (∼2400 K, Demory et al. 2016a), conditions on 55 Cnc e may allow for mineral clouds to form. Such cloud formation is likely only possible in a thin atmospheric region and requires strong vertical replenishment (Mahapatra et al. 2017). Hammond & Pierrehumbert (2017) also investigated the possibility of clouds and found that given the observed high equilibrium temperature and a range of partial pressures from Miguel et al. (2011), Na is unlikely to condense and form clouds, whereas SiO could potentially condense on the planet's nightside (both Na and SiO could arise from a dayside magma ocean; e.g., Schaefer & Fegley 2010; Miguel et al. 2011; Hammond & Pierrehumbert 2017).

Briefly, we note that while the effects of refraction were not considered in our model atmosphere calculations (see Section 4.1), it has been shown that refraction can act to mute spectral features in the lower atmosphere by creating a gray continuum similar to that produced by optically thick clouds (Bétrémieux & Swain 2018). For thin atmospheres around terrestrial exoplanets, the largest pressure that can be probed may indeed be the exoplanet's surface pressure; with increasingly thicker atmospheres, however, refraction may prevent observations of the lower atmosphere.

Bétrémieux & Swain (2018) calculate the refractive boundaries of several hot Jupiters and terrestrial exoplanets for various atmospheric compositions in order to assess the impact of refraction on observations. In most cases, these refractive boundaries are located at pressures >1 bar.

In the case of 55 Cancri e, which has a high equilibrium temperature (∼2400 K, Demory et al. 2016a), and the observations and models presented in this paper, refraction should not have a large effect on our results. However, future work (particularly of observations that probe deeper into the atmosphere) may benefit from including the effects of refraction in model calculations.

Figures 15–18 show the results of applying the full data-reduction process to the four nights of CARMENES observations used in this analysis. We chose to exclude several orders from all further analyses due to severe telluric contamination that prevented the blaze correction or the sysrem algorithm from being able to provide an adequate reduction; these are seen in regions where the standard deviation is much higher than the rest of the spectrum. Note that we also chose to exclude N₄ and N₆ from our analysis because of poor observing conditions (see Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 19–22 show the results of applying the full data-reduction process to the four nights of SPIRou observations. Again, we have excluded regions of severe telluric contamination from our analysis and indicated these in gray scale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The initial parameters used with Molecfit for every night are shown in Table 3. Further details on these parameters are available in Smette et al. (2015) and Kausch et al. (2015). We allowed Molecfit to fit for a Gaussian kernel that varies with wavelength and included fitting for the wavelength solution.

Table 3. Initial Parameters Used by Molecfit

Initial Parameters	Value	Notes
ftol	10−10	Relative χ2 convergence criterion
xtol	10−10	Relative parameter convergence criterion
Molecules	H2O, CO2, O3, CO CH4, O2	Molecules included in the fit
contn	2	Degree of coefficients for continuum fit
cont${}_{\mathrm{const}.}$	1	Initial constant term for continuum fit
nλ	2	Polynomial degree of the refined wavelength solution
ωGaussian	3.5	Initial value for FWHM of Gaussian in pixels
Kernel size	15	Size of Gaussian kernel in FWHM

Note. We allow Molecfit to fit for molecules with absorption features present in the wavelength range of our data (see Molecules row) and apply a second-degree Chebyshev polynomial fit to the wavelength solution (n_λ ). The continuum (cont_n ) is fit with a second-degree polynomial, and as the data are normalized, we set an initial constant term of 1 for the continuum, which is then refined through the fit. The instrumental profile is initially assumed to be a Gaussian with an FWHM of 3.5 pixels; however, this too is refined through the fit. Molecfit uses the Levenberg–Marquardt technique to quickly solve the least-squares problem (Smette et al. 2015); the χ² convergence parameter (ftol) and the parameter convergence criterion (xtol) of the Levenberg–Marquardt technique are set to 10⁻¹⁰.

Download table as:  ASCIITypeset image

To test whether residual noise left in the data after the application of the sysrem algorithm was affecting our final result, we followed a method similar to that of Esteves et al. (2017) by simulating a pure white noise data set generated to match the rms of each wavelength channel after applying sysrem (see, e.g., the bottom panels of the figures in Appendix A). We tried injecting both a model H₂O atmosphere and a model CO₂ atmosphere into this white noise data set and repeated our analysis routine as described in Section 4 for each model. This allowed us to assess whether the residual noise level in the data after applying the sysrem algorithm was the limiting factor in being able to recover H₂O and CO₂ models. The results of these white noise tests are shown in Figures 29 (H₂O) and 30 (CO₂) for the CARMENES data set, and Figures 31 (H₂O) and 32 (CO₂) for the SPIRou data set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all cases, we see that in a data set composed of pure white noise, we are unable to recover even the strongest H₂O and CO₂ models. This suggests that sysrem did not adequately remove stellar and telluric contamination in the data, and that the data themselves are simply too noisy to allow us to place strong constraints on the presence of H₂O or CO₂ in the atmosphere of 55 Cnc e. We note that H₂O and CO₂ contain a significant number of lines throughout the broad wavelength coverages offered by CARMENES and SPIRou; while these lines can in theory improve the strength of our detections, there will also be corresponding H₂O and CO₂ lines present in this wavelength range in the Earth's atmosphere that contaminate the data and pose a challenge to the sysrem algorithm. Future analyses would benefit from exploring other avenues of telluric absorption line removal. While we have presented a brief investigation of the efficacy of Molecfit on these data (see Sections 3.1.1 and 6.2, as well as Appendix C), we note that additional insight into the strengths of various telluric removal methods on SPIRou observations specifically would be useful. In the case of CRIRES data, we note that Ulmer-Moll et al. (2019) showed synthetic transmission methods to be preferable to a standard star method.

In addition to the white noise test described in Appendix D.1, we investigated whether or not a shift was present in the wavelength solution of the data. Such a shift could potentially impact our ability to recover an injected model. To test for this, we cross-correlated each spectrum of each night with the first spectrum of the night, and we used a centroid algorithm to fit for the peak of the cross-correlation function. We then compared the difference in measured peaks for each cross-correlation. This was done on a per-order basis.

In the case of the CARMENES data, we find that for each night the majority of orders do not experience a large shift in the wavelength solution. For most orders, the shift in the fitted centroid of the cross-correlation function is on average 0.04 pixels. For each night, we found that two orders experienced a slightly larger shift (on the order of ∼1 pixel); however, these were orders that had previously been discarded from the data due to excessive telluric contamination. The offending orders occur in the wavelength range between 1344.14 and 1400.24 nm; as can be seen in the figures presented in Appendix A, this region of the data is heavily contaminated by tellurics.

In the case of the SPIRou data, the shift in the fitted centroid of the cross-correlation function for each order is on average 0.03 pixels. We do not find that any orders experience a significantly larger shift.

In addition to the test described above, we also selected several individual telluric lines in each data set and fit these with a Gaussian, determining the centroid of each line by the Gaussian's expected value μ. For these individual lines, we also did not find a significant shift in the wavelength solution.

Finally, we note that our data-reduction procedure (described in Section 3) involves interpolating each frame to a common wavelength grid.

Overall, we do not think that a shift in the wavelength solution impacted our final results, given that the shift measured by the methods described above is significantly less than one pixel and less than the widths of the absorption lines themselves.

In addition to the tests described in the preceding sections, we investigated whether applying different numbers of iterations of sysrem affected our ability to recover injected models. While the methods we used to determine the optimal number of iterations were model-free (see Section 3.1), it is possible that sysrem had begun to overfit the data and remove signal from the injected model (see, e.g., the discussion in Brogi 2014). For the purposes of this test, we again used the water model with a VMR of 20% and a mean molecular weight of 2 amu (i.e., the model shown in Figure 28) and the CO₂ model with a VMR of 10% and a mean molecular weight of 2 amu (i.e., the model shown in Figure 27), because these models should be the easiest to detect.

The results are shown in Figures 33 and 34. The noise level decreases with increasing numbers of iterations of sysrem and appears to plateau at six iterations, as determined previously (see Section 3.1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, Figure 33 also shows that the strength of the injected water model decreases slightly after five iterations of the algorithm: while the injected model after five iterations surpasses the 3σ level, this is not the case for the injected model after six iterations. This indicates that the sysrem algorithm has begun to remove signal contributed by the model itself. In the case of CO₂, however, we find that the model is more readily detected after six iterations.

Although this test demonstrates that varying the number of iterations of sysrem can have an impact on our ability to recover an injected model, we have chosen not to change the number of iterations used in our analysis for several reasons. First, the number of iterations initially chosen was based on an average optimal result across all orders and was independent of the model or molecule being considered. We have chosen to implement a model-free method for determining the optimal number of iterations of sysrem so that the specifics of the models being injected are less likely to impact our results.

We also note that the change in strength between five and six iterations as seen in Figure 33 is not highly significant: while the injected model after five iterations slightly surpasses the 3σ level, the injected model after six iterations is close to (but not quite at) the 3σ level. While it is possible that the number of iterations of sysrem used in our analysis affected some of the results shown in Figure 8, these effects are small compared to the difference in the strengths of injected models in our study and those of Jindal et al. (2020), and we thus do not consider this to be the dominating factor in our ability to recover an injected model. Instead, our ability to recover an injected model appears to be limited by the overall noise level of the data, as demonstrated by the white noise test described in Appendix D.1. While sysrem does not appear to have adequately removed contributions from the Earth's atmosphere, it has also begun to remove the signal contributed by some injected models. Again, we believe that future analyses of data in the wavelength ranges covered by CARMENES and SPIRou would benefit from exploring alternate methods of removing stellar and telluric lines, including a detailed comparison between the capabilities of sysrem and Molecfit.