Feasibility of High-resolution Transmission Spectroscopy for Low-velocity Exoplanets

For our present exploratory study, we consider simulated planet spectra injected into high-resolution time-series observations of a known sub-Neptune, TOI-732 c (LTT 3780 c). The observed spectra were obtained using the IGRINS spectrograph at Gemini South (Gemini-S) as part of GO Program GS-2021A-Q-201 (PI: D. Valencia; Cabot et al. 2024). We consider the observations of one transit obtained on the night of 2021 February 23. TOI-732 c was discovered by the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015) and followed up with ground-based radial velocity observations with CARMENES and HARPS (Cloutier et al. 2020; Nowak et al. 2020). The system parameters are shown in Table 1. For the transit duration there are multiple values reported in the literature: 1.4 hr (Cloutier et al. 2020) and 1.8 hr (Nowak et al. 2020). We adopt the latter value of 1.8 hr to be conservative. There have been no atmospheric characterization studies of this planet to date, although several recent studies have considered simulated studies with JWST (Madhusudhan et al. 2021; Constantinou & Madhusudhan 2022).

Table 1. System Properties of TOI-732 c

Parameter	Value
P	${12.252131}_{-0.000064}^{+0.000072}$ days
T0	${2458546.8492}_{-0.0017}^{+0.0016}$ BJD
Rstar	0.382 ± 0.012 R⊙
Rp	2.42 ± 0.10 R⊕
Vsys	0.1954 ± 0.0005 km s−1
a	0.0762 ± 0.0034au
i	${89.08}_{-0.13}^{+0.11}$ °
e	${0.115}_{0.065}^{+0.07}$
T14	1.79 ± 0.21 hr
Kp	70 ± 20 km s−1

Note. Values adopted from Nowak et al. (2020), except V_sys and K_p, which are from Cloutier et al. (2020). We note that Cloutier et al. (2020) gives a value for T₁₄ of ${1.392}_{-0.050}^{+0.049}$ hr, but here we use the more conservative value from Nowak et al. (2020).

Download table as:  ASCIITypeset image

IGRINS is a high-resolution spectrograph mounted on the 8 m Gemini-S telescope at Cerro Pachón, Chile (Yuk et al. 2010; Park et al. 2014). IGRINS observes the NIR spectral range, covering ∼1.45–2.45 μm, split into 54 spectral orders, and with a resolution of R ∼ 45,000. The present time-series observations include 33 spectra (A/B pairs) obtained over 2.8 hr. This includes 21 ± 2 spectra during the transit, calculated using the transit duration given in Table 1. The different value for T₁₄ in Cloutier et al. (2020) would instead suggest that ${15}_{-0}^{+2}$ of the spectra are in transit. The air mass varied between 1.05 and 1.18, and the magnitude of the barycentric velocity correction varied between 3.2 and 2.8 km s⁻¹, over the course of the observations.

We follow the methods in Cheverall et al. (2023) to clean the spectra of bad pixels and outliers, with each spectral order treated independently. We begin with the reduced spectra of the observations made available through the Gemini/IGRINS archive. These are the time-series spectra divided by the spectrum of an A0V telluric standard star obtained using the same setting. We remove 24 orders (0–4, 24–37, and 49–53 inclusive) due to poor-quality data or lack of a well-defined continuum. This leaves 30 spectral orders for the remainder of the analysis. To normalize the spectra, we fit a second-order polynomial to the continuum of each spectral order and exposure. However, due to the relatively low temperature of the M-dwarf host star, stellar molecular features appear in the spectral continuum which in some cases are difficult to remove with a fitted polynomial. Despite this, with the exception of H₂O, most of the molecules typical of temperate planets are not expected to be present in the stellar photosphere, and therefore any uncorrected stellar spectral features should not interfere with our injection and recovery tests, the results of which we compare with and without injections. A more suitable normalization function would be desirable for real chemical inferences, however.

For ground-based observations of spectra in the NIR, the planet transmission signal is embedded in contaminating telluric and stellar features as well as correlated noise from other sources, such as the instrument. In order to isolate the planet signal from such contaminants, which are significantly stronger, we detrend our spectra using PCA, as is commonly done in other works (de Kok et al. 2013; Giacobbe et al. 2021; Lafarga et al. 2023). This number of principal components subtracted, or "PCA iterations applied," is calculated by optimizing the S/N from a noise-subtracted, differential cross-correlation function (CCF), ΔCCF (Holmberg & Madhusudhan 2022; Spring et al. 2022), in order to avoid the introduction of an optimization bias toward the signal for which we are searching (Cheverall et al. 2023). This is as opposed to optimizing for the detection S/N directly, which has the potential to introduce bias (Cabot et al. 2019; Spring et al. 2022; Cheverall et al. 2023). We calculate the optimum PCA iteration for all of the three models considered here, and find optimal values of three, five, and seven for NH₃, CH₄, and H₂O, respectively.

The planetary signal, which is quasi-stationary in wavelength, constitutes a lower-frequency signal in the time-series spectra than a typical high-velocity planet signal, which is shifting across the instrument pixels throughout the observing night. We find that the degradation of the planet signal is therefore significantly faster at lower velocities (Figure 1). It is hence reasonable to expect that fewer principal components may need to be subtracted in order to maximally recover low-velocity planet signals, since the subtraction of too many components will remove the planet signal alongside the quasi-stationary telluric and stellar lines.

Zoom In Zoom Out Reset image size

To compute model templates for the transmission spectra of the sub-Neptune TOI-732 c considered in this study, we follow a similar approach to that of Cheverall et al. (2023), which we briefly reiterate here. The transmission spectrum is modeled using the AURA atmospheric modeling code for exoplanets (Pinhas et al. 2018). The model computes line-by-line radiative transfer in transmission geometry assuming a plane-parallel H₂-rich atmosphere in hydrostatic equilibrium, over a pressure range of 10⁻⁷–100 bar. The chemical composition and temperature structure are free parameters in the model. We generate the spectra considering one molecule at a time, assuming no clouds/hazes and an isothermal temperature profile. We consider a nominal temperature of 300 K, close to the equilibrium temperature of TOI-732 c (Madhusudhan et al. 2021). The spectra are computed at high resolution (R ≳10⁵) in the NIR spanning the IGRINS spectral range (1.4–2.5 μm).

We consider opacity contributions due to prominent molecules expected in temperate H₂-rich atmospheres (H₂O, CH₄, and NH₃) and assuming a nominal mixing ratio of 10⁻⁴ for each molecule. The molecular cross sections were obtained following the methods of Gandhi & Madhusudhan (2017) using absorption-line lists from the following sources: H₂O (Barber et al. 2006; Rothman et al. 2010), CH₄ (Yurchenko & Tennyson 2014), and NH₃ (Yurchenko et al. 2011). We also include collision-induced absorption from H₂–H₂ and H₂–He (Borysow et al. 1988; Orton et al. 2007; Abel et al. 2011; Richard et al. 2012). As with the observed spectra, the model spectra are separated into orders and then normalized using a polynomial fit. This is then followed by convolution with the point-spread function of the instrument. While the absolute depths of absorption lines are lost after normalization, their relative depths and positions are conserved.

After detrending the spectra of unwanted telluric and stellar features, as described in Section 2.3, we are left with residuals containing the remaining planet signal. Each individual line in this planetary spectrum has S/N < < 1, so information from each needs to be combined in order for the signal to be characterized. Therefore, the residuals are typically cross-correlated with a Doppler-shifted model, shifted as a function of planetary velocity, to give a CCF as a function of radial velocity and orbital phase. This planetary radial velocity can be given by

$$\begin{eqnarray}&&{V}_{{\rm{p}}}={K}_{{\rm{p}}}\sin (2\pi \phi )+{V}_{\mathrm{sys}}-{V}_{\mathrm{bary}}+{V}_{\mathrm{offset}},\end{eqnarray} \tag{ 1 }$$

where K_p is the semi-amplitude of the planet's orbital motion, V_sys is the systemic velocity of the planetary system, V_bary is the barycentric velocity correction where a positive value indicates the Earth is moving toward the star, and V_offset accounts for any radial velocity offset, originating from sources such as planetary atmospheric winds.

The standard approach to extract a signal from the CCF is to calculate a detection S/N for each point in K_p−V_sys space. The signal is calculated by summing the CCF values along a trail defined by a given point in K_p−V_sys space, and then divided by a noise estimate obtained from the standard deviation of CCF values away from this trail, as described in Cheverall et al. (2023). A strong cross-correlation signal at the known orbital parameters of the planetary system would indicate a detection. The left-hand panel of Figure 2 shows the resulting S/N maps when the detrended residuals are cross-correlated with NH₃, CH₄, and H₂O models. There are no significant (S/N > 3) cross-correlation signals recovered at the planetary velocities at which we later inject model planet spectra (marked by crosshairs) to create synthetic data sets in Section 3, with only a tentative feature seen at V_sys ∼20 km s⁻¹ for NH₃.

Zoom In Zoom Out Reset image size

These are the basic techniques with which we complete injection and recovery tests in Section 3. This is then repeated using more in-depth Bayesian analysis of the spectra in Section 4.

In order to investigate the limits of high-resolution transmission spectroscopy for long-period planets, we simulate the spectra of such targets. Previous works using simulations of high-resolution spectra have often involved the generation of a synthetic data set, built from models of the component parts (Gandhi et al. 2020; Hood et al. 2020; Genest et al. 2022). In this work, rather than generate a synthetic data set, we use the real observations described in Section 2.1 and inject a high-resolution model planetary signal. This has the benefit of ensuring we inject synthetic planetary spectra into data sets containing realistic and comprehensive sources of noise. It does, however, limit our control and knowledge over the noise in the data.

Along with the telluric absorption and correlated noise from the instrument, the broadband stellar features present in the spectral observations of an M-dwarf star, as discussed in Section 2.2, also modulate the planetary spectral lines. To simulate this modulation, we initially inject our planetary spectra prior to normalization.

When simulating a data set, we must select which of the spectra will be in transit, i.e., the spectra into which we inject a model planet spectrum, and which will be out of transit. This selection is initially made symmetrically in phase as follows: For N in-transit and M out-of-transit spectra, the N spectra with the N smallest values of ∣ϕ∣ become the in-transit spectra, while the next smallest M exposures become the out-of-transit spectra for use in detrending. Any remaining exposures (in the case where N + M <33) are removed from the data set in this instance. Note that our phases are near-symmetrical about zero, with the central exposure near mid-transit, meaning that when an odd number of in-transit spectra are included in our simulated data set, there are an equal number of out-of-transit spectra on each side of the transit. Partial transits are investigated in Section 6.

We now complete injection and recovery tests with synthetic atmospheric model signals of NH₃, CH₄, and H₂O, which represent a range of telluric contributions. These models are individually injected, following the steps in Section 3.1, with the expected planetary parameters given in Table 1 (K_p of 70 km s⁻¹, with 21 in-transit and 12 out-of-transit spectra). We aim to recover these injected signals as described in Section 2, subtracting three, five, and seven principal components, respectively, to detrend in accordance with the ΔCCF metric.

The cross-correlation results are shown in Figure 2, with detection significances of 5.0, 4.2, and 2.1 found for NH₃, CH₄, and H₂O, respectively. The cross-correlation signals shown do not coincide with any of the positive background noise seen without injection; since the only difference between these two cases is the injection of a planetary signal prior to detrending, it is likely that the signals seen here are indeed recoveries of those injections. Similar injection and recovery tests are completed for the high-velocity case used in Figure 1, with results shown for comparison in Figure 12.

We find that we are able to recover these injections independently of the systematic velocity. This is demonstrated in Figure 2, where the value of V_sys at which each model is injected and recovered is varied for completeness, and the recovery of NH₃ at a V_sys of 0 km s⁻¹ eliminates a large V_sys as responsible for the separation of planet and telluric spectra. In order to test the robustness of the recovered signals, we verify that the detection S/N is somewhat stable against the number of principal components removed in detrending.

We additionally repeat this using the raw spectra without the A0V correction used here (Figure 13), and find that we continue to recover the low-velocity planet signals with S/N values consistent with before (5.0 for NH₃, 3.4 for CH₄, and 1.6 for H₂O). We also note that similar results are found when the signals are injected after normalization.

Our inability here to recover a significant (S/N > 3) H₂O signal may perhaps be due to the weaker spectral features of this molecule compared to NH₃ and CH₄. The signal strengths shown in Figure 1 are normalized by their initial values at zero PCA iterations, with the initial H₂O signal significantly weaker than that for NH₃. The reduced detection significance for this molecule is also seen when simulating a high-velocity injection (Figure 12), and no significant increase in S/N is observed when the signal is injected and recovered at artificially high systematic velocities, e.g., >300 km s⁻¹ , to decouple it from telluric features. The H₂O signal with features amplified by 2.5× is detected with S/N > 4 in Figure 4.

The relatively small phase range and radial velocity semi-amplitude of TOI-732 c means that the radial velocity change during the transit is very small, less than the velocity resolution of the instrument (see Table 2). However, we have shown that we are still able to recover injected signals of varying strengths for this planet using PCA-based detrending. We now probe the limits of low-radial-velocity-change planets, by injecting a model of NH₃ into the spectra at K_p = 0 km s⁻¹ while simulating zero barycentric velocity correction, and once again attempting to recover. After again subtracting three principal components in detrending, a signal with a S/N of 4.7 is found (Figure 3), despite there being no change in the radial velocity of the planet during the transit, i.e., the planet signal undergoes no wavelength shift relative to the telluric contaminants.

Zoom In Zoom Out Reset image size

Table 2. Typical Changes in the Radial Velocity of Planets During Their Transits

Planet	Kp	${\phi }_{\max }$	Δvp	References
	(km s−1)		(km s−1)
HD 189733 b	152.5	0.017	33	B16; A10; A19
HD 209458 b	145	0.018	33	G21; K07; A12
WASP-76 b	196.5	0.044	107	E20
TOI-270 d	70	0.004	3.5	G19; VE21
TOI-732 c*	70	0.003	2.6	N20; C20

Notes. For reference, a typical instrument pixel velocity V_px is around 1–2 km s⁻¹ for resolutions of R ∼ 50,000–100,000. We here compare the hot Jupiters studied in Cheverall et al. (2023; HD 189733 b, HD 209458 b, and WASP-76 b) to example sub-Neptunes (TOI-270 d and TOI-732 c). * indicates planet studied in this work. Note that contributions from the changing barycentric velocity correction are not considered here, and errors are not given as we are simply indicating the approximate pixel shifts for different regimes of planet. References: B16 (Brogi et al. 2016); A10 (Agol et al. 2010); A19 (Addison et al. 2019); G21 (Giacobbe et al. 2021); K07 (Knutson et al. 2007); A12 (Albrecht et al. 2012); E20 (Ehrenreich et al. 2020); G19 (Günther et al. 2019); VE21 (Van Eylen et al. 2021); N20 (Nowak et al. 2020); C20 (Cloutier et al. 2020).

Download table as:  ASCIITypeset image

We have shown that little-to-no relative change in the Doppler shift of the planet signal compared to that of the telluric spectrum does not prohibit PCA-based detrending. Therefore, it is possible to chemically characterize the atmospheres of temperate, longer-period planets, and the assumption that the target requires a distinguishable change in radial velocity during the transit is not necessarily true. However, despite being able to recover the signals, we are unable to constrain any information about the K_p of the planet in any of the cases shown in Figures 2 or 3.

For the methods described in Section 2.5, in order to determine the K_p of a cross-correlation signal we must measure the change in radial velocity of the planet during the course of the transit, Δv_p. This measurement is converted to K_p, assuming an accurate value for the maximum phase during transit, ${\phi }_{\max }$, for comparison with the known orbital parameters of the planetary system. ${\phi }_{\max }$ is calculated from the transit duration, T₁₄, and the period, P. We are, however, limited by the resolution of the spectrograph instrument, with the width of each wavelength channel corresponding to a given velocity shift, which we will henceforth refer to as the instrument pixel velocity, V_px. This propagates an intrinsic error into measurements of Δv_p and hence K_p. The number of instrument wavelength channels by which the planetary signal is Doppler shifted during the transit, which we define as the planet pixel shift, N_px, is given by

$$\begin{eqnarray}&&{N}_{\mathrm{px}}=\displaystyle \frac{{\rm{\Delta }}{v}_{{\rm{p}}}}{{V}_{\mathrm{px}}}=\displaystyle \frac{2{K}_{{\rm{p}}}\sin (2\pi {\phi }_{\max })-{\rm{\Delta }}{v}_{\mathrm{bary}}}{{V}_{\mathrm{px}}}.\end{eqnarray} \tag{ 2 }$$

The radial velocity changes of example exoplanets during their transits are shown in Table 2. The IGRINS (R ∼ 45,000) data we use here have an instrument pixel velocity of V_px ≈ 2 km s⁻¹. Typical hot Jupiters, such as those given in Table 2, may therefore have planetary pixel shifts N_px of ∼10–100 for this instrument. Other high-resolution spectrographs with resolutions R ∼ 50,000–100,000 may increase this pixel shift for a given planet. For example, CARMENES (R ∼ 80,000) observations of WASP-76 b (Landman et al. 2021; Sánchez-López et al. 2022) correspond to a planetary pixel shift of N_px ∼ 80, while CARMENES observations of HD 189733 b (Alonso-Floriano et al. 2019) correspond to N_px ∼ 25. We use a simulated N_px of 25 for the high-velocity comparisons in Figures 1 and 12 (corresponding to Δv_p ≈ 50 km s⁻¹ for this instrument). In such cases, one can obtain a constrained determination of the K_p of a planetary signal, although the finite instrument resolution is significant enough that prominent tails on signals in K_p space remain (see Figure 12). However, for sub-Neptune planets such as TOI-732 c, the change in radial velocity during the transit is typically much smaller. This is mainly driven by significantly smaller values for ${\phi }_{\max }$ due to the increased period of the planets. For this planet and instrument, accounting for the decrease in V_bary over the transit, the planet pixel shift is N_px ∼ 1.4 px, consisting of 1.3 px from the planet's orbit (Table 1) and 0.1 px from the barycentric velocity variation. We can therefore no longer constrain K_p and so any measurement of K_p is now redundant. We note that we are considering the number of wavelength channels crossed by the planet signal during transit, rather than the number of resolution elements. This is larger than the latter, typically by a factor of around 3, meaning it may be difficult to resolve the velocity shifts for planets with N_px ≲ 3, such as TOI-732 c.

The inability to constrain K_p could be an issue when aiming to confirm the robustness of a detection. Having only one useful parameter, V_sys, may not be sufficient to distinguish between spurious signals and a robust detection of a species in the atmosphere of the planet. For example, telluric and stellar features, and correlated noise from other sources, will likely produce similar unconstrained cross-correlation signals in K_p−V_sys space. To help us identify genuine planet cross-correlation signals from spurious noise, it would therefore be useful to include another parameter in our analysis.

As described in the previous section, it can be difficult to differentiate cross-correlation signals arising from low-velocity planets from those attributed to residual, quasi-static telluric and stellar signals. This is because we can no longer separate the planetary signal from these contaminants based on a well-constrained and large value for K_p characteristic of hot Jupiters. However, for planets whose change in radial velocity during the transit is small, it may still be possible to verify the origin of cross-correlation signals using a different metric. Whereas the planetary signal has a fixed transit duration, spurious signals in the CCF from, for example, residual telluric and stellar features may not exhibit the same phase dependence from being constrained to the transit.

To test this, we now treat ${\phi }_{\max }$ as a free parameter and explore over ${\phi }_{\max }$−V_sys space for constant K_p. All spectra are included in detrending, but for each point in ${\phi }_{\max }$ space the in-transit spectra are defined individually when summing the CCF in time. A S/N is then calculated for this set of in-transit spectra for each V_sys, eventually giving S/N as a function of ${\phi }_{\max }$−V_sys space, or N_in−V_sys space, where N_in is the number of in-transit spectra. This is done for constant K_p, here set to the known K_p of the planet in question (70 km s⁻¹; see Table 1). In Figure 4, the cross-correlation S/N across N_in−V_sys space is shown for each molecule considered. Unlike for K_p, the recovered signals for each of NH₃, CH₄, and H₂O are constrained at the N_in, or ${\phi }_{\max }$, of injection. This suggests that the signals arise only during the transit, and are therefore likely to be planetary in nature. As well as this metric helping to distinguish between real planetary signals and stellar/telluric cross-correlation features, it may also aid in providing an independent measurement of the transit duration of the planet, in cases where this may not be well known or constrained. For example, for TOI-732 c studied here, values for T₁₄ vary across the literature, and with significant errors (Cloutier et al. 2020; Nowak et al. 2020).

Zoom In Zoom Out Reset image size

However, a potential caveat of this method relates to the Rossiter–McLaughlin (R-M) effect (McLaughlin 1924; Rossiter 1924; Queloz et al. 2000). The modulation of stellar signals through the R-M effect may produce spurious peaks constrained in T₁₄ space, which will therefore mimic planetary signals under this metric. For molecules that exist in the M-dwarf photosphere as well as in the planet's atmosphere, such as H₂O, this metric could therefore be unreliable.

As discussed above, an important difference between the planet signal and the telluric/stellar residuals is that the planet signal is limited in time to just the in-transit spectra, whereas the others are present throughout our observations. In Figures 2 and 4, the planet signals recovered were injected into 21 of 33 spectra (in accordance with the transit duration listed in Table 1), with all 12 of the out-of-transit spectra used in the detrending. In the left-hand panel of Figure 5, we now trim our data set, in accordance with Section 3.1, to remove all of these out-of-transit spectra while the 21 in-transit spectra remain unchanged. In this manner, the injected NH₃ planet signal is also unchanged, with the only difference being that now the out-of-transit spectra are not available for the fitting of principal components in detrending. With no out-of-transit spectra, we are unable to recover the NH₃ planet signal, in direct contradiction to what is found in Figure 2. Likewise, in the right-hand panel of Figure 5, where the CH₄ planet signal is injected into all 33 available spectra prior to detrending, the recovered signal is also lost. Given these results, it follows that it is very likely that it is indeed the out-of-transit spectra which allow us to successfully use PCA-based detrending for low-velocity targets; without these spectra, the analysis is unable to recover the injected planetary signals. Although the planet signal is Doppler shifted over the transit by an amount less than an instrument resolution element, the out-of-transit spectra provide contrast against the telluric and stellar lines which prevent the planet lines being included in the principal components of the time-series spectra. Without out-of-transit spectra, the planet signal is included in the principal components of the spectra and can therefore not be isolated. Considering frequencies, such contrast or time inhomogeneity may have the effect of increasing the signal frequency of the planet signal relative to the telluric/stellar contaminants. This could mean that sufficient planet signal is able to survive the subtraction of principal components, or perhaps the application of a time-series high-pass filter, during detrending.

Zoom In Zoom Out Reset image size

We have therefore shown that it is indeed the contrast between the time-limited planetary signal and ubiquitous telluric and stellar features, provided by the out-of-transit spectra, that allows for the recovery of a low-velocity planet signal after detrending using PCA. In the high-velocity regime, the time-varying planet signal can still be recovered when the out-of-transit spectra are removed, albeit to a lower significance than in Figure 12 (e.g., a S/N of 5.6 for NH₃).

In Figure 1, it was shown that the planet signal degradation by the subtraction of principal components was faster for low-velocity planets than for warmer, short-period planets such as hot Jupiters. However, we here note that the degradation of the planet signal is not consistent across the transit. This differential signal loss across the transit is demonstrated in Figure 6. We find that the ingress and egress points of the transit correspond to higher-frequency components of the planet signal and are therefore degraded slower, and preserved for longer, after repeated subtraction of principal components from the spectra. On the other hand, the central part of the planet transit is included in the principal components of the spectra, and is therefore lost upon this subtraction. Therefore, after enough principal components are removed, the ingress and egress "edges" are the only surviving parts of the planet signal, and the available information about the planet is concentrated in these observations. A similar effect can be seen when applying a simple high-pass filter to a upside-down, transit-like, top-hat signal. This results in the loss of signal apart from at the edges, where the extra variation corresponds to a higher-frequency component of the signal.

Zoom In Zoom Out Reset image size

Our finding suggests that in the low-velocity regime it is particularly important to reprocess the model to reflect this differential erosion by PCA, with the in-transit signal being removed by application of PCA significantly more than the ingress/egress signal. This is especially true when considering any time- or phase-resolved spectroscopy of the planet, where constraints can be put on the variation of atmospheric conditions across different terminator viewing regions.

For the recovered signals shown in Figure 2, only the in-transit spectra were included when summing the CCF over time. To further demonstrate the differential erosion of the planet signal in phase shown here, we now repeat the recovery of these injections using different portions of the data set. We find that including only the first and last three spectra of the transit, i.e., those spectra just after/before the ingress/egress, in the S/N calculation returns a value S/N ≳ 3, whereas this threshold is not met when just the central 15 in-transit spectra are used, despite the increase in the number of in-transit spectra. We also calculate the S/N from the out-of-transit spectra only. In this case, we recover a strong anticorrelation signal (S/N = −4.7) at V_sys of 0 km s⁻¹, indicative of the presence of the out-of-transit, in-trail artifacts observed in Figure 6 and again suggesting that reprocessing is necessary (see Section 4.1.1). This unequal distribution of planet signal across the observations after detrending, characteristic of low-velocity planets, may be useful to consider when verifying real signals of such planets.

These effects seen here for low-radial-velocity-change planets may also mean that signals from such planets are potentially unsuitable for analysis with the Welch's t-test (Welch 1947), where a normal distribution in each of the in-trail samples is assumed. This metric would also likely be biased by the out-of-transit anticorrelated artifacts.

An alternative metric to cross-correlation by which the detection significance can be quantified is via a log-likelihood method (Brogi & Line 2019). To implement this, the spectra are cleaned, normalized, and detrended as before. Rather than be cross-correlated with a Doppler-shifted model spectrum, the residuals are then compared with a reprocessed model template using a log-likelihood framework. In this way, a grid of atmospheric models of different chemical abundances, temperature–pressure profiles, Doppler-shifts, and scale factors can be compared against the data for goodness of fit. The formulation of the log-likelihood log(L) for each trial model is described in the remainder of this section.

4.1.1. Model Reprocessing

4.1.2. Model Comparison

The resulting planet atmosphere model, reprocessed to account for the effects of detrending, can then be directly compared with the residual data to calculate a log-likelihood value at one point in K_p−V_sys−α space. While a cross-correlation to log-likelihood mapping can be used (Zucker 2003; Webb et al. 2022; Lafarga et al. 2023), we here calculate log(L) directly using the following equation (also see Brogi & Line 2019; Gibson et al. 2020):

$$\begin{eqnarray}&&\mathrm{log}(L)({K}_{{\rm{p}}},{V}_{\mathrm{sys}},\alpha )=-\displaystyle \frac{1}{2}\displaystyle \sum _{i,j}\left[{\left(\displaystyle \frac{{r}_{{ij}}-{m}_{{ij}}}{{\sigma }_{j}}\right)}^{2}+\mathrm{log}(2\pi {\sigma }_{j}^{2})\right],\end{eqnarray} \tag{ 3 }$$

where r_ij is the residual spectrum after detrending at wavelength channel j and observation i, m_ij is the Doppler-shifted and reprocessed model spectrum at the corresponding wavelength and observation, and σ_j is the error in the residuals for wavelength channel j. We have here assumed that for any given pixel the variance is contributed by the temporal variation of the flux throughout the observations. This is estimated for each wavelength channel using the median absolute deviation of the detrended residuals.

We repeat this comparison across model parameter space, with the resulting log-likelihood distribution used to calculate a Bayesian detection significance, as described in Section 4.2.

In addition to the log-likelihood distribution across parameter space for a model of a given chemical species, the log-likelihood value is calculated for a null model not containing that species. We here use a flat, zero-transit-depth spectrum (no planet transit: α = 0) as the null model. A Bayesian evidence, or marginal likelihood, is then calculated for each of the atmospheric and null models, from which a detection probability/significance can be determined. To calculate this Bayesian evidence, the log-likelihood distribution of the atmospheric model is integrated over the model parameters, modulated by their prior probability distributions:

$$\begin{eqnarray}&&{Z}_{m}={\int }_{\theta }\mathrm{log}(L)\pi (\theta )d\theta .\end{eqnarray} \tag{ 4 }$$

To obtain a Bayesian detection significance, the overall model evidence is first calculated as above. A Bayes factor, B, is then calculated as the ratio of the atmospheric model evidence to the null model evidence, Z_m/Z₀. This factor is converted to the detection significance, σ, of the atmospheric model relative to the null model using the following expression, from Welbanks & Madhusudhan (2021), based on work by Benneke & Seager (2013):

$$\begin{eqnarray}&&\sigma =\sqrt{2}.{\mathrm{erfc}}^{-1}\left[\exp \left({W}_{-1}\left(-\displaystyle \frac{1}{\mathrm{Be}}\right)\right)\right],\end{eqnarray} \tag{ 5 }$$

where erfc is the complementary error function and W₋₁ is the Lambert W function in its lower branch (k = −1).

The ability to integrate over parameters to form a marginal likelihood is useful in the low-velocity planet regime, where K_p is found to be unconstrained (Section 3.3). This may enable more consideration of the signal's degeneracy compared to cross-correlation, where in Section 3.4 we fixed K_p at its expected value when exploring N_in space.

We now apply this metric to our data set, as we did previously using cross-correlation, to recover an injected H₂O signal. We use this model planet signal as a conservative test case since it was recovered weakly with S/N = 2.1 using the cross-correlation method in Figure 2. We again detrend the spectra by subtracting seven principal components. To compare the spectral residuals to the model, we calculate log-likelihood values across the three-dimensional model space spanned by K_p, V_sys, and α.

The recovered signal is shown in Figure 7. In the top panel, the likelihood function is marginalized over just α (using uniform priors of 0.5–1.5) to give the model evidence and then Bayesian detection significance as a function of K_p and V_sys. In the bottom panel, we show the posterior probability distribution across V_sys−α space calculated from our likelihood distribution using emcee (Foreman-Mackey et al. 2013). In doing so, we use the above prior for α and uniform priors of 40–100 km s⁻¹ for K_p and 15–25 km s⁻¹ for V_sys. Values of V_sys = 20 ± 2 km s⁻¹ and α = 0.8 ± 0.2 are found. Clear preference is shown for an atmospheric model with parameters consistent with the injection. Marginalizing all the model parameters over the entire uniform prior volume, as described in Section 4.2, results in a detection significance of 2.9σ for H₂O compared with the null model. As in the cross-correlation case, no significant signals are recovered without an injection, demonstrating the successful PCA-based extraction of the injected planet signal from spectra of a low-velocity, temperate planet. Our results may therefore suggest that the Bayesian framework used here is more sensitive to weak planet signals than the cross-correlation metric used in Section 3, for example due to model reprocessing and the extra sensitivity to line depths and shapes. Once again, removing all 12 out-of-transit spectra prior to detrending means the injection is no longer recovered.

Zoom In Zoom Out Reset image size

PCA-based detrending uses changes/inhomogeneities in the time-series spectra, such as the change in the planetary radial velocity during the transit or the existence of transit edges, to differentiate the Doppler-shifting/temporary planetary lines from the quasi-static and persistent telluric and stellar lines. We here investigate how the quality of detrending, reflected by the recovered S/N of an injected signal, is dependent on the number of instrument wavelength channels by which the planetary signal is Doppler shifted during the transit, N_px, as introduced in Section 3.5. This is dependent on K_p, the maximum phase during transit ${\phi }_{\max }$, and the instrument pixel velocity V_px, as defined in Equation (2).

We inject and recover CH₄ models over a range of simulated K_p and ${\phi }_{\max }$, for constant V_sys and numbers of in-transit and out-of-transit spectra (−15 km s⁻¹, 16, and 12, respectively). The S/N of the recovered signal is plotted across this K_p−${\phi }_{\max }$ space in Figure 8, with N_px contours, calculated from these orbital parameters and the resolution of the IGRINS data used here, shown. This figure illustrates the strong dependence of the detrending on the planetary pixel shift, and shows how it can be useful to consider the compound parameter of N_px alone, rather than each of K_p, ${\phi }_{\max }$, and V_px individually, when predicting whether a given detection can be made. As expected, in the low-orbital-velocity regime, the recovery of a signal increases with increasing N_px. As we approach greater values of N_px (warmer planets with larger changes in radial velocity during transit), the increase in detection S/N begins to slow, after which more time variation in the planet signal leads to diminishing returns. We once again find that planetary signals can still be recovered when the change in the planet's radial velocity is small, including when it is smaller than that corresponding to an instrument resolution element, i.e., N_px ≲ 3.

Zoom In Zoom Out Reset image size

We have shown that it is possible to probe the atmospheres of slower-moving planets using PCA-based detrending due to the contrast provided by the out-of-transit spectra. We now investigate further how the number of out-of-transit spectra affects our ability to isolate a planet signal from the telluric/stellar contaminants using PCA-based detrending. As before, we do this by measuring the S/N at which various injected signals are recovered.

For a fixed number of in-transit spectra, we vary the number of out-of-transit spectra included in the detrending. Given that in Section 5.1 we demonstrate the importance of the planetary pixel shift, Figure 9 shows the regions of the space spanned by N_px and the ratio of out-of-transit spectra to in-transit spectra which are suitable for obtaining a detection in high-resolution transmission spectroscopy. A CH₄ model injected into 16 in-transit spectra is used here, and the S/N values presented in this figure are evaluated by averaging over V_sys space. We demonstrate that the recovery of a given signal increases when a greater fraction of out-of-transit spectra are included in the detrending. This is expected as the existence of the planet signal in a smaller fraction of the data creates more contrast with the tellurics and so means less of the signal is included in the principal components of the time-series spectra. The effect is to promote a more effective separation of telluric and planetary signals; even for very small values of N_px, atmospheric characterization can be achieved if sufficient out-of-transit spectra are included in detrending. The planetary pixel shift of TOI-732 c is marked, calculated using the values for K_p and ${\phi }_{\max }$ in Table 1, and assuming an instrument resolution equal to that of IGRINS (V_px ≈ 2 km s⁻¹). In the case of these simulations, we find that, averaged over V_sys, a TOI-732 c–like sub-Neptune planet signal can be recovered with S/N > 3 when the ratio of out-of-transit spectra to in-transit spectra is ≳0.8. This trend could therefore be an important consideration when planning observations of low-velocity planets.

Zoom In Zoom Out Reset image size

We now explore the erosion of a low-velocity NH₃ planet signal with each principal component subtracted from the spectra, as in Figure 1, for different out-of-transit to in-transit ratios (Figure 10). The number of in-transit spectra is again set as 16 in this test, with a higher-velocity signal also shown for comparison. As expected from Figure 9, the planet signal is lost more quickly when there are fewer out-of-transit spectra. In the case of no out-of-transit spectra and a subresolution change in radial velocity, the remaining planet signal is negligible after just one iteration. When more out-of-transit spectra are included in detrending, or even when there are no out-of-transit spectra but the change in the planet's radial velocity is larger, more of the planet signal evades inclusion in the principal components of the spectra. In these cases, the signal is still able to be recovered after a number of principal components have been subtracted to remove sufficient correlated noise.

Zoom In Zoom Out Reset image size

Sometimes only partial transits are available; conditions, for example air mass, may not permit the observation of a full transit plus both pretransit and posttransit spectra. In this case, the number of pretransit and posttransit observations will likely be unequal. There may be out-of-transit spectra available on only one side of the transit, or perhaps only a fraction of the transit itself was observed. In this section, we explore the case where there is an unequal number of pretransit and posttransit observations. Until now, we have injected signals into in-transit spectra and selected the out-of-transit spectra symmetrically in phase, as detailed in Section 3.1. Since our phases are symmetrical about zero, this approach has led to roughly equal numbers of out-of-transit spectra on each side of the transit.

We complete injection and recovery tests as in Figure 2, this time simulating the case where all out-of-transit spectra are either before or after the 21 in-transit observations (Figure 11). We find that while we are still able to recover signals from partial transits, it can be more difficult, demonstrated by the reduced S/N values found here for NH₃ (3.1 and 3.9) compared to the value from Figure 2 (5.0). Note the greater S/N in the second partial transit example (bottom panel, with all the out-of-transit observations after the transit), likely due to the higher air mass toward the end of our observations. In the case of a partial transit, the lack of out-of-transit spectra on one side of the spectrum removes an "edge" from the transit profile, further reducing the frequency of the planet signal and making it more difficult to separate from the telluric, stellar, and other sources of correlated noise (see Section 3.6). It therefore is preferable to have both sufficient pre- and posttransit observations.

Zoom In Zoom Out Reset image size

We note that our work here has focused solely on transmission spectroscopy. In emission spectroscopy, the typical existence of the planet signal in every observation means it is removed by PCA during detrending if the radial velocity change is too small. Despite this, our findings here suggest that it may be technically possible to detrend emission spectra in this regime if observations of the system include some of the secondary eclipse. These in-eclipse spectra, taken before/after/between observations in which the planet is visible, could then perhaps act as the contrasting "out-of-transit" spectra required for the detrending to work. In this case, the setup reduces to the case of partial transits, with "out-of-transit" observations taken only before/after the transit. While we note that temperate planets will not typically have detectable emission features, and therefore their lack of radial velocity change during observations is unlikely to be the main obstacle, this idea may find relevance in the study of young giant planets.