Contrast and Temperature Dependence of Multi-Epoch High-Resolution Cross-Correlation Exoplanet Spectroscopy

To set rough limits on the equilibrium temperature and photometric contrast required to make a detection, we simulate the thermal emission spectra of planets with temperature structures generated self-consistently for equilibrium temperatures 600 K, 1000 K, and 1400 K. We vary photometric contrast by scaling the planet models from 0.7 $\rm R_{J}$ to 1.3 $\rm R_{J}$ in steps of 0.1 $\rm R_{J}$. Note that we do not re-compute the planet model for each radius/photometric contrast to account for changes in gravity or temperature–pressure (T–P) profile. We simulate 10 epochs with $\rm S/N_{epoch}$ of 2500 and 25 and 50 epochs with $\rm S/N_{epoch}=1500$, for total S/N of 7900, 7500, and 10600, respectively. The variation in radius corresponds to an $L$-band photometric contrast range of $2.2\times 10^{-5}$ to $7.6\times 10^{-5}$ for $\rm T_{eq}=600$ K, $1.6\times 10^{-4}$ to $5.4\times 10^{-4}$ for $\rm T_{eq}=1000$ K, and $3.7\times 10^{-4}$ to $1.3\times 10^{-3}$ for $\rm T_{eq}=1400$ K. This allows the effect of variations in temperature given radius, variations in radius given temperature, and simultaneous variations in both radius and temperature to be distinguished and compared. A priori, we expect larger and hotter planets should be more easily detected for fixed stellar properties due to the smaller photometric contrast with the host star giving a better signal-to-noise on planet lines.

Previous applications of the 2D multi-epoch cross-correlation technique have found significant non-planetary structure in the observed $\log L(\rm K_{p})$ curves, well in excess of the errors estimated from jackknife tests (Piskorz et al., 2016, 2017, 2018; Buzard et al., 2020). These features are believed to arise from a combination of telluric residuals, features in the observed star/planet not represented in the template spectra, and the correlation between the planet template and observed stellar spectrum. In some cases (e.g. Lockwood et al., 2014; Buzard et al., 2020), these effects can be comparable in strength to the cross-correlation features caused by the planet, requiring careful identification and modeling to correctly identify the planet signal and complicating estimates of the detection strength. Proper echelle and cross-disperser settings can also help to minimize these effects by avoiding spectral regions prone to strong telluric residuals or windows where the stellar and planetary spectra are strongly correlated.

As our simulation framework does not include tellurics and uses the same model framework for both cross-correlation templates and simulating observed spectra, the off-peak structure (visible in Figure 2) should be dominated by the correlation between the planet template and the stellar spectrum. In order to remove structured off-peak correlation, a set of star-only observations containing no planet spectrum was simulated and the resulting $\log\bar{L}(\rm K_{p})$ curve subtracted from each of the simulated planet $\log L(\rm K_{p})$ curves to give a relative likelihood, $\log RL(\rm K_{p})$ (see equation 4). This has the effect of normalizing the likelihood surface by the correlation of the planet model with the observed stellar spectrum. The lack of significant structure in Figure 3 between -150 and 0 km $\rm s^{-1}$ compared with Figure 2 shows that the subtraction of the star-only simulation almost entirely eliminates the structured off-peak correlation in our simulations.

Simulations of the stellar spectrum were compared to observations of the combined star/planet spectrum in Buzard et al. (2020) and successfully reproduced some, but not all, observed non-planetary features in the log-likelihood surface, indicating that the near-perfect correction of structured correlation is not yet achievable in practice. A discussion of how this correction may be improved in observations is presented in Section 4. The inability to achieve good corrections for off-peak correlation in observations complicates the detection of smaller/cooler planets compared with simulations, and the results presented in Figure 3 and Table 1 likely overestimate the detectability of such planets with current pipelines. However, correcting the off-peak correlation enables a more robust comparison between different simulations, and offers insight on the possibilities with future pipeline improvements.

The $\log RL(\rm K_{p})$ curves for each $\rm T_{eq}$ and $\rm R_{p}$ are plotted in Figure 3, and detection likelihood ratios are listed in Table 1. Jackknife tests are used to estimate the shaded 1$\sigma$ error region. In addition to the total planet spectrum, we perform the cross-correlation with molecule-specific planet templates for H₂O, CH₄, CO, NH₃, H₂S, CO₂, and PH₃, in order to assess the impact of individual species on the detection. Only H₂O and CH₄ are robustly detected, and are plotted in Figure 3 along with the NH₃ and complete templates for comparison. The relative strength of the H₂O and CH₄ detections varies with simulated temperature, as expected from the equilibrium chemistry. The NH₃ template gives a featureless $\log RL(\rm K_{p})$ curve in the 1400 K and 1000 K models, but shows weak features in the 600 K case, though not significant in comparison with H₂O or CH₄. The H₂S template shows some weak features near the expected planet peak in the 1000 K models, but the detection is similarly not significant compared with H₂O or CH₄. The CO, CO₂, and PH₃ templates all produce featureless $\log RL(\rm K_{p})$ curves after correcting for the stellar contribution in all cases, consistent with the expected lack of $L$-band spectral features.

In the 600 K simulations, correlating with all molecules in the planet model results in a much stronger detection than would be expected from the individual H₂O and CH₄ cross-correlation functions. This does not appear to be the result of other molecules contributing substantially to the total detection. Rather, we note that omitting a molecule from the cross-correlation template – but not the simulated observation – effectively decreases the signal-to-noise of the observation. The features present in the simulated data but omitted from the template will still correlate with the template spectrum, but at values other than the planetary radial velocity, increasing the off-peak correlation and reducing the relative size of the planet peak. The impact is more pronounced than an increase in the Gaussian noise because the non-random structure produces off-peak structure that will add consistently when epochs are combined, whereas structure arising from random noise will average to zero across epochs.

Detection strength estimates for HRCCS detections have been made in a variety of ways. In NIRSPEC observations prior to the 2019 upgrade, Lockwood et al. (2014) reported a 6$\sigma$ detection in 5 epochs by comparing observed and synthetic log-likelihood functions. Piskorz et al. (2016) reported a strong detection based on 6 epochs by computing the Bayes factor, comparing a Gaussian planet detection with a linear nondetection in the log-likelihood space. Piskorz et al. (2017) reported a 3.7$\sigma$ detection in 7 epochs and Piskorz et al. (2018) reported a 3.8$\sigma$ detection in 6 epochs using a similar approach. Buzard et al. (2020) reported a 6.5$\sigma$ detection in 7 epochs, two of which were taken with the upgraded NIRSPEC, by fitting the full observed log-likelihood function with simulations.

Using the Bayes factor technique described in Piskorz et al. (2018) and the complete all-molecule planet template results in $>5\sigma$ detections for all 50-epoch simulations and all 25-epoch simulations except for two cases with $\rm T_{eq}=600$ K, $\rm R\leq 1.1R_{J}$. This is consistent with the higher resolution, increased spectral coverage, and increased number of epochs used in the simulations compared with prior NIRSPEC observations, as well as the correction for stellar correlation features clarifying the planet signal. The 10-epoch simulations range from $<2\sigma$ to $>5\sigma$, despite roughly the same total S/N as the 25-epoch simulations, indicating that many low-S/N observations yield more reliable detections than a smaller number of epochs with similar total S/N, consistent with the simulations from Buzard et al. (2020).

The correction of non-planetary correlation features allows us to compare the Bayes factor approach with a likelihood ratio, which we compute as:

Where $\log L(\hat{K}_{p})$ the the value of the log-likelihood function at the planet peak and $\log L(\bar{K_{p}})$ is the median value far from the planet peak. The resulting LR values are listed in Table 1. While these values can be used to compute a $\sigma$ detection confidence using Wilkes’ Theorem, doing so fails to account for the width of the planet feature in the log-likelihood space, and underestimates the detection confidence as a result. However, the likelihood ratio approach is easier to calculate for a large number of simulations, particularly for very well-detected cases, and still yields the correct relative detection strengths. These factors make the likelihood ratio a more useful measure of how detection confidence will vary with planet properties and observation strategy.

We assess the impact of planet equilibrium temperature on the cross-correlation detection by correlating 10 and 25 epochs of simulated observations with $\rm R=1.0\ R_{J}$, $\rm S/N_{epoch}=1500$, and varying planet equilibrium temperatures with models ranging from 600 K to 1400 K in steps of 100 K, covering the “warm/hot” range of exoplanet temperatures. This allows us to compute the two-parameter log-likelihood function $\log L(\rm K_{p},T_{eq})$, plotted in the left column of Figure 4 for $\rm T_{eq}=$ 1300 K, 1000 K, and 700 K. As before, we also compute the log-likelihood surface from a simulation of the stellar spectrum alone, which is then subtracted in order to remove structure arising from the correlation between stellar and planetary spectra (see equation 4). The resulting log relative likelihood surfaces, $\log RL(\rm K_{p},T_{eq})$, are plotted in the right column of Figure 4.

We then marginalize the two-parameter $\log RL(\rm K_{p},T_{eq})$ surface over the 10 km $\rm s^{-1}$ region surrounding the best-fit value of $\rm K_{p}$, roughly corresponding to the resolution of NIRSPEC, to obtain $\log RL(\rm T_{eq})$, plotted in black in the left panel of Figure 5. The size of the marginalization region has minimal impact on the final constraints obtained for simulations in which unwanted star/planet correlation features are effectively removed. We use a 10 km $\rm s^{-1}$ region for subsequent analysis in order to reflect the uncertainty in the best-fit $\rm K_{p}$ value while minimizing contamination from residual non-planetary correlation features. As expected based on the right column of Figure 4, the $\log RL(\rm T_{eq})$ curve is temperature-dependent but does not show a clear peak at the true value. We therefore use a modeling approach to attempt to constrain the temperature, analogous to the technique used by Buzard et al. (2020) to determine $\rm K_{p}$ in the presence of significant non-planetary features in the $\log RL(\rm K_{p})$ curve. In this case, $\rm K_{p}$ is already well-constrained (see Figure 4, right column), and we instead attempt model extraneous features in the $\log RL(\rm T_{eq})$ curve in order to better estimate $\rm T_{eq}$. A series of models is created with equilibrium temperatures from 600 K to 1400 K in steps of 100 K. The two-parameter log-likelihood surface is calculated for each model, and a simulated stellar spectrum is subtracted to remove off-peak correlation features. Finally, we marginalize the models over the same $\rm K_{p}$ range as was previously used for the “observed” spectrum. Both the “observed” slice we wish to constrain and these model slices are subtracted to have zero mean, accounting for the lack of Gaussian noise in the models, and plotted in the left panel of Figure 5.

While there is not a unique peak at the correct temperature, the $\log RL(\rm T_{eq})$ function still shows significant changes in its shape with model equilibrium temperature. Constraining the temperature therefore requires an additional step to compare the observed and expected $\log RL(\rm T_{eq})$. We compute the negative log-likelihood of the deviation between the “observed” curve and each model, under the assumption that the model points are Gaussian-distributed:

Where $O_{i}$ is the value of the observed log-likelihood function at equilibrium temperature $i$ in the parameter space and $M(\rm T_{eq})_{i}$ is the value of the model log-likelihood function with equilibrium temperature $\rm T_{eq}$ at temperature $i$. $\sigma(\rm T_{eq})_{i}$ is the error in the associated observed values determined from a jackknife test. To more easily compute best-fit values and confidence intervals, we convert the log-likelihood to a probability distribution $\rm P(\rm T_{eq})$ by exponentiating and normalizing the integral. This can be written as:

Where $C$ is a numerically-determined normalization constant so that $\rm P(T_{eq})$ integrates to unity. The value of $\rm T_{eq}$ which maximizes $\rm P(T_{eq})$ is the best-fit temperature of the “observed” planet and the 1$\rm\sigma$ errors on $\rm T_{eq}$ are computed as the interval enclosing an area of 0.68 around the best-fit value. These are reported in Table 2.

Meaningful constraints are obtained for $\rm T_{eq}\leq 1200$ K in both the 10 and 25 epoch cases. At higher temperatures, the fitting becomes unreliable, with the 1300 K and 1400 K cases effectively indistinguishable even in the 25 epoch simulations. The measured values are generally smaller than the true values, and an increased number of epochs results in negligible improvements to accuracy. Additionally, the reported precision based on jackknife estimates of the error in $\log RL(\rm T_{eq})$ appears to underestimate the true uncertainty. Comparing the true and measured values suggests the precision is $\sim$50–100 K, though the 100 K resolution of the model grid limits our ability to make a robust estimate.

At the same time, Figure 4 shows that the correct value of $\rm K_{p}$ is recovered even with significant mismatches in $\rm T_{eq}$ between observed and template spectra, so strong priors are not required in order to make an initial detection. While neglecting the day/night temperature mismatch will have a negative impact on the detection strength compared with models accounting for the full temperature structure, it does not produce an error in $\rm K_{p}$. Finally, we note that the detection strength is maximized when the cross-correlation template is cooler than the observed spectrum by $\sim 200$ K.

Piskorz et al. (2018) attempted to place a C/O ratio constraint on the hot Jupiter KELT-2Ab using 2D multi-epoch cross-correlation with data from Spitzer and pre-upgrade Keck-NIRSPEC. The large errors in this measurement are consistent with Figure 1 and Figure 3, which for a hot planet indicate the $L$-band spectrum and resulting cross-correlation function is dominated by H₂O features, with no measurable contribution from carbon-bearing species. The lack of carbon species lines at high temperature prevents measurement of the C/O ratio. However, the simultaneous detection of H₂O and CH₄ at equilibrium temperatures of approximately 1000 K in Figure 3 suggests a meaningful C/O constraint may be possible for cooler planets at the same wavelengths, assuming chemical equilibrium and known metallicity.

To assess the ability to constrain C/O in cooler planets, we simulate 25 epochs of planets with an equilibrium temperature of 900 K, radius 1 $\rm R_{J}$, $\rm S/N_{epoch}=1500$, and a range of C/O ratios, analogous to the simulations in Section 3.2 but varying C/O ratio instead of $\rm T_{eq}$. We perform a set of simulations with solar metallicity and another with $10\times$ solar metallicity to estimate the impact on the C/O constraint. The left column of Figure 6 plots the two-parameter log-likelihoods, $\log L(\rm K_{p},C/O)$ for simulations with C/O of 0.8, 0.5, and 0.2 at solar metallicity. As with the two-parameter temperature likelihood functions, a star-only simulation was subtracted to reduce structured off-peak correlation, producing the $\log RL(\rm K_{p},C/O)$ surfaces plotted in the right column of Figure 6 (see equation 4). The correct value of $\rm K_{p}$ is recovered regardless of C/O ratio, indicating that a even a large mismatch between template and observed C/O ratio will not result in a non-detection. In contrast with the $\log RL(\rm K_{p},T_{eq})$ surfaces, the $\log RL(\rm K_{p},C/O)$ surfaces appear to be uniquely peaked near the correct value of C/O, after correcting for features arising from the stellar spectrum.

As was the case with equilibrium temperature, the two-parameter likelihood surface changes with model C/O ratio, indicating the multi-epoch technique is sensitive to observed C/O, and the recovery of the correct $\rm K_{p}$ at all C/O values suggests strong priors on the C/O ratio are not required to make a detection. We use the same approach to characterize the C/O sensitivity as was previously used for equilibrium temperature. We begin by marginalizing the 2D $\log RL(\rm K_{p},C/O)$ surface over the 10 km $\rm s^{-1}$ region surrounding the best-fit $\rm K_{p}$, roughly corresponding to the resolution of NIRSPEC, to obtain $\log RL(C/O)$, which we plot in the left panel of Figure 7.

While the $\log RL(C/O)$ curves are generally peaked near the correct value, the modeling approach described in the previous section continues to provide more accurate constraints when the models closely replicate observations. We therefore correlate a set of noise-free models which we fit to the observed $\log RL(C/O)$ curve following equation 6, which we then convert to $\rm P(C/O)$ following equation 7. The right panel of Figure 7 plots the resulting $\rm P(C/O)$ function for a simulation with $\rm C/O=0.5$, analogous to the right panel of Figure 5.

Table 3 shows that a good constraint on C/O is achieved for input values below C/O = 0.7 from 25 epochs of $\rm S/N_{epoch}=1500$ and solar metallicity, and below C/O = 0.9 for $10\times$ solar metallicity. At higher C/O ratios, the planet spectra become too similar to reliably distinguish between different values for C/O from the $L$-band spectrum alone. At lower C/O ratios, the measured values tend to overestimate the true C/O, particularly in the higher-metallicity models. Despite this, the measured C/O values are accurate to better than 0.1, which should enable differentiation between substellar, approximately stellar, and superstellar C/O ratios. Higher metallicity leads to somewhat better constraints, particularly at high C/O ratios, as the increased heavy element content leads to stronger, better-detected H₂O and CH₄ lines. Both accuracy and precision are substantially improved in the 25-epoch simulations compared with 10-epoch simulations of similar total signal-to-noise.

Our final set of simulations assesses the role of spectral grasp and spectral resolution in both detection and atmospheric characterization. In order to compare instrumental properties directly, without considering the impact of the significant absorption features present in the $L$-band, we use “space-like” simulations in which ixels affected by saturated tellurics are not masked, unlike the simulations presented previously. We consider three cases with ten evenly-spaced epochs each. First, we consider a space-based NIRSPEC2 analog. The $\rm S/N_{epoch}$ is 2500, spectral resolution is $R\sim 35000$ and the spectral coverage is plotted in Figure 1. Second, we consider a space-based instrument with the same wavelength range as NIRSPEC2, but with double the resolution ($R\sim 70000$). To account for the increased dispersion, $\rm S/N_{epoch}$ is reduced to 1768. Finally, we consider a space-based instrument with the $R\sim 35000$ and $\rm S/N_{epoch}=2500$, but with the spectral grasp doubled to provide nearly continuous coverage from 2.9–3.8 $\mu$m. To maintain a consistent sampling of the line spread function, the latter two cases were run with 4096 pixels instead of 2048 for NIRSPEC2.

For each set of instrument parameters, we simulate a set of planets with $\rm T_{eq}=1400$ K and radii from 0.7 $\rm R_{J}$ to 1.3 $\rm R_{J}$, similar to Section 3.1, as well as a set with $\rm T_{eq}=900$ K, $\rm R=1.0\ R_{J}$, solar metallicity, and a range of C/O ratios, similar to Section 3.3. Following the same procedures as in those sections, we present detection strengths as likelihood ratios in Table 4 and C/O constraints in Table 5.

Examining first the impact on detection strength, we find that both grasp and resolution lead to a significant improvement in detection. This can be seen in Figure 8, which plots the relative likelihood curves for each set of instrumental factors for 10 simulated epochs with $\rm T_{eq}=1400$ K, $\rm T_{eq}=1000$ K, and $\rm T_{eq}=600$ K. Comparing with the $\rm T_{eq}=1400$ K, 10 epoch case in Table 1, all detections listed in Table 4 are significantly stronger than the corresponding entries in Table 1, indicating that telluric losses are significantly degrading ground-based observations. Both increased spectral resolution and increased wavelength coverage offer additional improvements in detection strength, with the double resolution case producing a $\sim 60$ percent improvement in the likelihood ratio compared with the “space-like” NIRSPEC2 and the double grasp case yielding a $\sim 70$ percent improvement. The improvement in detection strength appears to be somewhat dependent on planet temperature. While for the $\rm T_{eq}=1400$ K and $\rm T_{eq}=1000$ K planet models both increased resolution and increased grasp produce similar improvements in detection strength, the $\rm T_{eq}=600$ K models are much better detected in the increased grasp case.

Table 5 shows the impact of instrument parameters on the C/O measurement. While increased spectral range appears to offer some improvements in accuracy, particularly at high C/O, increased resolution does not have a significant impact on the measurements obtained. We can also estimate the impact of telluric features on the achievable C/O constraint by comparing the space-like NIRSPEC2 case in Table 5 with the values in Table 3. Both accuracy and precision are substantially worse when saturated telluric absorption features are removed, particularly at higher C/O, consistent with the large loss in effective spectral grasp due to telluric features.

Simulations of identical systems with varying number of epochs allow us to directly compare the results of different observation strategies. Table 1 lists detection strengths for simulations with 10 epochs at $\rm S/N_{epoch}=2500$ as well as simulations of 25 epochs with $\rm S/N_{epoch}=1500$. The total signal-to-noise ratio in these simulations combining all epochs is $\sim$7900 and $\sim$7500 respectively. Despite the slightly better total signal-to-noise, the 10-epoch case gives a substantially worse detection strength in all cases compared with the otherwise-equivalent 25-epoch simulations, corresponding to a factor of $\sim 2$ reduction in photometric contrast. Comparing the 25-epoch simulations with 50-epoch simulations, which have 40 percent greater total signal-to-noise, the 50 epoch case performs substantially better that would be expected from the increase in signal-to-noise, again corresponding to a factor of $\sim 2$ in photometric contrast. This suggests that confidence in the multi-epoch detection is much more dependent on the number of epochs combined than on either the per-epoch or total signal-to-noise. This is consistent with findings from Buzard et al. (2020) that spreading the same total signal-to-noise over an increasing number of epochs results in a stronger detection.

The improvement with increasing number of epochs suggests the optimal observing strategy for multi-epoch cross-correlation is to take a large number of relatively low signal-to-noise observations, similar to stellar radial velocity surveys. While stellar optical spectra can be corrected to the required precision for RV measurements by dividing by the spectrum of an A0 telluric standard star, the extremely low planet flux measured in multi-epoch cross-correlation would require a prohibitive amount of time spent telluric standards in order to prevent the uncertainty in the standard measurement from dominating over the planet signal.

Instead of telluric standards, multi-epoch cross-correlation observations since Lockwood et al. (2014) have made use of line-by-line atmospheric models such as RFM (Dudhia, 2017) or MOLECFIT (Kausch et al., 2014) to fit and divide out a model atmosphere. Beginning with Piskorz et al. (2016), this is followed by a principal component analysis (PCA) to correct errors in the model line profile function, changes in molecular abundances over the observation, and instrument flexure (Piskorz et al., 2016, 2018). The guided PCA requires a series of observations at varying airmass in order to identify residuals associated with tellurics, as these features should vary consistently with the airmass. Based on prior NIRSPEC observations, a time series covering at least $\sim 30$ minutes per target is necessary for guided PCA to be effective. This in turn requires the use of the few-epoch, high per-epoch signal-to-noise observing strategy simulated in the 10 epoch case when observing with large telescopes, with $\rm S/N_{epoch}\sim 2500$. The results of the simulations presented here suggest that the development of alternative telluric correction procedures which do not impose a minimum integration time should be a high priority for future multi-epoch cross-correlation pipelines. We also note that more sensitive instruments may shorten the time series required for effective PCA, but will also further reduce the integration time required to achieve the desired $\rm S/N_{epoch}$ for the many-epoch observing strategy.

Table 1 allows us to compare the impact of photometric contrast changes due to differences in equilibrium temperature with changes due to planet radius. Comparing otherwise-identical simulations with $\rm T_{eq}$ of 600 K and 1400 K, we find the likelihood ratio is only reduced by $\sim 20$ percent in the 600 K case, despite a 94 percent reduction in $L$-band photometric contrast compared with the 1400 K case. At the same time, a reduction in photometric contrast of $\sim 20$ percent due to planet radius results in a $\sim 20$ percent reduction in likelihood ratio at fixed equilibrium temperature, indicating that photometric contrast has a significant impact on detection at fixed equilibrium temperature. We note that these numbers assume cloud-free models for both temperatures, which is not necessarily appropriate for the 600 K case. The impact of clouds and hazes is discussed in more detail below. Additionally, we do not account for changes in the planet spectrum with varying radius. Changes in surface gravity and temperature–pressure profiles in smaller planets are likely to introduce changes in spectroscopic contrast in addition to the change in photometric contrast considered in these simulations.

The difference in the relationship between photometric contrast and detection strength from changes in radius compared with changes in equilibrium temperature emphasizes that multi-epoch cross-correlation is much more sensitive to spectroscopic contrast than photometric. On a per-photon basis, the $\rm T_{eq}=600\ K$ planet models are significantly easier to detect via multi-epoch cross-correlation, due to the presence of deep CH₄ features at cooler temperatures. However, the presence of additional species in the spectra of cooler exoplanets increases the difficulty of accurately matching the correlation template to the observed spectrum, a challenge which is not present in the simulations. This suggests anticipating the detectability of exoplanet systems with cross-correlation spectroscopy requires an understanding of the target’s thermochemical properties, in particular considering the potential presence and impact of clouds/hazes. Figures 4 and 6 indicate a detection can be made despite significant differences between the template and observed spectra, provided the template includes the species present in the observed spectrum. An initial detection with a poorly-matched model can be subsequently revised by varying the template to maximize detection strength.

Section 4.2 found that the 25-epoch simulations detected nearly all simulated planets with $>5\sigma$ confidence using the Bayes’ factor approach described in Piskorz et al. (2016). This suggests that using a many-epoch, low per-epoch signal-to-noise observing strategy, $L$-band multi-epoch cross-correlation using NIRSPEC is sensitive to planets with $\rm T_{eq}\geq 600K$ and $\rm R\geq 0.7\ R_{J}$ around a Sun-like host star. Assuming a Sun-like host star, this corresponds to planets within $\sim 0.2$ AU approximately Saturn-size and larger, which should enable the detection and characterization of warm Jupiters. While we do not simulate the change in orbital period with semi-major axis, this does not affect the detection provided the sampling of the orbital phase is fixed. We caution that the exact limits are likely to be strongly dependent on the exact wavelength range observed, the planet chemical composition, the presence of clouds/hazes, and how well features arising from the stellar spectrum can be removed from the $\log L(\rm K_{p})$ space.

Figure 4 indicates that while mismatches between the true temperature and the correlation model negatively impact the detection strength, the best-fit value of $\rm K_{p}$ is independent of the model equilibrium temperature over an 800 K range. This is particularly relevant to planning observations for tidally-locked planets such as hot Jupiters, which are known to have day/night temperature differences of several hundred Kelvin (e.g. Knutson et al., 2007; Wong et al., 2016; Komacek & Showman, 2016). Figure 4 indicates both dayside and nightside observations can be used effectively in the absence of longitudinally-varying clouds/hazes with a relatively minor impact on the detection strength, though the correlation model should use the lowest estimated temperature. Underestimating the planetary equilibrium temperature in the template spectrum consistently results in a better detection.

While the weak dependence on mismatches between the observed and model temperature is useful observationally, it suggests $L$-band observations have a limited ability to constrain thermal properties of the system. While Figures 4 and 5 do show variations with equilibrium temperature, the values in Table 2 show significant systematic errors. Additionally, the temperature constraint in the $L$-band arises primarily from the relative strengths of H₂O and CH₄ features, which will also depend on metallicity and C/O ratio. Multi-dimensional model grids to assess these degeneracies are computationally impractical at present, and will require improvements in the cross-correlation routine. The $L$-band temperature sensitivity is therefore unlikely to be useful in practice beyond providing rough constraints on the day/night temperature contrast for $\rm T_{eq}\leq 1200$ K, when CH₄ features are present in the spectrum. However, other NIR spectral features with stronger temperature dependencies, such as CO bandheads in the $K$ and $M$-bands, may offer much better constraints on 3D thermal properties, as was demonstrated in Beltz et al. (2020). With a large number of epochs, variations in the best-fit temperature between epochs at different orbital phase could enable thermal mapping of non-transiting planets, similar to the phase-curve mapping technique in Knutson et al. (2007). Future simulations will explore the wavelength-dependence of thermal constraints from multi-epoch cross-correlation and the impact of phase-dependent observed planet temperatures. We expect the temperature-dependence of the planetary spectrum will lead to constraints on planet properties varying strongly with equilibrium temperature, even for observations at fixed wavelength.

Figure 6 indicates that while the best-fit $\rm K_{p}$ is unaffected by a mismatched C/O ratio, multi-epoch cross-correlation successfully recovers the true value for $\rm T_{eq}=900$ K and C/O $\leq 0.8$, as shown in Figure 7 and Table 3. Higher metallicity offers marginally better performance at high C/O. The contrast with the weak C/O constraint obtained in Piskorz et al. (2018) is due to the presence of CH₄ in the $L$-band spectrum for $\rm T_{eq}\leq 1000$ K, which allows carbon and oxygen bearing species to be observed simultaneously. This again illustrates the dependence of multi-epoch cross-correlation detection and characterization on the wavelength range observed and the target’s atmospheric composition. While $L$-band observations can provide useful C/O constraints for warm planets, hot Jupiters require observations at different bands which contain carbon species (e.g. CO) features at high temperature. Useful C/O constraints are not obtained for simulated planets with $\rm T_{eq}\geq 1000$ K due to the lack of CH₄ lines at high temperatures, though for cooler planets differences in $\rm T_{eq}$ between the observed and template spectra do not significantly impact the C/O constraint beyond increasing the uncertainty to $\sim 0.1$.

We note that $L$-band C/O constraints require an equilibrium chemistry assumption, as only H₂O and CH₄ features are observed. Wallack et al. (2019) found that disequilibrium processes in cooler plants have a minor impact on molecular abundances compared with atmospheric metallicity and C/O ratio for planets with equilibrium temperatures around 1000 K, suggesting the equilibrium assumption is reasonable for the warm exoplanet populations. In cases where disequilibrium processes are expected to play a significant role, a CO detection is likely to be required, and disequilibrium processes must be included in the model planet spectrum. Large-grasp instruments which can detect multiple carbon and oxygen-bearing species simultaneously would allow large-scale disequilibrium processes to be identified through cross-correlation without a priori knowledge of the atmosphere.

The ability to measure C/O from $L$-band observations has the potential to clarify the formation processes for warm Jupiters. Öberg et al. (2011) found the C/O ratio of the gaseous and solid components of protoplanetary disks varies with distance. Gas beyond the H₂O snow line is carbon-enriched relative to the inner disk, while solids are carbon-depleted. If the gas accretion determines the final C/O ratio, planets which form beyond the snow line and migrate inwards should have significantly superstellar C/O (e.g. Öberg et al., 2011; Öberg & Bergin, 2016; Madhusudhan et al., 2017). Other work has suggested solid accretion dominates the final atmospheric C/O (e.g. Mordasini et al., 2016; Espinoza et al., 2017), resulting in substellar C/O values which drop further beyond the H₂O snow line. The relative importance of gas and solid accretion may also vary significantly with planet mass, leading to compositional differences between Jupiter and Neptune-size planets (Cridland et al., 2019). The precision of the C/O constraints in Table 5 suggest cross-correlation spectroscopy can distinguish between migration and in-situ formation scenarios, as well as between gas and solid-dominated accretion. The weak inclination dependence of cross-correlation techniques results in many more targets for which C/O measurements can be made compared with transit spectroscopy, increasing the potential sample size for studies of atmospheric C/O in intermediate-temperature ($\rm T_{eq}\sim 900$ K) planets.

These simulations suggest multi-epoch cross-correlation is capable of significantly more than has been demonstrated observationally. In part, this is because these simulations are based on the capabilities of NIRSPEC2 (Martin et al., 2018), which offers substantially improved wavelength coverage and resolution compared with the pre-upgrade instrument which was the basis for prior multi-epoch detections. As additional observations are taken with the upgraded instrument, future observational results should better match the simulations presented here. We also briefly note several factors which were not included in the simulations which impact observations. Stellar activity and atmospheric clouds/hazes are discussed in detail below.

Figures 4 and 6 suggest minor inaccuracies in the planetary spectral template temperature or C/O ratio will have minimal impact on the ability to make a detection. While inaccuracies decrease the strength of the detection, the correct value of $\rm K_{p}$ is recovered for $\rm T_{eq}$ between 600 K and 1400 K and C/O between 0.1 and 1.1, regardless of the true values. Furthermore, we demonstrate that correlating with a grid of planet models can provide constraints on planet properties which are not known a priori. The one case where differences between observed and template spectra causes significant issues with detection is when simulated observations with $\rm T_{eq}<1000$ K are correlated with $\rm T_{eq}\geq 1000$ models. This appears to be due to the absence of CH₄ features in warmer models. Provided the species which dominate the observed spectrum are represented in the cross-correlation template, errors in the template appear to have a minimal and identifiable impact on the detection. Errors in line shape or position in the template will lead to a broader, lower-amplitude peak in the cross-correlation function, and may lead to shifts in the best-fit $\rm K_{p}$ (Brogi & Line, 2019). We note that errors in linelists represent a significant additional source of modelling uncertainty, particularly for CH₄. Using an updated linelist, Gandhi et al. (2020) is unable to reproduce the CH₄ detection reported in Guilluy et al. (2019) for HD 102195 b, indicating the uncertainties in current linelists are large enough to impact both detection and C/O constraints with cross-correlation techniques. Future improvements in high-resolution linelists should reduce the impact of this model uncertainty.

These simulations also lack any residuals from the telluric removal procedure, unlike observations. Buzard et al. (2020) showed that such residuals are not necessary to reproduce major non-planetary features in the cross-correlation space, which arise from correlation between the stellar and planetary spectra. While the guided PCA approach should effectively remove all features present in a model telluric spectrum or which vary consistently over the observation, fixed differences between the observed line profiles and positions and the telluric model are likely to lead to residuals in observed spectra which are not included in this simulation framework, but which do contribute to off-peak structure in observations. Developing simulations which accurately include telluric residuals is an ongoing challenge, as it requires quantifying the deviations between observed and model telluric spectra and how such differences vary with airmass and observing conditions. The near-perfect correction of off-peak structure achieved in these simulations therefore represents a best-case scenario for telluric removal, and we emphasize the importance of accurate telluric correction for successful multi-epoch cross-correlation.

Table 4 and Figure 8 indicate improvements to either spectral resolution or spectral grasp lead to significant improvements in detection strength. While the increased coverage offers a somewhat greater performance improvement compared with increased resolution for warmer equilibrium temperatures, improved spectral grasp offers substantially better performance than increased spectral resolution for the $\rm T_{eq}=600$ K case. We believe the relative importance of resolution and grasp is likely to depend on both the observed wavelength and properties of the target planet, in particular the width of the observed planetary emission features. We leave a more detailed exploration of instrumental factors in detection strength considering a broader portion of the near-infrared spectrum to future work. Table 5 shows the atmospheric characterization from $L$-band features is much less dependent on instrument properties. Improving the C/O constraint is likely to require the ability to simultaneously detect additional species beyond H₂O and CH₄, necessitating a significant increase in spectral grasp rather than additional improvements in spectral resolution.

Several instruments have been proposed or are being developed to offer single-shot 1–5 $\mu$m coverage with spectral resolution greater than Keck-NIRSPEC, including GMTNIRS, CRIRES+, and IGNIS. In addition to improving the detectability of planets, large simultaneous wavelength coverage is likely to offer significant improvements in the ability to characterize planets by enabling the simultaneous detection of additional molecules. The large spectral ranges of these instruments will be especially valuable in cases of significant non-equilibrium chemistry, clouds or hazes, or a prior uncertainty in the atmospheric composition.

Finally, we note that the space-like NIRSPEC2 simulations in Table 4 result in a factor of $\sim$5 improvement in the likelihood ratio compared with the analogous 10-epoch simulations in Table 1 which remove portions of the spectrum affected by telluric absorption. A space-based high-resolution spectrograph covering the 1–5 $\mu$m spectral range would offer significantly greater capability to detect and characterize exoplanets through cross-correlation spectroscopy compared with ground-based facilities.

These simulations used a Sun-like star for the host stellar spectrum. For a given planet radius and equilibrium temperature, the relative brightness of the planet increases with decreasing stellar radius and temperature. From a purely photometric viewpoint, we would therefore expect to make stronger detections around later-type stars, as well as the ability to detect smaller and cooler planets. However, the results from Section 4.2 indicate multi-epoch cross-correlation is much more sensitive to spectroscopic contrast than photometric. Increased stellar activity and stronger molecular lines in late-type stars may cause spectroscopic contrast to diverge significantly from photometric contrast, impacting planet detectability in ways that are difficult to predict without dedicated modeling. A rigorous assessment of cross-correlation techniques around K and M type primary stars is currently limited by the lack of sufficiently accurate high-resolution spectral models for late-type stars.

The development of improved models for late-type stars would further expand the population accessible to characterization through multi-epoch cross correlation with existing instrumentation. Of particular interest is planets falling in the radius valley identified by Fulton et al. (2017) near $\rm 1.8\ R_{\oplus}$. While our simulations show such planets are not likely to be detectable with Keck-NIRSPEC around Sun-like stars, the increased relative brightness of exoplanets around M-dwarfs may enable detections, provided accurate stellar models are available and the atmospheric spectrum offers sufficient spectroscopic contrast. Cooler planets may pose additional challenges for telluric correction as the emission spectrum of the planet becomes more similar to the terrestrial spectrum. The ability to probe the atmospheric composition of these planets with multi-epoch cross-correlation could offer a new avenue to investigate the evolutionary processes affecting intermediate-mass, highly irradiated planets near the transition between rocky and gaseous compositions.