Reinvestigation of the Multi-Epoch Direct Detections of HD 88133 b and Upsilon Andromedae b

We reanalyze the multi-epoch direct detections of HD 88133 b and ups And b that were published in Piskorz et al. 2016 and Piskorz et al. 2017, respectively. Using simulations to attempt to reproduce the detections, we find that with the 6 and 7 $L$ band Keck/NIRSPEC epochs analyzed in the original works, the planets would not have been detectable unless they had unreasonably large radii. HD88133 and ups And both have fairly large stellar radii, which contributed to the difficulty in detecting the planets. We take this opportunity to consider how these planets may have been detectable with the small number of epochs originally presented by running simulations both with the upgraded NIRSPEC instrument and with near-zero primary velocities, as recommended by Buzard et al. 2021. While 7 $L$ band NIRSPEC2.0 epochs with near-zero primary velocities would have allowed a strong ($10.8\sigma$) detection of ups And b, many more than 6 $L$ band epochs would have been required for a strong detection of HD88133b, which could be due in part to both this system's large stellar radius and low stellar temperature. This work stresses the importance of careful analytic procedures and the usefulness of simulations in understanding the expected sensitivity of high-resolution spectroscopic data.


INTRODUCTION
Direct detection techniques that use radial velocity signatures from exoplanet orbital motion to detect their atmospheric thermal emission have become popular in the last decade. Two variations of such high-resolution techniques have been developed. Both aim to make a direct planetary detection through a measurement of the planetary velocity semi-amplitude, K p . One variation targets planetary systems during periods when the change in the planetary line-of-sight motion is the greatest, typically near conjunction. By observing the system at periods of maximum planetary line-of-sight acceleration, the data contain planetary signatures that shift across the instrument's resolution elements over the course of the night, and techniques like principal component analysis (PCA) can be used to tease apart the changing planetary signal and the stationary telluric and stellar signals. A onedimensional cross correlation routine can then be used to measure the planetary velocity. This technique was first introduced by Snellen et al. 2010 with VLT/CRIRES, and has since been applied to data from a range of instruments including Subaru/HDS (e.g., Nugroho et al. 2017), TNG/GIANO (e.g., Brogi et al. 2018;Guilluy et al. 2019), and CFHT/SPIRou (e.g., Pelletier et al. 2021).
The second variation of the high resolution technique, which is the focus of this work, instead limits observations so that the change in the planetary line-of-sight velocity is minimized, and the planetary spectrum does not shift across the detector during the course of an observation. This variation is more technically challenging because there is no longer a velocity variation that can be leveraged to separate the planetary and stellar spectra in a single epoch. After the data is telluric corrected, a two-dimensional cross correlation is relied upon to pull apart the stellar and planetary components. Since the planetary signal is so much fainter than the stellar signal, multiple epochs must be combined before the planetary signal becomes apparent. While technically challenging, this variation is the only currently viable highresolution method for studying the atmospheres of planets whose semi-major axes preclude both single epoch spectroscopic detection, because they move too slowly, and direct imaging with current adaptive optics capabilities, because they are too close to the star ( 0.1", e.g., Snellen et al. 2014). This gap includes planets in K dwarf and solar habitable zones. As the so-called multi-epoch technique is uniquely capable of directly studying the atmospheres of non-transiting planets in these systems, it deserves careful work and attention.
To date, the multi-epoch technique has mainly been applied to data from Keck/NIRSPEC, which is an echelle spectrograph that offered 4-6 orders in the K and L bands per cross disperser setting and R ∼ 25, 000 − 30, 000 before its upgrade in early 2019. The method was first applied to Tau Boötis b and was able to measure its K p as 111 ± 5 km/s (Lockwood et al. 2014), which was in good agreement with K p measurements from other techniques (e.g., 110.0 ± 3.2 km/s, Brogi et al. 2012). Subsequently, the non-transiting hot Jupiters HD 88133 b and Upsilon Andromedae b (ups And b) were detected at 40 ± 15 km/s (Piskorz et al. 2016) and 55 ± 9 km/s (Piskorz et al. 2017), respectively. These planets have yet to be studied via a different technique. Piskorz et al. 2018 then detected the transiting hot Jupiter KELT-2Ab with a K p of 148 ± 7 km/s, which was in good agreement with the transit measurement of 145 +9 −8 km/s (Beatty et al. 2012). Finally, Buzard et al. 2020 measured K p of the non-transiting hot Jupiter HD 187123 b to be 53 ± 13 km/s. This detection used simulations to identify sources of non-random noise and elucidate the true planetary detection. HD 187123 b has not to date been detected via another technique.
In this work, we look back on the multi-epoch detections of HD 88133 b (Piskorz et al. 2016) and ups And b (Piskorz et al. 2017). Piskorz et al. 2016 reported the Keplerian orbital velocity of HD 88133 b as 40 ± 15 km/s using 6 epochs of NIRSPEC L band data and 3 epochs of K band data. Piskorz et al. 2017 reported the Keplerian orbital velocity of ups And b as 55 ± 9 km/s using 7 epochs of NIRSPEC L band data, 3 epochs of K l band data covering the left-hand half of the NIRSPEC detector, and 3 epochs of K r band data covering the righthand half of the detector. In this work, we will focus on the L band data because the L band data provided the majority of the overall structure in both the HD 88133 b and ups And b detections.

STANDARD MULTI-EPOCH ANALYTIC APPROACH
To begin, we want to give a brief description of the multi-epoch analytic process. These approaches are explained in more detail in prior publications (e.g., Lockwood et al. 2014;Piskorz et al. 2016;Buzard et al. 2020). In brief, epochs of data are obtained from hot Jupiter systems over ∼2-3 hour periods during which the planetary signal is not expected to significantly shift compared to the wavelength scale of the detector. The two-dimensional echelle spectra are reduced, wavelength calibrated, telluric corrected, and run through a twodimensional cross correlation routine with appropriate stellar and planetary spectral models. Because the stellar signal is the major component of the data after telluric correction, the known stellar velocity, given by where v sys is the systemic velocity and v bary is the barycentric velocity, can always be correctly measured in each epoch. A cut along the known stellar velocity gives a one-dimensional cross correlation in terms of planetary velocity shift. With a very low contrast relative to the stellar signal, the planetary signal requires the combination of multiple epochs to become clearly detectable. To be combined, the cross correlations must first undergo two transitions. They must first be converted from functions of secondary velocity, which is dependent on orbital phase, to functions of a parameter independent of orbital phase, namely, the Keplerian orbital velocity. Second, they must be converted to log likelihoods. To convert them from functions of secondary velocity to Keplerian orbital velocity (K p ), we apply the equation, v sec (f ) = −K p (cos(f + ω st ) + e cos(ω st )) + v pri (2) where f is the planet's true anomaly at the observation time, ω st is the argument of periastron of the star's orbit measured from the ascending node (with the Z-axis pointing away from the observer, see Fulton et al. 2018), and e is the eccentricity.
To convert the cross correlations from v sec to K p space using Equation 2, we need stellar radial velocity (RV) parameters (e, ω st ) and true anomaly (f ) values at each epoch. Stellar radial velocity parameters can typically be found in the literature (e.g., Butler et al. 2006). We note that it is important that the stellar orbital parameters (e, ω st ) as well as those used to calculate f (t peri , P ) are pulled from the same literature source. This will be especially important for near-circular orbits where pericenter is not well defined because though references can set pericenter at vastly different points on the orbit, their other parameters (mainly t peri and ω) would then all be consistent to that chosen point of pericenter. A t peri and ω from different references could be referring to very different points on the orbit, and so could create a large error in the derived f s and secondary velocities.
The true anomalies can be calculated using the following equations, which are described in Murray & Dermott 1999. First, the mean anomaly (M ) is calculated from the observation time (t obs ), and the stellar radial velocity parameters, time of periastron (t peri ) and orbital period (P ). If the orbit under study can be assumed circular, the mean anomalies can be used in place of the true anomalies.
Then, the eccentric anomaly E can be calculated as follows, where e is the eccentricity. As this equation does not have a closed-form solution for E given M , E is calculated numerically.
Finally, the true anomaly f is calculated as We want to make a few important notes about Equations 2 and 3-5. First, the negative sign at the start of Equation 2, which is not present in the corresponding equation in Fulton et al. 2018, allows this equation to describe the planetary motion, rather than the stellar motion. At any given time, the planetary and stellar motions should have opposite signs. The addition of the negative sign would be equivalent to replacing ω st with ω pl (the argument of periastron of the planet's orbit measured from the ascending node) because ω pl = ω st + π (for a set direction of the Z-axis). Importantly, by defining K p this way, we specify that it must be a positive value. Second, it is important that the zero-point used to define the anomalies (t peri in Equation 3) is consistent with the offset used in Equation 2. In the equations as written, f is measured from pericenter, and adding ω st in Equation 2 brings the zero point from pericenter to the star's ascending node. The ascending node (for a circular orbit, when the star is at quadrature and moving away from the observer or when the planet is at quadrature and moving toward the observer) should have the largest negative v sec possible. The negative cosine ensures this is the case. Other works that have assumed a circular orbit (e.g. Piskorz et al. 2018;Buzard et al. 2020) use a positive sine equation with no added phase offset and with M s centered at inferior conjunction. This is also a valid approach because the offset between the M s and the sine equation are consistent.
Once the planetary cross correlations from different epochs are on the same K p axis, they must be converted into log likelihoods to be combined. There are a number of ways to do so (Zucker 2003;Brogi & Line 2019;Buzard et al. 2020). Here, we use the approach first introduced by Zucker 2003 and termed the "Zucker ML" approach by Buzard et al. 2020 to be consistent with Piskorz et al. 2016 andPiskorz et al. 2017. The Zucker ML method converts cross correlations to log likelihoods and combines them as, where s is the velocity shift, R i are the individual cross correlations, N i is the number of pixels per cross correlation, and N tot is the total number of pixels.
With this basis, we turn to a reinvestigation of previous multi-epoch works by Piskorz et al. 2016 andPiskorz et al. 2017.

HD 88133 B
Piskorz et al. 2016 used 6 Keck/NIRSPEC L band and 3 K band epochs to measure the planetary Keplerian orbital velocity of the non-transiting hot Jupiter HD 88133 b as 40 ± 15 km/s. To do so, they used orbital parameters from their own fit to stellar RV data. This stellar RV data set consisted of 55 RV points; 17 had been previously analyzed by Fischer et al. 2005 and the rest were new RV points taken by the California Planet Survey with HIRES at the W. M. Keck Observatory. Fitting these data with a Markov Chain Monte Carlo technique following Bryan et al. 2016, they reported the orbital parameters (velocity semi-amplitude K = 32.9 ± 1.03 km/s, period P = 3.4148674 +4.57e−05 3.1. Correction to Piskorz et al. 2016Results In Piskorz et al. 2016 there was a systematic error in the implementation of the time of periastron in Equation 3 that resulted in a 38.0% orbital offset in the mean anomalies. The anomalies used in the paper analysis and the corrected anomalies are listed in Table 1. In Figure 1, we show the corrected log likelihood curves along with the originally published curves, analyzed with both the SCARLET and the PHOENIX planetary models. This correction causes a drastic difference in the resulting log likelihood curve. In each subplot, the red curve was the published log(L) curve and the black dashed curve is the curve reproduced with the systematic offset. Interestingly, this curve exactly reproduces the SCARLET model function, but does not reproduce the PHOENIX model function. In fact, when the data are analyzed with the PHOENIX model, Piskorz et al. 2016 orbital parameters, and the offset, the resulting function appears much more similar in shape to the corresponding SCARLET result than was originally published in a Note, the values here are expressed from 0 to 2π, rather than from 0 to 1 as in Piskorz et al. 2016.
The M and f values reported in their Table 3 also differ from these values because while they too were affected by the systematic offset, they used Butler et al. 2006 orbital parameters rather than the newly fit parameters from Piskorz et al. 2016. Piskorz et al. 2016. This suggests that the two planetary models are much more comparable than was originally thought.
The blue and orange curves in both subplots show the corrected log likelihoods when the orbit is either treated as eccentric (blue) or assumed to be circular (orange). The similarity between these curves in both subplots shows that for low eccentricity orbits (here e ∼ 0.05) circular approximations do not greatly affect the shape of the resulting log likelihood surface.
We also note that when the corrected log likelihood curves show no significant peaks at positive values of K p . Remember that given how we defined the relationship between v sec and K p , only positive values of K p are physically meaningful. Negative values of K p would have the stellar and planetary radial velocity curves perfectly in phase rather than out of phase as they should be. This correction therefore implies that we cannot report a measurement of the Keplerian orbital velocity of HD88133b of 40 km/s from the six L band epochs presented in Piskorz et al. 2016.

HD88133 Simulations
We ran a few sets of simulations to determine what physical or observational factors would allow for the detection of HD88133b. We note that HD88133A has a rather large stellar radius of 2.20 R (Ment et al. 2018), meaning that this system has an even lower planet to star contrast than most hot Jupiters, though we cannot measure it directly because of the unknown planetary radius. Our first simulations ask how large the planetary radius would have to be for the planetary peak to be detectable in the six L band epoch observed (Section 3.2.2). Next, we investigate whether the same epochs taken with the upgraded NIRSPEC2.0 instrument would have been more successful than the NIRSPEC1.0 epochs (Section 3.2.3). Then, we consider the primary velocity. Buzard et al. 2021 found that epochs with near-zero primary velocities were more useful in damping down structured noise and revealing true planetary signatures than epochs with larger absolute primary velocities. We consider whether the observed epochs would have better revealed the planetary peak if they had smaller absolute primary velocities (Section 3.2.4).  In blue and orange are the corrected log likelihood curves considering an eccentric orbit (blue) and a circular orbit (orange). (B) Same as Panel A, except using a PHOENIX planetary model rather than a SCARLET planetary model. Note here that we were unable to reproduce the published curve. However, when the same orbital parameters are used, we get a much more similar result to that of the SCARLET models than was shown in the Piskorz et al. 2016. The two different planetary spectral model frameworks do not create as significant a difference as we thought.
For these simulations, we use a stellar model generated from the PHOENIX stellar spectral model grid (Husser et al. 2013) interpolated to an effective temperature of 5438 K, a metallicity of 0.330, and a surface gravity of 3.94 (Mortier et al. 2013).
We use the PHOENIX planetary spectral model from Piskorz et al. 2016. This modeled atmosphere does not have an inverted thermal structure in regions close to the molecular photosphere.
We generated the simulated multi-epoch data using the same framework initially presented in Buzard et al. 2020 and so a full description can be found there. As a quick summary, these simulations combine the stellar and planetary models on the planetary wavelength axis after scaling them by their relative surface areas and shifting them to the desired primary and secondary velocities during each epoch. The secondary velocities are calculated by plugging the f values in Table 1, This work, into Equation 2 using Piskorz et al. 2016 orbital parameters and a K p of 40 km/s. For these simulations, we assume a planetary radius of 1 R Jup (except in the planetary radius simulations) and a stellar radius of 2.20 R (Ment et al. 2018). After combination, the continuum is removed using a third order polynomial fit from 2.8 to 4.0 µm in wavenumber space. Then, the simulated data are broadened using the instrumental profiles fit to the data, interpolated onto the data wavelength axes, and the same pixels lost to saturated tellurics in the data are removed. The data, all taken before the NIRSPEC upgrade in early 2019, contain 4 orders per epoch which cover approx-imately 3.4038-3.4565, 3.2567-3.3069, 3.1216-3.1698, and 2.997-3.044 µm. Gaussian white noise is added in at the same level as in the data (total S/N per pixel = 5352).
The stellar model used to generate and cross correlate the simulated data differs from the stellar model used to cross correlate the observed data in Piskorz et al. 2016 in how its continuum is removed. The stellar model used for cross correlation in Piskorz et al. 2016 was continuum corrected with a second order polynomial fit from 2.0 to 3.5 µm in wavenumber space, while the model used here is corrected with a third order polynomial fit from 2.8 to 4.0 µm in wavenumber space. The method of stellar continuum correction actually has a large effect on the shape of the resulting log likelihood curve; when the data are cross correlated with a stellar model corrected by a 3 rd order polynomial fit from 2.8 − 4.0 µm in wavenumber space, the resulting log likelihood curve much more closely resembles the simulated curves (e.g. in Fig. 2). We use the 3 rd order, 2.8 − 4.0 µm approach in our simulations because this continuum correction method was validated in allowing common structure to be seen in the data and simulated log likelihoods of HD187123b (Buzard et al. 2020). Further, this approach resulted in a flatter looking stellar spectral model, and one that found the known stellar velocities in each epoch of data with higher likelihoods. We do note, however, that the seemingly strong dependence of the final log likelihood shape on the method of stellar model continuum correction is concerning and warrants deeper investigation.

Planetary Radius Simulations
Because there seems to be no clear peak at a positive K p in Figure 1, we use simulations to see how much larger the planetary radius would have to be for its peak to be distinguishable. For these simulations, we set the K p at 40 km/s and the stellar radius at 2.20 R (Ment et al. 2018), and step the planetary radius up from 1 R Jup to 4 R Jup in increments of 0.5 R Jup . Figure 2 shows the results of these simulations in the top panel. In the bottom panel, we show the detections that could be made if all of the structured noise (e.g. the R pl = 0 result) could be effectively removed from the results containing a planetary signal.
In this high S/N per epoch regime, we expect the contribution from structured noise to far outway the contribution from random noise (Buzard et al. 2020;Finnerty et al. 2021). The similarity between different radius simulations shows that this is still the case.
To quantify the strength of these detection peaks, we fit each with a Gaussian model and report the parameters in Table 2. In the non-corrected versions, the Gaussian model does not fit within one standard deviation of the input K p until the planetary radius reaches 2.5 R Jup , and the peak does not exceed 3σ until a radius of 3.5 R Jup . Even at this large radius, the Gaussian model does not clearly distinguish the planetary peak from the structured noise, which can be seen in the Gaussian center offset, large Gaussian width, and relatively low R 2 value. Much more promising detections could be made if there were a way to effectively remove the structured noise from the log likelihood results. Even at 1 R Jup , HD 88133 b could have been detected in the data with a significance over 3σ. These simulations take a num-ber of liberties, though, that are not yet applicable to real data. They consider no telluric contamination outside of pixels lost to saturated telluric absorption. They also assume that the stellar and planetary spectra in the data are perfectly matched to the stellar and planetary templates used for cross correlation. Thus, while the corrected simulations shown in the bottom panel of Figure 2 provide an optimistic view of the possible detections with the 6 particular NIRSPEC1.0 L band epochs presented in Piskorz et al. 2016, the uncorrected versions give much more realistic estimates.
HD 88133 b has a minimum mass of 0.27 ± 0.01 M Jup (Piskorz et al. 2016). With a radius of 3.5 R Jup , it would have a minimum density of 0.01 g/cm 3 . A growing classifcation of planets with exceptionally large radii for their masses, called "super-puffs," have low densities of 0.3 g/cm 3 (e.g. Cochran et al. 2011;Jontof-Hutter et al. 2014;Vissapragada et al. 2020). While HD 88133 b (M p sin i ∼ 85M ⊕ ) is too massive to be classified as a super-puff (M p 10−15M ⊕ , Piro & Vissapragada 2020), by comparison of its density, we can conclude that it is highly improbably the planet's radius would be as high as 3.5 R Jup . Indeed, hot Jupiter inflation can approach 2 R Jup (Thorngren & Fortney 2018), but has not been observed to exceed it to this extent.
Our simulations therefore confirm that HD 88133 b is not detectable from the six L band epochs of data presented in Piskorz et al. 2016. These radius simulations did, however, provide useful information in telling us that the planetary signal would need to be raised by about an order of magnitude (or the structured noise lowered by the same amount), to allow for a confident detection. We now turn to simulations to ask how that order of magnitude may be made up observationally rather than by altering parameters of the physical system like the planetary radius.

Upgraded NIRSPEC Simulations
The NIRSPEC instrument was upgraded in early 2019, after the Piskorz et al. 2016 publication. The upgrade would afford 6 usable L band orders per echelle/cross disperser setting as opposed to the 4 from NIRSPEC1.0. It would give twice as many pixels per order (2048 vs. 1024), a nearly doubled spectral resolution (∼41,000 vs. 25,000), and a ∼40% larger wavelength coverage per order (Martin et al. 2018).
We run NIRSPEC2.0 simulations with the same orbital phases and primary velocities in the original six NIRSPEC1.0 epochs to determine whether the instrument upgrade would have given the planetary signal the boost it needed to be detectable. These simulations are generated as described in Section 3.2.1, with the following exceptions. The six orders per epoch cover roughly 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. We broaden the simulated data to an average instrumental resolution of 41,000 and assume a S/N per pixel per epoch of 2860, or a total S/N per pixel of 7000 across the six epochs. The data wavelength axes, locations of saturated telluric pixels, and Gaussian kernals used to broaden the simulated data were taken from the 2019 Apr 3 and 2019 Apr 8 NIRSPEC2.0 data of HD187123 presented in Buzard et al. 2020. Figure 3 shows the results of the simulations which   Note.
-These simulations were all run with an input Kp of 40 km/s. Prior to fitting, these log likelihood results are subtracted by the mean of their values from -150 to 0 km/s. The Gaussian model is initiated with a 40 km/s center and 10 km/s standard deviation. The Gaussian peak height is reported over σ, which is measured as the standard deviation of the log likelihood structure more than 3∆Kp above or below the best-fit Gaussian center, where ∆Kp is the standard deviation of the best-fit Gaussian model.  consider imagine the HD88133 L band epochs had been taken with the upgraded NIRSPEC instrument. In light purple is the simulation with no planetary signal added and in darker purple is the simulation with a 1 R Jup planetary signal. While the likelihood at the input K p of 40 km/s is increased from the corresponding no planet log likelihood, it does not form a peak and would not constitute a detection. The six L band HD88133 epochs were positioned such that even with the upgraded instrument, they would not have enabled a planetary detection.

Near Zero Primary Velocity Simulations
Buzard et al. 2021 recently showed that, in the small epoch number limit, epochs taken when the primary velocity of a system is near zero are better at reducing structured noise and revealing planetary signatures than epochs taken during periods with larger absolute primary velocities. The majority of the structured noise that arises in the simulations presented here and in Buzard et al. 2021 results from correlation between the planetary spectral template and the stellar component of the simulated data. We thus suspect that the reduction of structured noise at a primary velocity of zero relates to the portion of the stellar spectrum masked by saturated tellurics when there is a minimal velocity shift between the two spectra. It could be that at this velocity shift, the stellar features that most strongly correlate with the planetary template are masked by saturated tellurics, and without them, the structured noise level decreases substantially. One must be careful in applying this prediction, though, since a primary velocity of zero would bring not just the stellar spectrum, but also the planetary spectrum, closer to the telluric rest frame. While our simulations assume perfect correction of non-saturated tellurics, any residual tellurics that make it through a correction routine could mask planetary features. In a small epoch number limit, an optimal routine might therefore include near-zero primary velocities (to reduce structured noise from star/planet correlation) and quadrature orbital positions (to maximize the planet/telluric velocity separation). With a much larger number of epochs, the structured noise from planet/star correlation may be reduced naturally by the many different combinations of primary and secondary velocities and the usefulness of near-zero primary velocity epochs may be lessened.
We can ask whether, with just the six epochs on HD88133, near-zero primary velocities might have helped. HD88133 has a primary velocity of zero twice a year: in late February and in mid-August. It would also be accessible from Keck in late February. For the following simulations, we assume epochs had been taken with the same orbital phases (f ) as the original data, but in late February when v pri = 0 km/s. The original data epochs had primary velocities of 17.4, 18.1, 8.1, 16.2, 25.7, and 19.5 km/s. We run these simulations with both the NIRSPEC1.0 and NIRSPEC2.0 configurations. Figure 4 shows the results of the 0 primary velocity simulations with NIRSPEC1.0 results in the top panel and NIRSPEC2.0 results in the bottom panel. In each, we show a pure structured noise simulation (light purple), or simulation with no planetary signal in the simulated data, so that the planetary peak in the simulation with the planetary signal (darker purple) can be distinguished from the structured noise.
We first consider the NIRSPEC1.0 simulation. While the 1 R Jup planetary signal definitely shows a larger peak here than when analyzed with the original primary velocities (Fig. 2), it still does not constitute a very strong detection. We can think of a number of reasons for this. Buzard et al. 2021 showed that near-zero primary velocity epochs could raise the detection significance on average ∼ 2 − 3× over random primary velocity epochs. From our radius simulations, we estimate an order of magnitude is needed. The gain from near-zero primary velocities then may not be sufficient. HD88133A has an effective temperature of 5438 K (Mortier et al. 2013), putting it on the cooler end of host stars considered by Buzard et al. 2021. Cooler host stars showed a stronger preference for near-zero primary velocity epochs, which means in this case, we might expect a bit more than a 3× increase. On the other hand, here, we consider a K p of 40 km/s, smaller than the 75 km/s K p used for most of the simulations in Buzard et al. 2021. A smaller K p brings all of the secondary velocities closer in magnitude to the primary velocity; when the primary velocity is 0 km/s, the secondary velocities are closer to 0 km/s, and the planetary spectrum is closer to the telluric rest-frame. That the near-zero primary velocity approach brings the planetary spectrum closer to the telluric frame when combined with a smaller K p could detract from its advantage over a more random set of primary velocities. Regardless of how these factors work out, Figure 4 confirms that a near-zero primary velocity observing strategy could not have made up for the order of magnitude needed for a strong detection of HD88133b with the orbital phases of the six L band NIRSPEC1.0 epochs obtained and presented in Piskorz et al. 2016.
The simulations considering near-zero primary velocity epochs taken with the upgraded NIRSPEC instrument, shown in the bottom panel of Figure 4, show the most promising chance of detection. There is a peak centered at K p = 40 km/s. A Gaussian fit to the simulated result (dark purple) with no prior information reports a measurement of 22 ± 20 km/s, with a height of 3.2σ. If this result were from real data, and we were able to assign the peak at ∼17 km/s to noise rather than the planetary signature through either fits with simulations (e.g., Buzard et al. 2020)   assumed here, we could expect a more refined fit.

UPSILON ANDROMEDAE B
Piskorz et al. 2017 reported the direct detection of upsilon Andromedae b at a K p of 55 ± 9 km/s using 7 epochs of Keck/NIRSPEC L band data, 3 epochs of K band data covering the left-hand side of the detector, and 3 K band epochs covering the right-hand side of the detector. For this work, we will again just consider the L band epochs. Piskorz et al. 2017Results Piskorz et al. 2017 approximated the orbit of ups And b as circular, and so reported mean anomaly M values, rather than true anomaly f values, because they would be the same in the circular limit. We find, however, that there was a systematic error in the calculation of the secondary velocities that stemmed from a mismatch between the zero points used in Equations 3 and 2. This resulted in a net error comparable to mean anomalies roughly -3.3% offset from their true values. Table 3 lists the mean anomalies used in Piskorz et al. 2017 and the corrected anomalies measured from pericenter, calculated using orbital parameters from Wright et al. 2009. Figure 5 shows how these offsets affect the resulting log likelihood curve from the seven epochs of ups And NIR-SPEC L band data. The originally published log likelihood curve is in red and is reproduced in black dashed. The corrected log likelihood curves are shown in blue (eccentric orbit) and orange (circular orbit). As in the case of HD 88133 b, we can see here that for low eccentricity orbits, there is no benefit to considering an eccentric orbit rather than assuming a circular one.

Ups And Simulations
We were interested in running similar simulations to those run for HD88133 in Section 3.2 to see whether Piskorz et al. 2017

Generation of Simulated Data
We generate ups And simulated data as described in Section 3.2.1. For these simulations, we use a stellar model generated from the PHOENIX stellar model grid (Husser et al. 2013) interpolated to an effective temperature of 6213 K, a metallicity of 0.12, and a surface gravity of 4.25 (Valenti & Fischer 2005). We assume a stellar radius of 1.7053529 R (Gaia Collaboration et al. 2018) and a planetary radius of 1.0 R Jup unless otherwise stated. The simulated data are continuum corrected with a 3 rd order polynomial fit from 2.8 to 4.0 µm in wavenumber space.
We rotationally broaden the stellar model with a stellar rotation rate of 9.62 km/s (Valenti & Fischer 2005) and limb darkening coefficient of 0.29 (Claret 2000). The FWHM of the instrumental kernels of NIRSPEC1.0 and NIRSPEC2.0 are about 12 and 7.3 km/s, respectively, so while rotational broadening makes little difference to data from NIRSPEC1.0, it would have a larger effect on data from the upgraded NIRSPEC instrument. The stellar spectral model used to analyze the ups And L band data in Piskorz et al. 2017 was not from the PHOENIX grid. Instead, they used a model similar to that described in Lockwood et al. 2014. It was generated from a recent version of the LTE line analysis code MOOG (Sneden 1973) and the MARCS grid of stellar at-mospheres (Gustafsson et al. 2008). Notably, individual abundances were set by matching observed lines for elements that were well measured by NIRSPEC. While tests run on both Tau Boo and ups And NIRSPEC1.0 L band data show that stellar models generated this way are able to measure the true stellar velocities at each epoch with higher likelihoods than PHOENIX stellar models, because these models are generated without a continuum, they cannot be used to generate simulated data. Further, because they rely on fits to NIRSPEC1.0 observational data, they could not be used to generate simulated data outside of the NIRSPEC1.0 order wavelengths. Therefore, we are limited to the PHOENIX stellar model.
We use a planetary model from the SCARLET planetary spectral modeling framework (Benneke 2015) without a thermal inversion. The planetary model used in the original work was also from the SCARLET framework, but did include a thermal inversion. We decide to run simulations with a non-inverted planetary model because most recent hot Jupiter atmospheric studies are finding non-inverted thermal structures (e.g., Birkby et al. 2017;Piskorz et al. 2018;Pelletier et al. 2021).
We use orbital positions f from the final column of Table 3, Wright et al. 2009 orbital parameters, and a K p of 55 km/s in Equation 2 to determine the secondary velocities at each epoch. The primary velocities at each epoch are -49.7, -49.4, -48.9, -30.5, -29.6, -25.4, and -39.2 km/s. Gaussian noise is added at the level of the data (total S/N per pixel = 18192).

Planetary Radius Simulations
As for HD88133b, we first run simulations with an increasing planetary radius. Figure 6 shows the results of these simulations with the planetary radius increasing from 1.0 to 4.0 R Jup in increments of 0.5 R Jup . Table 4 report the parameters from Gaussian fits to the log likelihood curves. While Gaussian models can reliably measure a peak centered around the input K p , the R 2 values show that a Gaussian model would not be justified until at least a planetary radius of 3.5 to 4.0 R Jup . While transiting hot Jupiters have been observed with radii approaching 2 R Jup (e.g., KELT-26 b, Rodríguez Martínez et al. 2020), it is improbable that ups And b would have a radius larger than that. These simulations, therefore, do not suggest that, with a reasonable radius, ups And b could be well detected in the 7 original NIRSPEC1.0 L band epochs.

Upgraded NIRSPEC Simulations
The ups And NIRSPEC2.0 simulations are set up the same way as the HD88133 NIRSPEC2.0 simulations with one exception. Because ups And (K = 2.9) is much brighter than HD88133 (K = 6.2), we assume a S/N per pixel per epoch of 9000, or a total S/N per pixel of 23800, across the 7 epochs. At the average S/N of 6530 per pixel in the NIRSPEC1.0 data, we were already well into the regime where structured noise far outweighs white noise, so anything more should make little to no difference to the results. Figure 7 shows a clear peak at the input K p of 55 km/s. It does, coincidentally, fall at the same position as a structured noise peak (in light purple), suggesting that its significance could be overestimated. Any other   Note. -These simulations were all run with an input Kp of 55 km/s. Prior to fitting, these log likelihood results are subtracted by the mean of their values from -150 to 0 km/s. The Gaussian model is initiated with a 55 km/s center and 10 km/s standard deviation. The Gaussian peak height is reported over σ, which is measured as the standard deviation of the log likelihood structure more than 3∆Kp above or below the best-fit Gaussian center, where ∆Kp is the standard deviation of the best-fit Gaussian model.  value of K p would result in a weaker peak that would need to be distinguished, through some mechanism, from the noise peak at ∼ 55 km/s. With the input K p at 55 km/s, a Gaussian model reports a fit at 57 ± 7 km/s with a height of 2.1σ.
This result is encouraging in that it implies that NIR-SPEC2.0 would have allowed a multi-epoch detection of ups And b with the exact seven epochs presented in Piskorz et al. 2017 even with a planetary radius of 1 R Jup .

Near Zero Primary Velocity Simulations
Because of its relatively large systematic velocity of -28.59 km/s (Nidever et al. 2002), ups And never reaches a primary velocity of 0 km/s. The nearest its primary velocity gets to zero is -2 km/s in late January/early February every year. During this time of year, it would be accessible from Keck during the first few hours of the night, setting, in early February, at around 7 UT. This would optimistically allow for an hour and a half on target after telescope focusing on a good night. Because ups And is a very bright source, enough S/N could be achieved to enter the regime where structured noise, rather than random noise, dominates very quickly. PCAbased telluric correction approaches, like those used in Piskorz et al. 2016 andPiskorz et al. 2017, require enough observation time to witness variation in the telluric spectrum. We run the following simulations assuming that the time before ups And sets would be enough to witness telluric variation sufficient to be picked up by PCA or that the data could be well telluric corrected by another approach. Then, these simulations are run with either the NIRSPEC1.0 or NIRSPEC2.0 set up and with the same orbital phases (f ) as were in the original data, but with primary velocities at each epoch of -2 km/s.   ther atmospheric characterization (e.g., Finnerty et al. 2021).

DISCUSSION
In this work, we reevaluated the multi-epoch detections of HD 88133 b (Piskorz et al. 2016) and ups And b (Piskorz et al. 2017), correcting for errors in the estimated orbital positions at the observation times. Unfortunately, we find that the data is insufficient to report planetary detections or measurements of K p in either case. Multi-epoch detections with small data sets have always been a risk; stellar radial velocity measurements are now made with tens, if not hundreds, of epochs.
HD 88133 b and ups And b were two particularly difficult planets to target. Both orbit very large stars, resulting in planet-to-star contrasts lower than those of typical hot Jupiters. At high resolution, the planet/star photospheric contrast provides an upper bound for the spectroscopic information content and thus gives only a first check on how easily detectable planets may be. Predictions of the line-to-continuum, or spectroscopic, contrasts, which can be significantly lower than photospheric contrasts, would provide a more useful guide to direct detection studies; but are highly dependent on the nature of atmospheric chemistry, in particular whether hazes are present, and thermal structure, meaning that model predictions will have large uncertainties. Here we consider photospheric contrasts as a first glimpse into why the spectroscopic detections of HD 88133 b and ups And b may have been so elusive. Measuring photospheric contrast as the mean ratio of planetary flux to stellar flux across the L band, we estimate planet-to-star contrasts for HD 88133 b and ups And b as 2.3 × 10 −4 and 2.5×10 −4 , respectively, with an assumed radius of 1 R Jup for each planet. By comparison, HD 187123 b, studied in Buzard et al. 2020, has an expected L band contrast of 1.4 × 10 −3 and Tau Boötis b, studied in Lockwood et al. 2014, has an expected contrast somewhere between 1.1-1.5×10 −3 , depending on whether or not water is included in its model spectrum (see Pelletier et al. 2021 for a discussion into water on Tau Boo b). Each of these planets has a contrast nearly an order of magnitude larger than do HD 88133 b and ups And b. Our radius simulations show that in each case, a planetary radius of 3 R Jup would have allowed for a strong detection. As planetto-star contrast increases with R 2 pl , a radius of 3 R Jup would increase their contrasts nearly the order of magnitude needed to be comparable to HD 187123 b and Tau Boo b.
We considered whether the upgrade to the NIRSPEC instrument (Martin et al. 2018) or the near-zero primary velocity observing strategy presented by Buzard et al. 2021 would have allowed for enough reduction in structured noise to reveal these low contrast signals. The combination of both would offer a stronger chance of detection in both cases. Near-zero primary velocity epochs obtained with NIRSPEC2.0 would have allowed for K p measurements of HD 88133 b as 22 ± 20 km/s with a height of 3.2σ (input K p = 40 km/s) and of ups And b as 56 ± 8 km/s with a height of 10.8σ (input K p = 55 km/s). Several factors could explain why ups And b could be much more strongly detected under these conditions. It was observed with 7 epochs while HD 88133 b was observed with 6. Additionally, while both stars have large radii, ups And A is not quite as large as HD 88133 A (1.7053529 vs. 2.20 R ). Perhaps most importantly, ups And A has a higher effective temperature than HD 88133 A (6213 vs. 5438 K). The cooler a stellar effective temperature, the more complex its spectrum will be. Therefore, HD 88133 A would have a more complex spectrum that would allow for more correlation between the stellar component of the data and the planetary spectral model used for cross correlation that gives rise to the structured noise in the final log likelihood results. Additionally, ups And A not only has a less complex spectrum, but also one that is rotationally broadened. The stellar rotational broadening would also work to lessen the degree of correlation between the planetary model and the stellar component of the data. A smaller factor could be the total S/N. Ups And is a much brighter system, and so was simulated with a higher total S/N of 238000 compared to 7000 for HD88133b. Both cases are in a regime where structured noise outweighs white noise, though, so this is not likely a major contributor. Collectively, these factors all work to make ups And b more easily detectable than HD 88133 b. Though, even ups And b was below the detection limit with a small number of NIRSPEC1.0 L band epochs.
The upgrade to the NIRSPEC instrument will provide a significant advantage to multi-epoch planetary detection, due to its increases in both resolution and spectral grasp (Finnerty et al. 2021). However, in these difficult cases, it might be worthwhile to consider new observational approaches altogether. Since white noise does not appear to be the limiting factor in these multi-epoch studies, we could consider an observational campaign on a smaller telescope (for example, UKIRT) that could dedicate more nights to this work. Another approach could be to consider an instrument like IGRINS or GMTNIRS which could simultaneously afford both a higher spectral resolution than NIRSPEC and a wider wavelength coverage, both of which are beneficial for planet detectability (Finnerty et al. 2021). Many of these instruments are optimized for shorter wavelengths (< 2.5µm) than we have observed with NIRSPEC. As such, careful work into the optimal observing strategies as well as instrument settings will be crucial in the multi-epoch approach's journey beyond NIRSPEC.
Ultimately, we want to stress the importance of using simulations in multi-epoch work. Simulations are essential for understanding the origin and structure of the expected noise in high resolution data, considering both white noise and any structured noise that may arise in cross correlation space. They can offer realistic estimates for the overall sensitivity of the data beyond expectations from a shot noise limit. As such, simulations can and should be used in many ways, e.g. for planning observations (e.g., Buzard et al. 2021), for identifying and reducing sources of structured noise (Buzard et al. 2020), and for evaluating approaches of atmospheric characterization (e.g., Finnerty et al. 2021).
Despite its challenges, the NIRSPEC multi-epoch approach has been used to characterize planetary atmospheric structure. Piskorz et al. 2018 combined the multi-epoch detection of KELT-2A b with Spitzer secondary eclipse data. They found that the multi-epoch data provided roughly the same constraints on metallicity and carbon-to-oxygen ratio as the secondary eclipse data. Further, while the secondary eclipse data provided a stronger constraint on f , the stellar incident flux which is a rough measure of energy redistribution, the multiepoch data constrained it to low values, with a 50% confidence interval at 1.26. As models with f 1.5 show a temperature inversion, this indicates that using NIRSPEC1.0 multi-epoch data alone, Piskorz et al. 2018 were able to determine that KELT-2A b has a noninverted thermal structure in the regions probed by ∼3 µm data. Finnerty et al. 2021 used NIRSPEC2.0 L band simulations to look more deeply into the atmospheric constraints that could be made with multi-epoch data and found that warm Jupiters' (T eq ∼ 900 K) carbonto-oxygen ratios could be constrained enough to differentiate between substellar, stellar, and superstellar values. While planetary detection using the multi-epoch approach can a challenging pursuit, once the true planetary peak has been identified, the approach holds potential for detailed atmospheric characterization.

CONCLUSION
In this work, we present and correct errors in the multiepoch detections of HD 88133 b (Piskorz et al. 2016) and ups And b (Piskorz et al. 2017). Unfortunately, we find that the original NIRSPEC1.0 L band data presented (6 epochs for HD 88133 b, 7 for ups And b) are insufficient for planetary detections. We run simulations to determine what would have been required for confident detections. Ups And b could have been strongly detected (10.8σ) if its seven L band epochs had been taken with the upgraded NIRSPEC instrument and following the near-zero primary velocity observing strategy presented by Buzard et al. 2021. HD 88133 b would be more difficult to detect, because of its larger stellar radius and lower stellar effective temperature, and would likely have required more, carefully planned, data. This has been a very difficult result to discover and to publish. C.B. wants to acknowledge and thank D.P. for her mentorship. She is an incredible scientist, researcher, and mentor, and his journey through graduate school would not have been possible without her support. As Judith Butler wrote, "There was and remains warrant ... to distinguish between self-criticism that promises a more democratic and inclusive life for the movement and criti-cism that seeks to undermine it altogether." We hope this publication falls into the former category and promises a more democratic and inclusive life for the multi-epoch cross-correlation approach moving forward.
We want to acknowledge Katie Kaufman's contribution to the ups And NIRSPEC data reduction along with support from NASA XRP grant NNX16AI14G. We thank our anonymous reviewer for their kind and close reading of our manuscript.