Primary Velocity and Orbital Phase Effects on Planetary Detectability from Small Epoch Number Data Sets

Cross correlation analyses of high resolution spectroscopic data have recently shown great success in directly detecting planetary signals and enabling the characterization of their atmospheres. One such technique aims to observe a system at multiple epochs and combine the measured planetary radial velocities from each epoch into a measurement of the planetary Keplerian orbital velocity $K_p$, constituting a direct detection of the planetary signal. Recent work has shown that in few-epoch ($\sim$5) data sets, unintended structure can arise at a high level, obscuring the planetary detection. In this work, we look to simulations to examine whether there are ways to reduce this structured noise in few-epoch data sets by careful planning of observations. The choice of observation date allows observers to select the primary (stellar) velocity - through a set systemic velocity and chosen barycentric velocity - and the planetary orbital phase, and so we focus on the effects of these two parameters. We find that epochs taken when the primary velocity is near zero, and the stellar lines remain relatively fixed to the telluric rest-frame, greatly reduce the level of structured noise and allow for much stronger planetary detections, on average more than twice the significance of detections made with epochs using randomly selected primary velocities. Following these results, we recommend that observers looking to build up high-resolution multi-epoch data sets target nights when their system has a near-zero primary velocity.


INTRODUCTION
As thousands of exoplanets are being discovered through indirect methods such as transit and radial velocity surveys, astronomers have begun to consider how to follow-up on these detections and measure the planets' atmospheric properties, especially the presence and relative abundances of molecular species, the atmospheric pressure/temperature profiles, and the nature of winds and planetary rotation. High resolution cross correlation spectroscopy has been introduced as an effective technique to directly detect planets' thermal emission and begin to characterize their atmospheres (e.g., Brogi et al. 2012;Birkby et al. 2013;Lockwood et al. 2014). High resolution cross correlation spectroscopy works by allowing researchers to disentangle planetary and stellar radial velocities. By fitting the planetary radial velocities with an equation for the orbital motion, observers can constrain the amplitude of that motion, called the planetary Keplerian orbital velocity, K p . With prior knowledge of the stellar mass and the stellar Keplerian orbital velocity K * from optical radial velocity measurements, we can constrain the true mass and orbital inclination of the planet. Further, planetary models with different assumptions about various atmospheric properties can be cross correlated against the data and the resulting strength of the planetary detection can be used to understand the true nature of the planetary atmosphere.
Two approaches have been applied to constrain K p from the planetary radial velocities. In one, observers target a system at times when the line-of-sight planetary acceleration is largest, e.g. near inferior/superior conjunction (e.g., Snellen et al. 2010;Brogi et al. 2012Brogi et al. , 2013Birkby et al. 2017;Guilluy et al. 2019). Then, by observing for many (∼5-7) hours in one stretch, they can watch the planetary signal shift systematically with respect to the fixed stellar and telluric reference frames, as its lineof-sight velocity changes. The Keplerian velocity K p is measured through a fit to this changing planetary radial velocity. While this technique has been very effective for hot Jupiters, it requires a significant change in the lineof-sight orbital velocity (tens of km/s) over the course of a single continuous observing sequence. This will preclude its application to longer period planets, including those is the nearest M star habitable zones. Another approach, first introduced by Lockwood et al. (2014), limits wall-clock observing times to ∼2-3 hours to obtain measurements that do not allow the planetary signal to cross detector pixels. After building up a data set of several such measurements around the planet's orbit, the data can be fit to constrain K p . We will call this approach the "multi-epoch" approach. Buzard et al. (2020) used the multi-epoch approach to detect the thermal emission from the hot Jupiter, HD187123b. In that work, we introduced a simulation frame that could account for a portion of the structured noise that arose in the K p detection space. With these simulations, we showed how beneficial many-epoch (∼20) data sets are over few-epoch (∼5) data sets, even if they share the same total S/N. This was because in the fewepoch simulations structured noise resulting from unwanted correlation between the planetary spectral template and the observed stellar signal was well in excess of the shot noise readily obtainable with Keck on bright stars. Large number epoch data sets afford more variation between the stellar and planetary signals that works to beat down this source of structured noise. However, it can be difficult to build up to such large data sets with highly over-subscribed telescopes such as Keck, especially with the current generation of high resolution echelle spectrographs that have modest instantaneous spectral grasp in the thermal infrared, such that significant integration times are needed per epoch. In this work, we investigate whether there is a way to provide more efficient detections with few-epoch data sets by carefully selecting which nights we choose to observe. To do so, we consider the effects of the primary (stellar) velocities and orbital phases at each epoch.
At any given observation time, the stellar velocity in a system will be determined by the systemic velocity, the barycentric velocity from the component of the Earth's orbital motion in the direction of the system, and the radial velocity caused by the planetary tug on the star. The current NIRSPEC does not have the velocity precision necessary to resolve the radial velocity caused by a planet; while NIRSPEC L band velocity precision is ∼ 3.2 km/s, the stellar RVs caused by even the massive hot Jupiters are 0.1-0.2 km/s. Assuming a constant systemic velocity, then, the barycentric velocity is variable and can be chosen by when the observing night is scheduled. We will consider the primary velocity at a given observation time t obs , v pri (t obs ) = v sys − v bary (t obs ). (1) The planetary velocity is similarly comprised of a dynamical radial velocity (which, unlike the stellar RV, is large enough to be resolved by NIRSPEC), the systemic velocity, and the barycentric velocity. The magnitude of the planetary gravitational radial velocity signature is K p , and depends on the planetary and stellar masses, and the orbital inclination, semi-major axis, and eccentricity. The orbital phase, M , of the planet will determine the magnitude of the planetary radial velocity given K p , and will therefore determine the planetary velocity relative to the stellar velocity. If the radial velocity parameters have been determined through optical, stellar RV measurements, M can be calculated at any time through the equation, where T 0 is the time of inferior conjunction and P is the orbital period. As a function of t obs , M can also be chosen with careful observation scheduling. With these two parameters, and assuming a circular orbit, the secondary, or planetary, radial velocity can be described as v sec (t obs ) = K p sin(2πM (t obs )) + v pri (t obs ).
Since the primary velocity (v pri ) and orbital phase (M ) can both be selected by the choice of observing nights, we set out to understand how different combinations of primary velocity and orbital phase epochs affect the detectability of planetary Keplerian orbital velocities, K p , for a modest, and readily obtainable, 5 epoch data set.
The rest of this work is organized as follows. In Section 2, we describe the planetary spectral models used in these simulations, how the simulations are generated, and how they are analyzed. In Section 3, we consider the effectiveness of different groupings of orbital phases and primary velocities. We examine whether the magnitude of K p affects the results of these primary velocity/orbital phase simulations in Section 4. In Section 5, we attempt to see whether a combination of NIRSPEC data epochs agree with these simulation results. In Section 6, we consider the primary velocity and orbital phase effects on larger data sets. Finally, we discuss some implications of these results in Section 8 and conclude in Section 9.

Spectral Models Used
For these simulations, we used a spectral model generated from the PHOENIX stellar spectral modeling framework (Husser et al. 2013). We interpolated the effective temperature, metallicity, and surface gravity to those of the sun-like star HD187123 (T eff = 5815 K, [Fe/H] = 0.121, and log(g) = 4.359; Valenti & Fischer (2005)).
Our planetary thermal emission model was generated from the SCARLET framework (Benneke & Seager 2012Benneke 2015;Benneke et al. 2019a,b) at R = 250, 000. This framework computes equilibrium atmospheric chemistry and temperature structure assuming a cloud-free atmosphere with a solar elemental composition, efficient heat redistribution, and an internal heat flux of 75 K. We assume a solar metallicity and C/O ratio. The SCARLET model framework includes molecular opacities of H 2 O, CH 4 , HCN, CO, CO 2 , NH 3 , and TiO from the ExoMol database (Tennyson & Yurchenko 2012;Tennyson et al. 2020), molecular opacities of O 2 , O 3 , OH, C 2 H 2 , C 2 H 4 , C 2 H 6 , H 2 O 2 , and HO 2 (HITRAN database by Rothman et al. (2009)), alkali metal absorptions (VALD database by Piskunov et al. (1995)), H 2 broadening (Burrows & Volobuyev 2003), and collisioninduced broadening from H 2 /H 2 and H 2 /He collisions (Borysow 2002). The atmosphere does not have an inverted thermal structure in regions close to the molecular photosphere.

Generation of Simulated Data
In this work, we generated simulated data sets following the framework introduced by Buzard et al. (2020). In short, the stellar and planetary models are scaled by assumed stellar and planetary radii squared and shifted to velocities determined from Equations 1 and 3. Next, the stellar spectrum is interpolated onto the planetary wavelength axis, and the two models are added. The stellar continuum is removed with a third order polynomial fit to the combined spectrum from 2.8 to 4 µm in wavenumber space. The spectrum is broadened with a Gaussian kernal fit to real NIRSPEC data. Finally, the spectrum is interpolated onto a NIRSPEC data wavelength axis, saturated telluric pixels from the data (where tellurics absorb more than about 40% of the flux) are masked, and Gaussian noise is added. The masking of saturated telluric removes about 40% of the data. Non-saturated tellurics are assumed to be perfectly corrected.
For these simulations, we assume a 1 R Jup planet and a 1 R star. Unless otherwise stated, these simulations approximate post-upgrade L band data. The upgraded NIRSPEC instrument was first available in early 2019 (Martin et al. 2018). Across the L band, it doubled the number of pixels per order (from 1024 to 2048), increased the number of usable orders on the chip (from 4 to 6), and nearly doubled the spectral resolution (from ∼25,000 to ∼40,000). The Gaussian kernals used to broaden the simulated data and wavelength axes with their corresponding locations of saturated telluric pixels were taken from the 2019 Apr 3 and 2019 Apr 8 NIRSPEC data of HD187123 presented in Buzard et al. (2020). Each epoch has six orders, covering wavelengths of approximately 2. 9331-2.9887, 3.0496-3.1076, 3.1758-3.2364, 3.3132-3.3765, 3.4631-3.5292, and 3.6349-3.6962 µm. The average instrumental resolution is about 41,000. We applied the average S/N per pixel from the 2019 Apr 3 and 8 data of 2860 to each epoch, which resulted in a total S/N per pixel of 6390 for the 5 epoch simulations.

Analysis of Simulated Data
The simulated data sets are analyzed analogously to the data presented in past multi-epoch detection works, e.g., Piskorz et al. (2018), Buzard et al. (2020). A twodimensional cross correlation, TODCOR as described in Zucker & Mazeh (1994), is used to measure the stellar and planetary velocities at each epoch. Segments of data (e.g., orders and pieces of orders after saturated telluric pixels are masked) are cross correlated separately and converted to log likelihood functions in order to be combined. In this work, we use the Zucker (2003) log(L) approach to convert cross correlations to log likelihoods; the "Zucker log(L)" approach is described and differentiated from the "Zucker ML" approach in Buzard et al. (2020). This approach converts cross correlations to log likelihoods as where R is the two-dimensional cross correlation and n is the number of pixels in the data segment. Once the two-dimensional cross correlation of each epoch is converted to a two-dimensional (stellar and planetary velocity shifts) log likelihood surface, the planetary log likelihood cut is taken from the measured stellar velocity. The measured stellar velocity is always consistent with the expected stellar velocity from Equation 1. The planetary log likelihood curves at each epoch are converted from v sec to K p space by Equation 3, and added. This planetary log likelihood vs. K p curve is calculated from −150 ≤ K p ≤ 150 km/s.

PRIMARY VELOCITY SIMULATIONS
To study how primary velocities and orbital phases affect planetary detectability, we generate sets of simulated data with different combinations of primary velocities and orbital phases. These data sets all consider five epochs and have K p set at 75 km/s. For these simulated data sets, we allow v pri to range from -30 to 30 km/s, the rough maximum variation given by the Earth's orbital velocity, v bary . We create five different groupings of primary velocities: (1) a most blue-shifted v pri sample, in which the primary velocities at all five epochs are pulled from a uniform distribution from -30 to -28 km/s; (2) an even v pri sample, in which the five primary velocities are evenly spaced from -30 km/s to 30 km/s; (3) a most red-shifted v pri sample, in which the five epochs are pulled from a uniform distribution from 28 to 30 km/s; (4) a near-zero v pri sample in which the five primary velocities are pulled from a uniform distribution from -2 to 2 km/s; and (5) a random v pri sample. For the evenly spaced v pri sample, the five epochs are pulled from uniform distributions covering: -30 to -29 km/s, -16 to -14 km/s, -1 to 1 km/s, 14 to 16 km/s, and 29 to 30 km/s. These slight variations in the most blueshifted, most red-shifted, even, and near-zero primary velocity groups better resemble actual observations that could be scheduled than if all five large v pri epochs had exactly 30 km/s, for example. For the randomly sampled primary velocity group, we choose the Earth's orbital position from the uniform distribution from 0 to 2π. These positions are then converted to barycentric velocities assuming the maximum barycentric velocity is 30 km/s. This results in a bimodal barycentric velocity distribution that is relatively uniform through the central velocities and increases significantly towards ±30 km/s. We consider a systemic velocity of 0 km/s, so the resulting random primary velocity distribution has the same shape as the barycentric velocity distribution (v pri = −v bary ).
If the systemic velocity were non-zero, the primary velocity distribution would be shifted and the probability would increase towards its maximum (v sys − min(v bary )) and minimum (v sys − max(v bary )) values. We discuss how this may affect the results of the random primary velocity simulations in Section 8. Realistically, systems are not observable from the Earth for the full year. A pull from half of the Earth's orbit (e.g., v bary = 30 cos(0) to 30 cos(π)) results in the same random primary velocity probability distribution, so we use this moving forward. We contemplate further effects of target accessibility in Section 8.
We split combinations of orbital phases, M , up into three groups: (1) all five epochs near conjunction, (2) all five epochs near quadrature, and (3) five epochs evenly spaced around the orbit. The five near-conjunction epochs are pulled randomly from the uniform distributions, 0 ± 0.02 (inferior conjunction) and 0.5 ± 0.02 (superior conjunction), and the quadrature epochs are pulled from the uniform distributions, 0.25 ± 0.02 and 0.75 ± 0.02. The evenly spaced M epochs have one epoch pulled from similarly wide uniform distributions centered on each of 0.05, 0.25, 0.45, 0.65, and 0.85. This is only one example of an evenly distributed set of orbital phases, and we expand the analysis to include other combinations of orbital phases in Section 3.1. Figure 1 shows the results of these combinations of primary velocities and orbital phases, with the five primary velocity groups taking up different subplots, and the three orbital phase groups shown in different colors. For each primary velocity subplot, near quadrature orbital phases are shown in dark blue, evenly spaced orbital phases in light purple, and near conjunction phases in green.
Several notable results stand out. First, there is very little structure in any of the simulations that considered all five epochs near conjunction. These simulations show no convincing detections of K p . This trend makes sense because when the planet is near either inferior or superior  conjunction, it will have little to no line-of-sight velocity difference from its star. Regardless of what K p is, the planet line-of-sight velocity at conjunction will simply be equal to the primary velocity. Conjunction epochs on their own are not useful for constraining the Keplerian orbital velocity through the technique that aims to measure stationary planetary velocities at multiple epochs. This is in contrast to cross correlation techniques that aim to measure changing planetary velocities (e.g., Snellen et al. 2010;Brogi et al. 2012); they actually prefer conjunction epochs, during which the planetary acceleration is the largest. It is also useful to note that these techniques that target changing planetary velocities would also require higher spectral resolution than the technique that targets stationary planetary velocities. The simulations with evenly spaced orbital phases and orbital phases near quadrature do have a peak at K p in all of the primary velocity groups, but there is often large off-peak structure at the same magnitude, if not larger, than the true peak. Both Buzard et al. (2020) and Finnerty et al. (2021) found that at the S/N used in these simulations (> 2500 per pixel per epoch), shot noise has very little effect on the log likelihood surface, and the off-peak structure is due instead to non-random, structured noise. This structured noise is caused by correlation between the planetary spectral model template and the stellar features in the simulated data. Finnerty et al.
(2021) removed this structure in their simulations by subtracting a stellar-only log likelihood curve, which comes from simulated data generated with no planetary signal and then cross correlated in the same two-dimensional way with a stellar and a planetary spectral model. This approach can nearly eliminate all of the off-peak structure in simulations, but it would not be as effective on data due to a variety of factors including mismatches between the real and model stellar and planetary spectra and imperfect removal of tellurics from the data. We therefore do not try to remove this off-peak structure. Instead, we look for combinations of primary velocities and orbital phases that can reduce it by design.
Interestingly, we see in Figure 1 that the simulations with primary velocities around 0 seem to show stronger peaks at K p relative to the noise than for the other primary velocity groups. This appears to be true for both the evenly spaced orbital phases and the near quadrature orbital phases. While these simulations consider a planetary spectral model without a thermal inversion, simulations generated and analyzed with a planetary model with an inversion showed a similar trend in that nearzero primary velocity epochs produced the strongest detections.
In order to investigate whether this trend that epochs taken when the primary velocity of the system is near zero provide stronger detections is more broadly true, we generated five-epoch data sets within each of the five primary velocity sampling groups, but with orbital phases pulled from a uniform distribution from 0 to 1, i.e. the full orbit. These data sets are likely more representative of real data sets that could be obtained from systems of interest too. While the barycentric velocity changes over the course of an Earth year, the orbital phase changes on the time frame of the planet's year. For hot Jupiters, the orbital phase changes significantly from night to night, making it difficult to obtain multiple epochs with the same orbital phase (especially when trying to schedule nights which will provide useful epochs for multiple targets). The primary velocity, on the other hand, will be approximately the same on a monthly timescale. So, simulations with set primary velocities and randomly picked orbital phases might be a good approximation to data sets that could be easily obtained.
We generate 100 sets of five epoch simulations with randomly chosen orbital phases for each of the five primary velocity groups. We define a parameter, γ, to quantify each combination of orbital phases, as follows.
where N is the number of epochs. For epochs at quadrature, sin(2πM ) = ±1, and for epochs at conjunction, sin(2πM ) = 0. Defined this way, γ = 0 if all of the epochs are at conjunction and γ = 1 if all of the epochs are at quadrature. After analyzing each of the 500 simulations, we attempt to fit Gaussians to the resulting normalized log likelihood curves. We first normalize the curves by subtracting the mean of the curve from -150 to 0 km/s. We choose to normalize by the mean of this range because the way that we define K p enforces that it must be a positive value, meaning that the log likelihood curve at negative values of K p must be completely due to structured noise, and not to any real planetary signal. Because the simulated data is generated with a K p of 75 km/s, we fit Gaussians with an initial mean of 75 km/s, σ of 10 km/s, and height equal to the normalized log likelihood value where K p = 75 km/s. The results of these Gaussian fits are shown in Figure 2. In the top row, we plot the height of each Gaussian fit over the standard deviation of the curve from -150 to 0 km/s. We use this as a means to show how much stronger the true planetary peak is than the structured noise. In the bottom row are the standard deviations of the Gaussian fits. We plot a blue point for each simulation for which a Gaussian could be fit within 1σ of the true value of K p , 75 km/s, and a red point at 0 for both of the Gaussian parameters when it could not. The stars are the Gaussian fits to the simulations shown in Figure 1.
These results confirm that primary velocity near 0 km/s will generally allow for stronger detections of the planetary signal and more confident measurements of K p . The planetary peak was detected in all 100 of the nearzero primary velocity cases, but only 91 of the even primary velocity cases, 64 of the most red-shifted primary velocity cases, 90 of the most blue-shifted primary velocity cases, and 86 of the random primary velocity cases.
Further, the heights of the Gaussians fits to the nearzero primary velocity cases relative to the noise are much larger, on average, than for any of the other primary velocity cases.
Interestingly, we see no obvious relationship between the γ for a simulation and its peak height over the noise for any of the primary velocity groups. We suspect that at the larger values of γ, near quadrature, the planetary peak becomes resolved, leading to a larger height, but the noise structure also becomes narrower and of larger amplitude, so the increase in the peak height and noise level balance each other out. With the lack of dependence on γ, we can consider the mean and standard deviation of the peak heights over noise. Considering only those simulations in which the planetary peak was detectable, the peak height over noise was 6.2 ± 1.9 for the nearzero primary velocity case, while it was 2.6 ± 1.0 for the even primary velocity case, 2.3 ± 1.0 for the most redshifted primary velocity case, 1.9 ± 0.7 for the most blueshifted primary velocity case, and 2.7±1.0 for the random primary velocity case. The simulation results show a significant amount of scatter around these averages, even at a single value of γ, as seen in Figure 2. We found no significant relationship between Gaussian height over noise and the mean orbital phase, standard deviation of the orbital phases, or the standard deviation of the | sin(2πM i )| values, though. We also ran 100 simulations with the same orbital phases (γ = 0.67) and primary velocities (near-zero) to see how much of the scatter could be explained by white noise. The Gaussian heights from these simulations had a standard deviation of only 0.3, less than the scatter in any of the five primary velocity groups. This implies that with S/N of 2860 per pixel per epoch, the structured noise related to the combination of orbital phases dominates over random Gaussian noise. This is consistent with findings from Buzard et al. (2020) and Finnerty et al. (2021).
The lower panels of Figure 2 show the widths of the Gaussian fits. The gray dashed line in each subplot shows the average velocity precision of the upgraded NIRSPEC data on which these simulations were based, at 3.1 km/s. The widths of the near-zero primary velocity simulations show an interesting trend. At large values of γ, the widths are quite small and do not have much variation. The widths increase toward intermediate γ values in both magnitude and degree of variation. The representative near-conjunction simulation, shown by the blue star near γ = 0, has by far the largest width. The trend in width magnitude reflects the fact the conjunction epochs have very little constraining power on K p while quadrature epochs are the most effective for constraining K p . The degree of variation in Gaussian width at intermediate values of γ as opposed to large or small values can also be explained. While there is only one combination of epochs each that will give a γ value of 0 (all at conjunction) or 5 (all at quadrature), there are many different combinations of epochs that could result in an intermediate value of γ. For instance, the γ values of 1000 sets of 5 epochs with orbital phases pulled from a uniform distribution from 0 to 1 form an approximately Gaussian shape with a mean of 0.64 and a standard deviation of 0.14. With more cases at intermediate values of γ, there will be more variation as some of them will provide better constraints on K p than others. The four primary velocity groups other than the nearzero group do not show the same strong relationship between the Gaussian widths and γ. The most blue-shifted and most red-shifted primary velocity groups do show some evidence of a corner where the widths increase below a certain γ value but this behavior is not nearly as strong as the trend in the near-zero primary velocity case. We suspect that the higher planetary peaks in the near-zero group are better fit by a Gaussian, meaning that the widths are more representative of the true planetary detection peak than for the other primary velocity groups. This can be corroborated by the mean R 2 value of the detected planetary peaks in each v pri group, which compares the goodness of the Gaussian fit to that of a horizontal line at the mean of the simulated log likelihood curve. The mean R 2 values of the near-zero, most blue-shifted, even, most red-shifted, and random primary velocity groups are 0.69, 0.25, 0.34, 0.24, and 0.36, respectively. The especially low R 2 values of the four primary velocity groups other than the near-zero one reflect the high levels of structured noise. They also support our conjecture that the widths of those four groups show less dependence on γ than the near-zero primary velocity group because the Gaussian fits are not accounting for the planetary peak structure as accurately.
Since neither the most blue-shifted (most negative) nor the most red-shifted (most positive) primary velocity groups were able to strongly detect the planetary signal, going forward, we will refer to both as the "largest absolute primary velocity group." Doing so allows us to focus on the magnitude of velocity separation between the stellar signal and the telluric frame, rather than the direction in which the stellar signal has moved.
Collectively, these results show that a set of 5 epochs with randomly selected orbital phases will have the best chance of showing a strong detection of the planet if they are taken during times when the system's velocity is canceled out by the Earth's velocity in the direction of the system. The closer these epochs are to quadrature, the better the data will be able to constrain the value of K p . Further, we would expect that obtaining data from both quadrature positions (M = 0.25 and 0.75) would be better for constraining K p than data at just one quadrature position because having data at both quadrature positions would give us access to different velocity shifts relative to the telluric frame and so collectively more complete wavelength coverage of the planetary spectrum. Finnerty et al. (2021) showed that a larger spectral grasp can drastically increase detection significance; doubling the grasp increased the significance by nearly a factor of 2. These predictions should be useful for the planning of future multi-epoch observations.

MAGNITUDE OF K P
All simulations presented in Section 3 considered data sets generated with a K p of 75 km/s. We showed that near-zero primary velocity epochs allow for the strongest planetary detections. However, we might expect to encounter a challenge, especially as K p decreases, with setting the primary velocity to 0, or in other words, allowing very little velocity shifting between the stellar and telluric spectra. If the stellar spectrum is not velocity shifted relative to the telluric spectrum, as K p decreases, the planetary lines will not be able to stray much from the telluric spectrum either.
The value of K p is set by both the semi-major axis and the orbital inclination. A decrease in K p due to a larger semi-major axis would be accompanied by a colder planetary effective temperature, while a decrease in K p solely due to a smaller inclination would not affect the planetary temperature. While either decrease in K p would limit the separation between planetary and telluric lines in near-zero primary velocity epochs, cooler planetary atmospheres could further complicate the issue. While hot Jupiter L band spectra are dominated by water features, their water is much hotter ( 1000 K) than water in the telluric spectrum (∼ 300 K), resulting in very different spectral line shapes, position, and relative contrast. However, as the planetary effective temperature decreases, its spectrum will more and more resemble that of the Earth's. Then, with neither very different temperatures altering the shape of the planetary spectrum from the telluric spectral shape or much velocity shifting off of the telluric spectrum, primary velocity near-zero epochs may no longer be as useful. Such cool planets will require large epoch number data sets to be detected.
To test how varying K p affects our simulation results, we generate 100 simulations with near-zero primary velocities and randomly selected orbital phases with K p values of 37.5, 75, and 150 km/s. We maintain a common planetary effective temperature in order to examine the effects on hot Jupiter detectability with a changing orbital inclination, not on the detectability of planets with different effective temperatures due to different semi-major axes. In analyzing these simulated data sets, we calculated the planetary log likelihoods vs. K p from −250 ≤ K p ≤ 250 km/s to allow for sufficient parameter space to robustly constrain the 150 km/s detections. Figure 3 shows the Gaussian heights relative to the noise and the Gaussian widths of these simulations. Of all 300 simulations, only one of the K p = 37.5 km/s cases was unable to fit the peak within 1σ.
We used a Kolmogorov-Smirnov (KS) test to determine whether there was any statistical difference between the Gaussian heights over noise and hypothesis that the 100 γ values of each of the K p cases were pulled from the same distribution. We know the null hypothesis to be true in this case; all 300 γ values were pulled from the same distribution, the conversion of the 5 M values uniformly pulled from 0 to 1 through Equation 5. This result then helps to validate the use of the KS test.
The KS p-values between the Gaussian heights over noise were 0.0030 (37.5 vs. 75 km/s), 0.031 (37.5 vs. 150 km/s), and 0.89 (75 vs. 150 km/s). The Gaussian width p-values were 0.031 (37.5 vs. 75 km/s), 0.069 (37.5 vs. 150 km/s), and 0.56 (75 vs. 150 km/s). The peak heights and widths of the 75 and 150 km/s cases can both be assumed to be pulled from the same parent distributions. On the other hand, KS tests reject the hypotheses that the 37.5 km/s heights are pulled from the same parent distribution as the 75 km/s heights at the 3.0σ level and from the same parent distribution as the 150 km/s heights at the 2.2σ level. They reject a common parent distribution between the 37.5 and 75 km/s Gaussian widths at the 2.1σ level and between the 37.5 and 150 km/s widths at the 1.8σ level. While these levels of statistical rejection of the null hypotheses are mostly in the "weak" to "moderate" support categories (e.g., Gordon & Trotta 2007), when considered alongside the means and standard deviations of the Gaussian heights of each distribution, they do start to show weaker planetary detectability at lower values of K p . The 37.5, 75, and 150 km/s sets of simulations have average Gaussian peak heights over the noise level of 5.5±2.4, 6.2±1.9, and 6.2 ± 2.3, respectively. While the 37.5 km/s peak heights were moderately lower than the 75 km/s and 150 km/s peak heights here, they are still on average significantly higher than the 75 km/s peak heights measured from the primary velocity groups not near zero, as described in Section 3.
For hot planets, we find that near-zero primary velocity epochs allow for stronger detections than other combinations of primary velocities even as K p gets quite small. This result may be challenged as we look to cooler planetary atmospheres which will have a higher degree of spectral similarity to our own telluric spectrum.

COMPARISON TO DATA
We were interested to test whether our prediction that epochs with primary velocities near zero would give stronger detections than other samples of primary velocities would hold up against previous NIRSPEC observations. NIRSPEC has been used to obtain multi-epoch detections of exoplanets dating back to 2011; the first set of which were published by Lockwood et al. (2014) (Tau Boötis b). Unfortunately, there are not enough NIRSPEC epochs for any one system to be able to test the effectiveness of different primary velocity groupings. Therefore, in order to test our predictions, we combine epochs from different targets. From our archive of NIR-SPEC observations, we compile the five epochs with the primary velocities nearest zero and the five epochs with the largest absolute primary velocities (which happen to all be in the "most blue-shifted"/most negative category). These epochs are from Tau Boo b, HD187123b, 51 Peg b, and KELT2Ab. All of these planets are on orbits that can be approximated as circular.
In order to combine all epochs we need to perform a change-of-base so that the epochs reflect a single Keplerian line-of-sight orbital velocity K p . We denote all true parameters from the different systems without a prime, and all parameters of the fictitious combined system with a prime ( ). The primary (v pri ) and secondary (v sec ) velocities are encoded in the data and so cannot be altered. For a single system, v pri is variable because of the changing barycentric velocity in the direction of the system, but here, the variability in v pri can account for both changing barycentric velocities, and the different systemic velocities of the different targets. K p must be large enough to account for all the values of v sec − v pri ; we set it to 150 km/s. Then, rearranging Equation 3 for the secondary velocity, we get Table 1 gives the true parameters (K p , P , T o , v sys ) from the four target systems. Table 2 gives the information about the specific dates we are considering. Because these systems have different expected stellar and planetary spectra, we use different spectral templates to cross correlate each epoch of data, and then use the new change-of-base M values to combine the log likelihood curves generated from the two dimensional cross correlations. The data reduction and stellar and planetary spectral template used for cross correlation for each of the sources is described in the Appendix. The five log likelihood curves that make up the two primary velocity groups are shown in Figure 9. Aside from the different planetary and stellar spectral models used for each epoch, there are a few other differences between these combinations of data epochs and the predictions for the near-zero and largest absolute primary velocity groups in Section 3. First, all 10 of the data epochs were taken prior to the NIRSPEC upgrade, while the simulations considered post-upgrade NIRSPEC specifications. These differences affect the total S/N, the instrument resolution, wavelength coverage, and wavelength range covered (see Appendix). Second, the primary velocity groups are not defined as strictly here as they were in Section 3. This is simply due to the availability of data epochs. While the near-zero primary velocity group in Section 3 chose primary velocities from -2 to 2 km/s, the data near-zero primary velocity group have primary velocities ranging from -11.9 to 1.3 km/s. The simulated largest absolute (most blue-shifted/most negative) primary velocity group pulls v pri values from -30 to -28 km/s, while the data largest absolute primary velocity group ranged from -59.2 to -38.6 km/s.
In Section 4, we found no statistical difference between the Gaussian heights relative to the noise or Gaussian widths of the sets of simulations with K p values of 75 vs. 150 km/s. Therefore, the fact that these data are set up with a K p of 150 km/s should not be one of the factors differentiating these results from the K p = 75 km/s simulations. Figure 4 shows the normalized log likelihoods of the five epochs with primary velocities near zero and the five epochs with the largest absolute primary velocities. When fit with Gaussians in the same way as the analysis shown in Figures 2 and 3, the primary velocity near-zero case can be fit by a Gaussian at 157 ± 15 km/s with a height over the noise of 1.7, while the largest absolute primary velocity case fit gives a value of 180 ± 34 km/s with a height over the noise of 1.4. References.
- (1)  Note. -We assume circular orbits for all of these targets. We calculate γ values for the data combinations of epochs as 0.28 for the near-zero primary velocity group and 0.75 for the largest absolute primary velocity group. While the near-zero primary velocity case has a lower γ value, its fit is slightly more accurate and higher relative to the noise than the largest absolute primary velocity case.
We ran simulations with the exact primary velocities and M s (and γs) from the data near-zero and largest absolute primary velocity groups to determine (1) if the primary velocity groups still showed a similar trend when they were not defined as strictly as in the simulations in Section 3, (2) if the near-zero primary velocity epochs were more effective with pre-upgrade NIRSPEC settings as well as with post-upgrade NIRSPEC settings, and (3) if this trend in the data is based on the different combinations of epoch orbital phases or can be assigned to primary velocity differences. For these simulations, we only consider a single planetary and stellar spectral model, defined in Section 2.1, for each epoch. Thus, we do not expect them to appropriately reproduce the offpeak structure in Figure 4, but they should allow us to answer the questions listed above.
The results of the simulations are shown in Figure 5, with the primary velocities and orbital phases of the data in the top panel. We show these simulated log likelihood curves with separate y-axes because the structured noise in the largest absolute primary velocity simulation is at a much higher level than that in the near-zero primary velocity case. Gaussian curves find fits for the near-zero v pri and M s (in maroon) of 138 ± 81 km/s with a height relative to the noise of 4.4 and the largest absolute v pri and M s (in orange) of 36 ± 79 km/s with a height of 1.1, which would not be considered a detection as it is more than 1σ away from 150 km/s. These simulations then do indeed agree that the data near-zero primary velocity epochs have a better chance of detecting the planet. This validated that the near-zero v pri epochs are more effective with both pre-and post-upgrade NIRSPEC and that they are still preferable to larger absolute value primary velocity epochs even if not as strictly defined to −2 ≤ v pri ≤ 2 km/s. We do note that the simulations predict a much larger improvement going from the largest absolute primary velocity case to the near-zero primary velocity case (1.1 to 4.4) than was seen in the data (1.4 to 1.7). Recall that the simulated data sets are generated assuming that all non-saturated tellurics can be perfectly corrected. This is likely not the case in the real data. Because the planetary velocities are closer to the telluric frame in the near-zero primary velocity group (both because of the near-zero primary velocities themselves and because of the smaller γ, see Table 2), the planetary spectral lines are closer to the corresponding telluric lines than in the largest absolute primary velocity group. The L band wavelengths covered are dominated by water features at hot Jupiter temperatures, and while this water is much hotter than telluric water, the closer overlap between its features and the imperfectly corrected telluric water features could be responsible for hampering the near-zero primary velocity  Figure 4. The top panel has the correct grouping of primary velocities and orbital phases and the bottom panel swaps the primary velocities and orbital phases to test whether the difference in detection strengths can be said to be mostly from the different primary velocity groups or whether the combination of orbital phases also had a large effect on the detection strengths. Unlike Figures 1-3 and 6, these simulations are for pre-upgrade NIRSPEC data. case detection significance in the real data.
Using these simulations, we next test whether the improvement of the primary velocity near-zero epochs over the largest absolute primary velocity epochs was in fact due to the primary velocity differences, or if it was made by the different combinations of orbital phases. To do so, we ran simulations with the primary velocities and orbital phases of the data swapped. The results of these swapped simulations are in the lower panel of Figure 5. The simulation with the near-zero primary velocities but orbital phases of the largest absolute primary velocity epochs (in dark green) can be fit as 146 ± 9 km/s, with a relative height of 3.3 and the simulations with the largest absolute primary velocities but the near-zero orbital phases (in light green) was fit as 145 ± 10 km/s and a height over noise of 1.0. These simulations would both qualify as detections, but the one with near-zero primary velocities is stronger. Table 3 lists the Gaussian param-eters of the two data and four simulated log likelihood curves. Figure 2 saw no real trend in the Gaussian heights relative to the noise as a function of γ in any of the primary velocity groups. This is seen in the data simulations as well. The near-zero primary velocity data epochs had orbital phases corresponding to a γ of 0.28, while the largest absolute primary velocity data epochs had orbital phases corresponding to a γ of 0.75. The data near-zero primary velocities gave a Gaussian height of 4.4 when combined with the γ = 0.28 orbital phases vs. 3.3 with the γ = 0.75 epochs. The largest absolute data primary velocities showed a height relative to the noise of 1.0 with the γ = 0.28 epochs and was not detected with the γ = 0.75 epochs. These results support our finding that the primary velocities of the data epochs had a stronger effect on the detection strength than the positions of the orbital phases. It could be that the velocity separation given by epochs far from conjunction will be important though, especially for near-zero primary velocity epochs, when residual telluric features cannot be perfectly corrected from the data. Figure 2 also showed a fairly significant correlation between increasing γ and decreasing Gaussian width in the near-zero primary velocity group. This too is found in these new data simulations: the near-zero primary velocities found a Gaussian width of 81 km/s for a γ = 0.28 data set and 9 km/s for γ = 0.75.

NUMBER OF EPOCHS
The simulations thus far have all considered 5 epochs of data. We were interested to see what increasing the number of simulations would do to the detection strengths of the different primary velocity groups. To do this, we compared the random and near-zero primary velocity groups with 5, 10, and 20 epochs. We maintain a constant total S/N per pixel in a simulation across all of the epochs. The S/N per epoch then decreases with increasing epoch number from 2860 (5 epochs), to 2020 (10 epochs), to 1430 (20 epochs). Figure 6 shows the results of these epoch number simulations, with near-zero primary velocity epoch simulations shown in light blue and random primary velocity epoch simulations shown in green. These Gaussian results are plotted with respect to γ.
Of the six groups (near-zero and random primary velocity groups with 5, 10, and 20 epochs), all simulations were able to detect the planetary signal except 14 of each the 5-and 10-epoch random primary velocity simulations. The mean Gaussian heights over noise of the detected planetary signals for the near-zero primary velocity groups are 6.2 ± 1.9, 8.2 ± 2.3, and 8.9 ± 2.3, for the 5, 10, and 20 epoch simulations, respectively. For the random primary velocity epochs, the means are 2.7 ± 1.1 for the 5 epoch case, 2.6 ± 1.1 for the 10 epoch case, and 2.8 ± 0.8 for the 20 epoch case. Interestingly, we see that the two populations actually appear to diverge as the number of epochs increases, rather than converge as we had expected. The random primary velocity simulations maintain a similar Gaussian height over noise as the number of epochs increases. KS statistics tell us that the 5 and 10 epoch random primary velocity results can be assumed to have been pulled from the same distribution with very high confidence (p = 0.89). The 20 epoch random v pri group can be assumed to be pulled from a different distribution from each of the 5 and 10 epoch random v pri groups at a 3.7σ confidence level. We can see in Figure 6, and in the reported mean, that the 20 epoch simulations have much less variance in peak height over noise than the 5 and 10 epoch random v pri populations. Visually, and through the reported means, unlike the random primary velocity cases, the near-zero primary velocity heights over the noise increase from 5 to 10 epochs, and then seem to stabilize from 10 to 20 epochs. While the 5 and 10 epoch near-zero simulation heights can be said to be from different parent distributions at a 5.4σ confidence level, the 10 and 20 epoch near-zero distributions are only distinct at a 2.3σ confidence level. As we will address later, these values do not account for the fact that in addition to being generated from simulations with different numbers of epochs, these heights are apparently pulled from different γ distributions. If this were accounted for, we would expect even more commonality between the near-zero heights from 10 and 20 epoch simulations.
If we compare the Gaussian heights over the noise level from the random and near-zero primary velocity groups at 5, 10, and 20 epochs, respectively, we find that twotailed p-values measured from the Kolmogorov-Smirnov statistic decrease from 5.1 × 10 −30 with 5 epochs to 2.8 × 10 −40 with 10 epochs to 9.5 × 10 −44 with 20 epochs. While the heights of planetary detections resulting from the random and near-zero primary velocity groups can always be said to be pulled from statistically distinct populations, the level at which this claim can be made increases by orders of magnitude as the number of epochs in the simulation increases. The larger jump from the 5 to 10 epoch p-values versus the 10 to 20 epoch p-values reflects the leveling off of the peak heights over noise from 10 to 20 epochs.
As in Figure 2, the near-zero primary velocity groups widths, at every number of epochs, show a decreasing trend in both magnitude and variability with increasing γ. The random primary velocity group widths do not show this trend as strongly, likely indicative of the Gaussian models not fitting the planetary peak structure as well as in the near-zero primary velocity case.
Another thing we see in Figure 6 is that as the number of epochs increases, the distribution of γ values narrows. In fact, while the mean γ distribution across 1000 5-epoch simulations would be 0.64 ± 0.14, for 1000 10and 20-epoch simulations, it would be 0.64 ± 0.10 and 0.64 ± 0.07, respectively. We might then expect the distribution of random primary velocities to be narrowing with increasing epoch number too. To see this directly, we can define a parameter β that quantifies the combination of primary velocities in a simulation, similarly to how γ quantifies the combination of orbital phases.
We divide by 30 to normalize by the maximum absolute primary velocity of our simulation. Then, the mean of 1000 β values for 5, 10, and 20 epoch simulations would be 0.63 ± 0.14, 0.64 ± 0.09, and 0.64 ± 0.07 as well. This could explain why the random primary velocity simulations do not benefit from more epochs to the same extent that the near-zero primary velocity epochs do. As the number of epochs increases, the β distribution from which the random primary velocities are drawn is increasingly pulled away from the optimal near-zero case. While the heights of the Gaussian peaks from the nearzero primary velocity group increase from 5 to 10 epoch simulations, the 10 epoch random primary velocity simulations would have a slightly worse placement of primary velocities than the 5 epoch case. The larger number of epochs and worse placement balance each other out so that the Gaussian heights remain comparable at different epoch numbers.
From these analyses, we have seen that the average planetary detection from near-zero primary velocity epochs grew by a factor of 1.3 from 5 to 10 epoch simulations and nearly leveled out from 10 to 20 epoch simulations. The planetary detectability from random primary velocity epochs was nearly comparable at 5, 10, and 20 epochs, though there were no non-detections with 20 epoch simulations and a 14% non-detection rate at 5 and 10 epochs. As we saw above, the near-zero primary velocity epochs can be better modeled by a Gaussian than the random primary velocity epochs and as a result, their Gaussian widths show a stronger dependence on γ.

STELLAR PROPERTIES
We next investigated how the properties of the host star affect the optimal primary velocity observing strategies. Could it be that the near-zero primary velocity observing strategy works well with the 5815 K stellar model assumed because this star has strong lines corresponding to strong telluric features that are removed by telluric masking when these spectra are aligned? If so, will the near-zero primary velocity observing strategy be as effective for other stellar temperatures? To test the generalizability of the near-zero primary velocity approach, we ran 100 near-zero and 100 random primary velocity simulations with five stellar models ranging from 5200 to 7500 K, to cover the F and G spectral types. We increased the stellar radius along with the temperature, but maintained a constant metallicity and surface gravity. We also used the same planetary model in each case. Results from these simulations are outlined in Table 4 and shown in Figure 7.
These simulations reveal some interesting trends. They show that the near-zero primary velocity approach becomes even more beneficial when targeting hot Jupiters around late G-stars with lower temperatures, but less so for hot Jupiters around hotter stars. This could be due to a number of factors. Cooler stars have much more complex spectra which could allow for more improvement from well-chosen alignments of stellar and telluric features. They also have more spectral similarity with both their planets and the Earth. Water signatures arising in cooler stars could add to the need for carefully chosen, near-zero primary velocity epochs. Hotter stars, on the other hand, with fewer lines, will not be as affected by velocity shifts relative to the telluric frame.
We also see quite a similarity between the heights derived from the random primary velocity cases across stellar temperatures and radii. This shows that there is a well balanced trade-off between more complex stellar spectra at lower stellar temperatures and lower planet/star contrast at higher stellar temperatures.
These simulations indicate that the primary velocity trends observed in this work will be increasingly important to the study of hot Jupiters around cooler stars. We encourage future work into how more appropriate planet populations for each host stellar temperature and radius inform optimal high resolution, cross correlation observing strategies.

Applicability of Near-Zero v pri Observing Strategy
We have seen that epochs during which the primary velocity is near zero, or, equivalently, the systemic velocity is canceled as much as possible by the barycentric velocity in the direction of the system so there is very little velocity separation between the telluric and stellar spectra, provide the strongest planetary detections. Not all systems will have periods during which the primary velocity is near zero however. The magnitude of the barycentric velocity, determined by a system's right ascension and declination, must be large enough to cancel out its systemic velocity. The second Gaia data release DR2 (Gaia Collaboration et al. 2016 published radial velocities averaged over 22 months from 7,224,631 stars. These reported radial velocities were all from sources brighter than G RV S = 14 (the flux measured in the Radial Velocity Spectrometer G band); with a fraction of transits where the source was detect as having a double-lined spectrum less than 0.1 (to remove detected double-lined spectroscopic binaries); with an uncertainty on the radial velocity below 20 km/s; and a spectral template used to derive the radial velocity with an effective temperature from 3550 to 6900 K (Gaia Collaboration et al. 2018). This was a substantial and collaborative effort and many researchers contributed to this impressive radial velocity data set (e.g., Cropper et al. 2018;Sartoretti et al. 2018;Soubiran et al. 2018;Katz et al. 2019).
Of the 7,224,631 stars with radial velocities reported in Gaia DR2,3,209,212,or 44.4%, have combinations of locations and systemic velocities that will allow for a near-zero primary velocity at some point during the year. This then suggests that our predicted optimal nearzero primary velocity strategy will then be applicable to nearly half of the planetary systems in the sky.

Factors Influencing the Random Primary Velocity
Distribution In this work, to define the random primary velocity distribution, we pulled uniformly from the Earth's orbit, converted the Earth's position to a barycentric velocity assuming an orbital motion of 30 km/s, and then converted to primary velocity assuming a systemic velocity of 0. This resulted in a bimodal primary velocity distribution that was relatively flat through the center and rose quickly in probability toward ±30 km/s (Figure 8). It was this distribution from which we pulled "random" primary velocities.
This distribution could be different for different systems depending on their right ascensions, declinations, and systemic velocities, though. Their RAs and DECs will determine the component of the Earth's orbital motion that is in the line-of-sight to the system. While 30 km/s is about the Earth's actual orbital velocity, the barycentric velocity in the direction of a system would only vary from 30 to -30 km/s if the system were precisely on the Earth's orbital plane. Anywhere else and the range of possible barycentric velocities would shrink. A smaller range of primary velocities -if the systemic velocity was still 0 km/s and so the primary velocity range still centered around 0 km/s -might allow for slightly stronger detections than the random primary velocity simulations shown in this work. This would be because the random primary velocity distribution would have a smaller horizontal extent, so the values pulled would be nearer the optimal 0 km/s. Systems with non-zero systemic velocities would also have a differently shaped random primary velocity distribution than the one used in this work. The center of the distribution would now be around the systemic velocity, rather than zero, and the sharp increases in probability would be at the minimum (v sys − max(v bary )) and maximum (v sys − min(v bary )) ends of the distribution. We would expect systems with systemic velocities canceled by the maximum or minimum possible barycentric velocities along their line-of-sight to show the strongest detections possible with a random primary velocity sampling strategy. This would be because one of the sharp increases in towards the edges of the primary velocity distribution would be at 0 km/s, meaning that a random primary velocity strategy would offer the most near-zero primary velocity epochs if the system had this configuration of systemic and barycentric velocities (set by its RA and DEC) than any other alignment.
Another consideration for the random primary velocity observing strategy will be the location of the telescope. The location of the telescope will set if and when 2.3 ± 0.9 3.2 × 10 −4 3.6 a The p-value and nσ come from two-tailed KS tests run between the near-zero and random primary velocity simulations with the same stellar properties. nσ refers to the level at which we can reject that the two groups are pulled from the same parent distribution. Below a p-value of ∼ 10 −16 , the confidence level that the populations are distinct approaches ∞.  Table 4. "Random" primary velocity distribution generated from 10,000 uniform pulls of Earth's orbital position.
during the year a system is observable, and so will cut out portions of the random primary velocity distribution corresponding to barycentric velocities that arise during times the system is not observable from that telescope. The primary velocities removed could include a range around zero or a range of the largest possible absolute primary velocities. If the primary velocities around 0 were removed, we would expect a weaker detection, while if a range of the largest absolute primary velocities were removed, we would expect a stronger detection.
We do note that these predictions are all based on which selection of random primary velocities would give the most near zero. We would still expect a dedicated near-zero primary velocity observing strategy to provide the strongest detection because it would not be diluted by any of the suboptimal non-zero primary velocity epochs that could arise from a random primary velocity sampling strategy regardless of systemic velocity, magnitude of the barycentric velocity variation, or telescope location. If the combination of the telescope location and systemic and barycentric velocities are such that there is no period of near-zero primary velocities, we would recommend targeting the smallest absolute (nearest zero) primary velocities as we saw provided stronger results in our pre-upgrade NIRSPEC simulations of Section 5.

Random Orbital Phases
In this work, we found that the combination of orbital phases did not have a large effect on the height over noise of the planetary detection. Our simulations considered cloud-free models and did not vary the planetary spectrum as a function of orbital phase, to account for day-to night-side differences for tidally locked planets, though. If day-to night-side differences were considered, we would expect the day-side orbital phases (0.25 ≤ M ≤ 0.75), which should have higher effective temperatures, to allow for stronger detections (Finnerty et al. 2021).
While clouds have presented a challenge to lowresolution transmission spectroscopy, thermal emission spectra of the same planets show strong molecular lines (e.g., Crouzet et al. 2014;Morley et al. 2017). Gandhi et al. (2020) recently showed that high resolution transmission spectroscopy could be used to detect water and other trace species, namely CH 4 , NH 3 , and CO, in cloudy atmospheres with a modest observing time from a ground-based telescope. In high resolution emission spectra, clouds could decrease the line contrast by shifting the continuum to higher altitudes and lower tempera-tures, rather than by blocking stellar rays below the cloud tops as they do in transmission spectra. By decreasing line contrasts, clouds would make the planet more difficult to detect through cross correlation analysis. If, in tidally-locked atmospheres, the clouds are mainly constrained to the night-side (e.g., Demory et al. 2013;Parmentier et al. 2016), day-side epochs would be even more preferable.
For longer period planets that are not tidally-locked, neither day-to night-side temperature differences nor night-side clouds would uniformly degrade one set of orbital phases over another. Additionally, neither day-to night-side differences nor the presence of clouds should affect our predictions for the optimal primary velocity observing strategy.
Importantly, the fact that random orbital phases are sufficent, at least for non-tidally locked atmospheres, indicates that a robust detection could be made with only a fraction of an exoplanet's orbital period. Short period planets could be well detected with a selection of orbital phase epochs taken over a period when v pri is near zero. With the much more quickly varying planetary orbital phase relative to Earth's orbital phase, these periods during which v pri is near zero should offer a range of dayto night-side planetary epochs. Longer period planets could be targeted at multiple stretches when v pri is near zero, each offering a different selection of orbital phases (as long as the orbital period is not highly commensurate with that of the Earth). Such observing strategies could be easily obtainable and should lead to strong (nontransiting) planetary detections.

Wavelength Dependence and Atmospheric
Characterization Further, while we investigated ways to strengthen detections of planetary emission through the recovery of K p in this work, ultimately, we would be interested in constraining various planetary atmospheric properties, such as the presence and relative abundances of various molecular species and the natures of the atmospheric thermal structure, winds, and planetary rotation. Previous work has found that since there are no spectral lines from major carbon-bearing species in the L band of hot Jupiter atmospheres, this data alone is not sufficient to constrain their atmospheric C/O ratios (e.g., Piskorz et al. 2018;Finnerty et al. 2021). Such measurements may be possible for warm Jupiters (T eff ≈ 900 K) from L band data alone however. At cooler effective temperatures, sufficient methane can be expected under equilibrium conditions to be detectable in L band data. With both methane and water appearing, L band data can provide constraints on the C/O ratios of warm Jupiters (Finnerty et al. 2021).
To make these C/O constraints for hot Jupiters would likely require additional epochs in the K or M bands, where prominent carbon monoxide bandheads exist. In this work, we found that both pre-and post-upgrade NIRSPEC L band simulations were better able to detect planetary signals with near-zero, rather than with random, primary velocity epochs. The pre-and postupgrade simulations differ in both number of orders per epoch and order wavelength coverage, with no overlap between the wavelengths covered. The fact that both still preferred a near-zero primary velocity epoch strat-egy implies that these predictions are not completely wavelength dependent and we expect that they should hold for L band observations in general. We encourage more simulation work to determine how widely generalizable these predictions will be both at other NIRSPEC bands (specifically K and M ) and across the large instantaneous spectral grasp promised by upcoming and proposed instruments such as GMTNIRS (1.1-5.3 µm) and IGNIS (1-5 µm). Data covering these multiple bands would allow us to detect carbon monoxide as well as water in hot Jupiter atmospheres and allow for constraints on their atmospheric C/O ratios.

CONCLUSION
In this work, we aimed to determine how to best strengthen planetary detections and reduce structured noise in few epoch data sets with careful observing strategies. The two key parameters that can be selected with the choice of observing nights are the primary velocity (because of the variable barycentric velocity) and the planetary orbital phase. We found that epochs taken during nights when the primary velocity of the system is near 0 km/s, so that there is very little relative velocity shifting of the stellar and telluric reference frames, will provide the strongest planetary detections. With a random selection of planetary orbital phases, these nearzero primary velocity epoch simulations produce planetary peaks more than two times higher relative to the noise than simulations generated with randomly selected primary velocities. Further, for near-zero primary velocity epochs, the closer their orbital phases are to quadrature, the better the constraints on K p will be. Following these results, we recommend that observers looking to build up multi-epoch near-IR high resolution data sets target, first, epochs with near-zero primary velocities, and second, epochs with orbital phases near quadrature to get the best constraints on the planetary detection.
In this work, we demonstrated how greatly the combinations of primary velocities and orbital phases can affect a planetary detection. Moving forward, careful attention should be paid to planning observations, and all few epoch data sets should not be assumed to have an equal probability of detecting a planet. Following these predications, observations taken from upcoming multi-echelle instruments, such as GMTNIRS and IGNIS, during periods when the primary velocity of a system is near zero could provide both robust detections of exoplanets in a fraction of their orbital periods and constraints on their atmospheric composition.
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
We would also like to thank an anonymous reviewer who pointed us in interesting new directions and by doing so improved the content of this paper.

Notes on Epochs from Individual Sources
All of the data used in Section 5 is L band data from the pre-upgrade NIRSPEC instrument. Each epoch has 4 orders, covering approximately 2.9962-3.0427, 3.1203-3.1687, 3.2552-3.3058, and 3.4026-3.4554 µm. The average spectral resolution is 20,000, and the total S/N across the 5 epochs is about 4100.

HD187123
The reduced data and PHOENIX stellar and SCARLET planetary spectral models used here were those presented in Buzard et al. (2020).

KELT2A
The reduced data and PHOENIX stellar and planetary models used here were those presented in Piskorz et al. (2018). For the planetary model, we used the best fitting ScCHIMERA model, which had a metallicity (log z) of 1.5, a C/O ratio of 0.5, and an incident solar flux f of 1.0. This parameter f accounts for day-night heat transport and an unknown albedo by scaling a wavelength-dependent incident stellar flux (from a PHOENIX stellar grid model).
Defined this way, model atmospheres with f 1.5 show a temperature inversion.

Peg
The 51 Peg epochs were reduced in the same way as the other epochs (e.g., Piskorz et al. 2018;Buzard et al. 2020), and telluric corrected through a Molecfit (Kausch et al. 2014) guided principal component analysis. We used a PHOENIX stellar spectral model interpolated to an effective temperature of 5787 K, a metallicity of 0.2, and a surface gravity of 4.449 (Turnbull 2015). The planetary spectral model we use was generated from the SCARLET framework. It does not have an inverted thermal structure, which was suggested is appropriate by Birkby et al. (2017).

Tau Boo
The Tau Boo data used here were processed using a Molecfit initial telluric model followed by PCA to remove residual tellurics. The stellar and planetary spectral model used here were the ones used in Lockwood et al. (2014). The stellar model was not from the PHOENIX framework. Rather, it was generated from the LTE line analysis code MOOG (Sneden 1973) and the MARCS grid of stellar atmospheres (Gustafsson et al. 2008). Individual elemental abundances were set through fitting to well measured lines in the NIRSPEC data. See Lockwood et al. (2014) for a full description of the stellar spectral model generation.  Table 1. These curves, converted to Kp space and summed, make up Figure 4. In each subplot, the red dashed line corresponds to the primary velocity at that epoch and the black dashed line corresponds to the vsec given by a K p of 150 km/s at that each. If the fictitious combined system were face-on, with a Kp of 0 km/s, the black dashed line would coincide with the red dashed line. If, on the other hand, it were edge-on, vsec would fall on the other end of the white range of possible planetary velocities. Here we have the maximum value of Kp arbitrarily set to 230 km/s.