# Evidence for the Direct Detection of the Thermal Spectrum of the Non-Transiting Hot Gas Giant HD 88133 b

## Abstract

We target the thermal emission spectrum of the non-transiting gas giant HD 88133 b with high-resolution near-infrared spectroscopy, by treating the planet and its host star as a spectroscopic binary. For sufficiently deep summed flux observations of the star and planet across multiple epochs, it is possible to resolve the signal of the hot gas giant’s atmosphere compared to the brighter stellar spectrum, at a level consistent with the aggregate shot noise of the full data set. To do this, we first perform a principal component analysis to remove the contribution of the Earth’s atmosphere to the observed spectra. Then, we use a cross-correlation analysis to tease out the spectra of the host star and HD 88133 b to determine its orbit and identify key sources of atmospheric opacity. In total, six epochs of Keck NIRSPEC L band observations and three epochs of Keck NIRSPEC K band observations of the HD 88133 system were obtained. Based on an analysis of the maximum likelihood curves calculated from the multi-epoch cross correlation of the full data set with two atmospheric models, we report the direct detection of the emission spectrum of the non-transiting exoplanet HD 88133 b and measure a radial projection of the Keplerian orbital velocity of 40 15 km/s, a true mass of 1.02, a nearly face-on orbital inclination of 15, and an atmosphere opacity structure at high dispersion dominated by water vapor. This, combined with eleven years of radial velocity measurements of the system, provides the most up-to-date ephemeris for HD 88133.

###### Subject headings:

techniques: spectroscopic — planets and satellites: atmospheres^{16}

## 1. Introduction

Since the discovery of 51 Peg b (Mayor & Queloz, 1995), the radial velocity (RV) technique has proven indispensable for exoplanet discovery. Hundreds of exoplanets have been revealed by measuring the Doppler wobble of the exoplanet host star (Wright et al., 2012), principally at visible wavelengths. To first order, the RV method yields the period and the minimum mass () of the orbiting planet. In order to complete the characterization of a given exoplanet, one would want to measure its radius and constrain its atmospheric constituents. Traditionally, this information is accessible only if the planet transits its host star with respect to our line of sight via transmission or secondary eclipse photometry. Successes with these techniques have resulted in the detections of water, carbon monoxide, and methane on the hottest transiting gas giants (Madhusudhan et al., 2014). These gas giants orbit their host stars in days, are known as hot Jupiters, and have an occurence rate of only 1 in the exoplanet population (Wright et al., 2012). Broadband spectroscopic measurements of transiting hot Jupiter atmospheres are rarely able to resolve molecular bands, let alone individual lines, creating degeneracies in the solutions for atmospheric molecular abundances.

High-resolution infrared spectroscopy has recently provided another route to the study of exoplanet atmospheres, one applicable to transiting and non-transiting planets alike. Such studies capitalize on the Doppler shift between the stellar and planet lines, allowing them to determine the atmosphere compositions, true masses, and inclinations of various systems. With the CRyogenic InfraRed Echelle Spectrograph (CRIRES) at the VLT, Snellen et al. (2010) provided a proof-of-concept of this technique and detected carbon monoxide in the atmosphere of the transiting exoplanet HD 209458 b consistent with previous detections using Hubble Space Telescope data (Swain et al., 2009). By detecting the radial velocity variation of a planet’s atmospheric lines in about six hours of observations on single nights, further work has detected the dayside and nightside thermal spectra of various transiting and non-transiting hot Jupiters, reporting detections of water and carbon monoxide, as well as the presence of winds and measurements of the length of day (Snellen et al., 2010; Rodler et al., 2012; Snellen et al., 2014; Brogi et al., 2012, 2013, 2014, 2016; Birkby et al., 2013; de Kok et al., 2013; Schwarz et al., 2015). With HARPS, Martins et al. (2015) recently observed the reflected light spectrum of 51 Peg b in a similar manner, combining 12.5 hours of data taken over seven nights when the full dayside of the planet was observable.

Lockwood et al. (2014) studied the hot Jupiter tau Boo b using Keck NIRSPEC (Near InfraRed SPECtrometer), confirmed the CRIRES measurement of the planet’s Keplerian orbital velocity, and detected water vapor in the atmosphere of a non-transiting exoplanet for the first time. NIRSPEC was used to observe an exoplanet’s emission spectrum over 2-3 hours each night across multiple epochs, in order to capture snapshots of the planet’s line-of-sight motion at distinct orbital phases. In combination with the many orders of data provided by NIRSPEC’s cross dispersed echelle format and the multitude of hot Jupiter emission lines in the infrared, Lockwood et al. (2014) achieved sufficient signal-to-noise to reveal the orbital properties of tau Boo b via the Doppler shifting of water vapor lines in its atmosphere.

Here, we continue our Keck NIRSPEC direct detection program with a study of the emission spectrum of the hot gas giant HD 88133 b, a system that allows us to test the brightness limits of this method and develop a more robust orbital dynamics model that can be applied to eccentric systems. In Section 2 we present new (stellar) radial velocity (RV) observations of HD 88133 and an updated ephemeris. In Section 3 we outline our NIRSPEC observations, reduction, and telluric correction method. In Section 4 we describe the cross correlation and maximum likelihood analyses, and present the detection of the thermal spectrum of HD 88133 b. In Section 5 we discuss the implications of this result for the planet’s atmosphere and for future observations.

## 2. HIRES Observations and RV Analysis

The RV measurements of HD 88133 have been made under the purview of the California Planet Survey (CPS; Howard et al. 2010) with the HIgh Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) at the W.M. Keck Observatory. Seventeen RV measurements of HD 88133 were published in an earlier study (Fischer et al., 2005), and here we extend that data set to 55 individual RV measurements, having a baseline of eleven years (see Table 1). RV data are reduced with the standard CPS HIRES configuration and reduction pipeline (Wright et al., 2004; Howard et al., 2010; Johnson et al., 2010). Doppler shifts are recovered by comparison to an iodine absorption spectrum and a modeling procedure presented in Butler et al. (1996) and Howard et al. (2011). Processed RV data are shown in Figure 1.

Julian Date | Radial Velocity | |
---|---|---|

(- 2,440,000) | (m/s) | (m/s) |

13014.947812 | -21.97 | 2.03 |

13015.947488 | 23.44 | 2.06 |

13016.952546 | 20.55 | 1.91 |

13044.088461 | 21.71 | 1.63 |

13044.869410 | -24.07 | 1.46 |

13045.843414 | -31.17 | 1.42 |

13046.081308 | -19.97 | 1.52 |

^{17}

The RV data are fit with a Markov Chain Monte Carlo technique following Bryan et al. (2016). Eight free parameters (six orbital parameters: the velocity semi-amplitude , the period of the orbit , the eccentricity of the orbit , the argument of periastron , the true anomaly of the planet at a given time , and the arbitrary RV zero point ; a linear velocity trend ; and a stellar jitter term as in Isaacson & Fischer 2010) having uniform priors contribute to the model and are simultaneously fit to the data . We initiate the MCMC chains at the values published by Fischer et al. (2005), allowing the chains to converge quickly and avoiding degeneracies in and . The likelihood function is given by

(1) |

where is the instrument error and is the stellar jitter. The stellar jitter term is added in quadrature to the uncertainty value of each RV measurement. Best-fit orbital elements indicated by this analysis, as well as other relevant system parameters, are included in Table 2, and the best-fit velocity curve is shown in the Figure 1. This represents a substantial improvement to the ephemeris originally published in Fischer et al. (2005). We combine these values with the velocities derived by our NIRSPEC analysis described in Sections 3 and 4 to break the degeneracy for HD 88133 b.

Property | Value | Reference |
---|---|---|

Stellar | ||

Mass, | 1.18 0.06 | (1) |

Radius, | 1.943 0.064 | (2) |

Effective temperaure, | 5438 34 K | (1) |

Metallicity, | 0.330 0.05 | (1) |

Surface gravity, | 3.94 0.11 | (1) |

Rotational velocity, | 2.2 0.5 | (1) |

Systemic velocity, | -3.45 0.119 km/s | (3) |

K band magnitude, | 6.2 | (4) |

Velocity semi-amplitude, | 32.9 m/s | (5) |

RV zero point, | 3.08 m/s | (5) |

Velocity trend, | -0.0013 m/s/yr | (5) |

Stellar jitter, | 4.68 | (5) |

Planetary | ||

Indicative mass, | 0.27 | (5) |

Mass, | 1.02 | (6) |

Inclination, | 15 | (6) |

Semi-major axis, | 0.04691 0.0008 AU | (7) |

Period, | 3.4148674 | (5) |

Eccentricity, | 0.05 | (5) |

Argument of periastron, | 7.22 | (5) |

Time of periastron, | 2454641.984 JD | (5) |

Phase uncertainty, | 6.34 | (5) |

References. – (1) Mortier et al. (2013); (2) Torres et al. (2010); (3) Chubak et al. (2012) (4) Wenger et al. (2000) (5) From HIRES measurements presented in Section 2; (6) From NIRSPEC measurements presented in Sections 3 and 4; (7) Butler et al. (2006)

## 3. NIRSPEC Observations and Data Reduction

We pursue multi-epoch observations having planetary features shifted with respect to the star’s spectrum and the Earth’s atmosphere in order to develop techniques that can eventually be used on more slowly moving planets nearer the habitable zone. For nearly synchronously rotating hot Jupiters, near-infrared emission from the dayside is likely to be most readily detectable. This strategy is fundamentally different from that used at CRIRES to date since we do not allow the planet’s signal to move across several pixels on the detector during one night of observations.

### 3.1. Observations

Data were taken on six nights (2012 April 1 and 3, 2013 March 10 and 29, 2014 May 14, 2015 April 8) in L band and three nights (2015 November 21, 2015 December 1, and 2016 April 15) in K band in an ABBA nodding pattern with the NIRSPEC instrument at the W.M. Keck Observatory (McLean et al., 1998). With a 0.4x24 slit NIRSPEC has an L band resolution of 25,000 (30,000 in K band). Individual echelle orders cover from 3.4038-3.4565/3.2567-3.3069/3.1216-3.1698/2.997-3.044 m in the L band and from 2.3447-2.3813/2.2743-2.3096/2.2085-2.2422/2.1464-2.1788/2.0875-2.1188/2.0319-2.0619 m in the K band.

A schematic of the planet’s orbit during our observations is given in Figure 2 assuming the best-fit orbital parameters from our HIRES RV analysis in Section 2. Details of our observations are given in Table 3.

Date | Julian Date | Mean anomaly | True anomaly | Barycentric velocity | Integration time | S/N |
---|---|---|---|---|---|---|

(- 2,440,000 days) | (2 rad) | (2 rad) | (km/s) | (min) | ||

L band (3.0 - 3.4 m) | ||||||

2012 April 1 | 16018.837 | 0.89 | 0.86 | -20.96 | 140 | 1680 |

2012 April 3 | 16020.840 | 0.48 | 0.48 | -21.66 | 140 | 2219 |

2013 March 10 | 16361.786 | 0.29 | 0.33 | -11.59 | 180 | 2472 |

2013 March 29 | 16380.726 | 0.84 | 0.80 | -19.70 | 150 | 1812 |

2014 May 14 | 16791.796 | 0.18 | 0.22 | -29.27 | 120 | 1694 |

2015 April 8 | 17120.835 | 0.50 | 0.50 | -23.01 | 160 | 2938 |

K band (2.0 - 2.4 m) | ||||||

2015 November 21 | 17348.129 | 0.06 | 0.68 | 29.95 | 60 | 2701 |

2015 December 1 | 17358.117 | 0.96 | 0.62 | 29.25 | 60 | 2823 |

2016 April 15 | 17166.300 | 0.54 | 0.53 | -29.15 | 80 | 2466 |

### 3.2. Extraction of 1-D Spectra

We flat field and dark subtract the data using a Python pipeline à la Boogert et al. (2002). We extract one-dimensional spectra and remove bad pixels, and calculate a fourth-order polynomial continuum by fitting the data to a model telluric spectrum after the optimal source extraction. For L band data, the wavelength solution (described to the fourth order as , is pixel number and , , , and are free parameters) is calculated by fitting the data to a model telluric spectrum. However, since telluric lines are generally weaker near 2 m, the wavelength solution for K band data is calculated by fitting the data to a combination of model telluric and stellar spectra (given that the stellar relative velocity is well known from optical data, and is later confirmed by the cross-correlation analysis described in Section 4.3). We use an adapted PHOENIX model for our model stellar spectrum in the K band, as described in Section 4.1 (Husser et al., 2013). We show one order of reduced L band spectra in the top panel of Figure 3. We fit an instrument profile to the data and save it so that we may apply it to our stellar and planetary models. This instrument profile is similar to the formulation given in Valenti et al. (1995) and is parameterized as a central Gaussian with four left and four right satellite Gaussians, all with variable widths. Often, the best-fit widths of the third and fourth left and right satellite Gaussians are zero.

### 3.3. Telluric Correction with Principal Component Analysis

In this work, we depart from our traditional methods of division by an atmospheric standard (typically an A star, c.f. Boogert et al. 2002) and/or line-by-line telluric correction (modeling atmospheric abundances with the TERRASPEC software; Bender et al., 2012) in favor of a more automated and repeatable technique: principal component analysis (PCA). This efficient method of telluric correction was also implemented by de Kok et al. (2013) in their reduction of CRIRES data on HD 189733 b. PCA rewrites a data set in terms of its principal components such that the variance of the data set with respect to its mean or with respect to a model is reduced. The first principal component encapsulates the most variance; the second, the second most, etc. Over the course of an observation, the telluric components should vary the most as the air mass and atmospheric abundances change, and the planet lines should remain constant. Note that we observe for only 2-3 hours at a time and at a lower resolution than CRIRES, so the planet lines do not smear. For a typical Keplerian orbital velocity of a hot Jupiter, we would have to observe for at least four hours to smear the planet lines across the NIRSPEC detector pixels. To remove the telluric lines and any other time-varying effects, we aim to isolate only the strongest principal components.

To perform the principal component analysis^{18}

(2) |

where is the row number, is the column number, is either the mean of or a telluric model, and is the standard deviation of the values in column . We guide our principal component analysis with a telluric model for (rather than the mean of ) that uses baseline values for water vapor^{19}^{20}^{21}

(3) |

A singular value decomposition of the covariance matrix is then performed to find the principal components:

(4) |

where contains the left singular vectors (or the eigenvectors), is a diagonal matrix of the singular values (or the eigenvalues), and contains the right singular vectors. The first three eigenvectors, or principal components, for a 3.12-3.17 m order taken on 2013 March 29 are shown in Figure 3. The first component recovers changes in the total systemic signal with air mass; the second encapsulates changes in abundances of telluric species, resulting in adjustments to line cores and wings; and the third describes changes to the plate scale. Higher order components typically reflect instrumental fringing and other small effects.

We reconstruct the time-varying portion of each AB spectrum by combining the first principal components of the data set, given by , where is the first columns of and is the x upper-left part of . Rank- data, , can be built as

(5) |

To produce a telluric-corrected spectrum, , each row of is divided by its corresponding un-mean normalized row in :

(6) |

Finally, we collapse the rows of data in and clip regions of substantial telluric absorption (75%). Depending on the order, anywhere between 30 and 60% of the data is removed by clipping out strong features. This results in a single high signal-to-noise spectrum for each night of observations. The final telluric-corrected spectrum for 2013 March 29 is given in the final panel of Figure 3.

The version of PCA described here diverges from the approach outlined in de Kok et al. (2013) which used PCA to determine the eigenvectors making up the light curves in each spectral channel. Our formulation calculates the eigenvectors comprising each observed spectrum. We also guide our principle component analysis with a telluric model. The equivalent for de Kok et al. (2013) would have been to guide the PCA with vectors for air mass, water vapor content, etc.

For data taken in L band, we find that PCA can reliably correct for the Earth’s atmosphere for all orders of data. However, for K band data, we cannot effectively remove the dense, optically-thick telluric forest of CO lines in the order spanning from 2.03 - 2.06 m, and we omit this wavelength range from subsequent analysis.

It is essential to determine the efficacy of PCA in removing telluric signatures and further ensure that we are not removing the planet’s signal as well. Since 99.9% of the variance is explained by the first principal component, we find that the following results are roughly consistent when different numbers of principal components are removed from the data. We calculate the percent variance removed by each component and, if we assume the planet signal is on the order of 10 of the total signal, then for most cases we would have to remove about fifteen components to delete the planet signal. We have experimented with incrementally removing up to ten principal components as input to the subsequent analysis; the results and conclusions presented below use data with the first five principal components removed.

## 4. NIRSPEC Data Analysis and Results

We run a two-dimensional cross-correlation analysis (TODCOR algorithm; Zucker & Mazeh 1994) to find the optimum shifts for the stellar and planetary spectra entwined in our telluric- and instrument-corrected data. This requires accurate model stellar and planet spectra.

### 4.1. Model Stellar Spectrum

Our synthetic stellar spectrum is a PHOENIX model (Husser et al., 2013) interpolated to match the published effective temperature , surface gravity , and metallicity of HD 88133 listed in Table 2. As HD 88133 has a 5 km/s, instrumental broadening dominates, and we convolve the stellar model with the instrumental profile calculated in Section 3.2.

### 4.2. Model Planetary Spectrum

We have computed the high-resolution thermal emission spectrum of HD 88133 b using both the SCARLET (Benneke, 2015) and PHOENIX (Barman et al., 2001, 2005) frameworks. An example of one order of our L band planet models is shown in Figure 4. Both models compute the thermal structure and equilibrium chemistry of HD 88133 b given the irradiation provided by the host star. Models are computed for a cloud-free atmosphere with solar elemental composition (Asplund et al., 2009) at a resolving power of R 250K, and assume perfect heat redistribution between the day and night sides. The model spectra are subsequently convolved with the instrumental profile derived in Section 3.2 (Figure 4). We find consistent results for both models despite minor differences in the molecular line lists used.

The SCARLET model considers the molecular opacities of , , , HCN, , and and from the high-temperature ExoMol database (Tennyson & Yurchenko, 2012), and , , , , , , , and from the HITRAN database (Rothman et al., 2009). Absorption by the alkali metals (Li, Na, K, Rb, and Cs) is modeled based on the line strengths provided in the VALD database (Piskunov et al., 1995) and -broadening prescription provided in Burrows & Volobuyev (2003). Collision-induced broadening from and collisions is computed following Borysow (2002).

Unlike SCARLET, PHOENIX is a forward modeling code that converges to a solution based on traditional model atmosphere constraints (hydrostatic, chemical, radiative-convective, and local thermodynamic equilibrium) for an assumed elemental composition. The PHOENIX model uses similar opacities as SCARLET (for example, most of the latest linelists from ExoMol and HITRAN). Additional line data for metal-hydrides come from Dulick et al. (2003) (and references therein). Broadening of alkali lines follows Allard et al. (2003). These differences are of only minor importance for this study because SCARLET and PHOENIX use the same water linelist and water opacity dominates the spectral features across the spectral range of our observations (Figure 4).

Based on the effective temperatures of the planet and the star, the photometric contrast (defined as the ratio of the planet flux to the stellar flux) is on the order of 10. This is also a rough upper bound for the spectroscopic contrast . Since the cross-correlation analyses described in Section 4.3 are not very sensitive to contrast ratios, varying the value of does not change our conclusions on the radial velocity of the planet, and thus the system inclination (Lockwood et al., 2014). However, the specific value of does affect our conclusions on the composition and structure of the planet’s atmosphere. That is, the overall velocity structure of the cross-correlation surface is not much affected by , though the size and structure of the final maximum likelihood peak near the planet’s signature will be.

### 4.3. Two-Dimensional Cross Correlation

Each order of data for each epoch is cross-correlated with the model predictions to determine the cross-correlation function (CCF) using the TODCOR algorithm (Zucker & Mazeh, 1994). This calculation results in a two-dimensional array of correlation values for shifts in the velocities of the star and planet.

In testing our models, we find that the correlation coefficient between our stellar and planet models for the two orders spanning from 2.31 - 2.38 m is generally an order of magnitude higher than the same correlation coefficient in any other order in the L or K band, and so we omit them from this study. This behavior is due to the strong correlation between stellar and planetary CO at R=30,000. HD 88133 has an effective temperature of 5438 K and its CO band (and especially the CO band head at 2.295 m) is extremely prominent. The K band data analysis that follows is only performed on the three orders spanning from 2.10 - 2.20 m. The main absorbers in this region include carbon dioxide and water vapor.

For each night’s observation, we combine the CCF’s of each order with equal weighting to produce the maximum likelihood curves shown in Figure 5. Lockwood et al. (2014) showed that the cross-correlation function is proportional to the log of the likelihood. At each epoch, we can easily retrieve and confirm the stellar velocity, as shown by the single strong peak in panel A of Figure 5. The stellar velocity is dominated by the barycentric velocity at the time of observation and the systemic velocity of HD 88133. We are insensitive to the reflex motion of the star, which is on the order of 0.01 km/s.

The retrieval of the planet velocity is more complex, as evidenced by the multiple peaks in panels B-J of Figure 5, and requires combining the data from multiple epochs. Though there are many peaks and troughs in the maximum likelihood curves produced for each night’s observation (Figure 5), only one peak per night represents the properly registered correlation of the data with the model planetary spectrum. This is not to say that most of the maximum likelihood peaks in Figure 5 are spurious; rather, they are the results of unintended systematic structure in the cross-correlation space.

### 4.4. Planet Mass and Orbital Solution

The orbit of HD 88133 b is slightly eccentric, and we calculate the velocity of the planet for a given epoch as a function of its true anomaly according to

(7) |

where is the planet’s Keplerian orbital velocity, is the longitude of periastron measured from the ascending node, is the eccentricity of the orbit, and is the combined systemic and barycentric velocities at the time of the observation. We test a range of orbital velocities from -150 to 150 km/s in steps of 1 km/s, and in turn test a range of planet masses and orbital inclinations. Figure 6 shows the maximum log likelihood versus the planet’s Keplerian orbital velocity.

When the six maximum likelihood curves produced from L band data (panels B-G in Figure 5) are combined with equal weighting according to Eq. 7, we see that the most likely value for the radial projection of the planet’s Keplerian orbital velocity is 41 16 km/s. These error bars are calculated as in Lockwood et al. (2014) by fitting a Gaussian to the likelihood peak and reporting the value of , assuming all the points on the maximum likelihood curve are equally weighted. We deduce that the peak in the likelihood curve based on the PHOENIX model at 60 km/s does not represent the planet’s velocity, and we prove this in Section 5.1. The same calculation applied to the three nights of K band data (panels H-J in Fig. 5) suggests that is 32 12 km/s.

The combination of all nine nights of data yields = 40 15 km/s. It is this value of that we use to calculate the secondary velocity curve shown in the bottom panel of Figure 2. The peak near = 40 km/s is consistent between the two planet models cross-correlated with the data. Here, given the full suite of data, we calculate error bars on the individual points with jackknife sampling. One night’s worth of data is removed from the sample and the maximum likelihood calculation is repeated. The standard deviation of each point amongst the resulting eight maximum likelihood curves is proportional to the error bars on each point on the maximum likelihood curve.

We note that error bars calculated by jackknife sampling and shown in Figure 6 are merely an estimate. In fact, for the Gaussian fit, the reduced chi-squared value (chi-squared divided by the number of degrees of freedom) is 0.1, indicating that the error bars are overestimated. This can be explained by the fact that there is high variance between jackknife samples, driving a high standard deviation and therefor large error bars. To examine this behavior further, we fit a Gaussian distribution (indicating the presence of a planetary signal) and a flat line (indicating no planetary signal) and determine the significance of the signal. As in Kass & Raftery (1995), we define the Bayes factor to be the ratio of likelihoods between two models, in this case the likelihood of the Gaussian distribution compared to the likelihood of the straight line. must be greater than 10 for a model to be very strongly preferred.

For the Gaussian distribution compared the the straight line, is nearly 430, indicating the the Gaussian approximation to the signal at 40 km/s is significantly stronger than a flat line. Even if our error bars are overestimated, they are likely not overestimated by a factor of 100. Combining with the parameters given in Table 2, a of 40 15 km/s implies that the true mass of HD 88133 b is 1.02 and its orbital inclination is 15.

Note that the values of implied by the most likely value of often, but do not always, correspond with peaks in the maximum likelihood curves for each night, as indicated by the vertical grey dashed lines in Figure 5. Especially for nights having a small line-of-sight velocity, planetary lines may be lost in the telluric and/or stellar cross correlation residuals.

## 5. Discussion

### 5.1. Tests of the Orbital Solution

We first check our detection of HD 88133 b’s emission spectrum at a radial velocity projection of 36 km/s by varying the spectroscopic contrast uniformly with orbital phase. We tested nine values of from 10 to 10, and find the maximum likelihood peak near 40 km/s is robust for 10.

We create a “shuffled” planetary model by randomly rearranging chunks of each planetary atmosphere model. If the maximum likelihood peak at 36 km/s is real, then cross correlating our data with a “shuffled” planetary model (which has no coherent planet information) should show little to no peak at the expected . And, indeed, the data-“shuffled” planetary model cross-correlation shows no peak at 40 km/s while the peak at 60 km/s remains for the PHOENIX model (see Panel D of Figure 6).

We also check our results by varying the orbital elements of the system. We obtain roughly the same values for Keplerian orbital velocity (within the error bars) for various combinations of eccentricities down to 0 and arguments of perihelion within 20 of the reported value. Our results are most sensitive to these orbital elements as they affect the calculations of the true anomaly and secondary velocity, and therefore the positions of the dashed vertical lines in each epoch’s maximum likelihood curve shown in Figure 5. Even with a different ephemeris, so long as the dotted vertical lines are near their current peaks, we obtain a comparable final result for the Keplerian orbital velocity of the system.

We note that as HD 88133 b’s orbit is slightly eccentric, it is not truly tidally locked. Calculations following Hut (1981) suggest that the planet spins about 10% faster than synchronous. Our strategy prefers that the planet be tidally locked so that as much of the planet’s (dayside) emission as possible is captured when the planet has a high line-of-sight velocity relative to the star.

Finally, we have reprocessed the tau Boo b data published by Lockwood et al. (2014) with the methods presented in Section 3. The specific departure from Lockwood et al. (2014) is the use of principal component analysis to correct for telluric features in the data (Section 3.3). Our analysis of five nights of L band data recovers a projected Keplerian orbital velocity of 121 8 km/s, a mass of 5.39, and an orbital inclination of 50 for tau Boo b. This is in good agreement with Lockwood et al. (2014) (=1115 km/s) as well as Brogi et al. (2012) (=1103.2 km/s) and Rodler et al. (2012) (=11511 km/s), thereby validating our PCA-based methods.

### 5.2. Observation Notes

The peak in log likelihood produced solely from L band data (panel A of Figure 6) is an order of magnitude larger than the peak produced solely from K band data, though the values of preferred by each data set are consistent. Though the single channel signal-to-noise ratio for our data is greater in the K band than in the L band, we only observe HD 88133 in K for three nights and only use three orders of data in our final cross-correlation analysis whereas we observe in L for six nights and use all four orders of data. Furthermore, six nights of L band observations on Keck yields an aggregate shot noise of 800,000 for all epochs and all wavelength bins, suggesting that the detection of a planet having a spectroscopic contrast down to 10 is feasible.

Additionally, we attempt to detect the planet’s CO band near 2.295 m, and find that effectively separating the stellar spectrum from the planet’s is a complex process. Future observations should consider avoiding the CO band head and focusing on the CO comb (low- to moderate- angular momentum P and R branches) itself, particularly the regions between the stellar lines where shifted planetary CO should be present. These intermediate regions have v 60 km/s, which is certainly sufficient for the detection of shifted planetary CO and especially so given the high resolution of the next generation of cross dispersed infrared echelle spectrographs, including iSHELL and the upgraded CRIRES and NIRSPEC instruments.

### 5.3. The Spectrum of HD 88133

We took the opportunity during our reprocessing of the tau Boo b data to evaluate the required accuracy of the stellar model. The original Lockwood et al. (2014) analysis was performed with a stellar spectrum using a MARCS solar model with adjustments made to specific line parameters. We re-run the analysis with a PHOENIX model for tau Boo. The correct stellar velocity at each observation epoch is still recovered and the final result for the planet’s remains unchanged. This suggests that, for the sake of detecting the planet’s spectrum and thus the planet’s velocity, a detailed model of the star is not necessarily required; however a refined stellar spectrum will be critical for learning about the planet’s atmosphere in detail.

### 5.4. The Atmosphere of HD 88133 b

Acquisition of both L and K band NIRSPEC data provides a unique opportunity to constrain the atmosphere of HD 88133 b. Generally, the spectroscopic contrast is dominated by water vapor in the L band and by carbon monoxide and other species in the K band. Ultimately, the shape of the log likelihood peak in Figure 6 will provide information about the structure of the atmosphere (e.g., Snellen et al. 2010 and Brogi et al. 2016) and even information about the planet’s rotation rate (e.g., Snellen et al. 2014 and Brogi et al. 2016). We do not presently consider orbital phase variations in atmosphere dynamics, radiative transfer, and the resulting spectral signatures from the day to night sides of the hot Jupiter. We plan to examine detailed properties of HD 88133 b’s atmosphere in a future study.

## 6. Conclusion

We report the detection of the emission spectrum of the non-transiting exoplanet HD 88133 b using high resolution near-infrared spectroscopy. This detection is based on the combined effect of thousands of narrow absorption lines, predominantly water vapor, in the planet’s spectrum. We find that HD 88133 b has a Keplerian orbital velocity of 40 15 km/s, a true mass of 1.02, and a nearly face-on orbital inclination of 15.

Direct detection of hot Jupiter atmospheres via this approach is limited in that it cannot measure the absolute strengths of molecular lines, relative to the photometric contrast. Thus, this method will yield degeneracies between the vertical atmospheric temperature gradients and absolution molecular abundance ratios, but the relative abundances of species should be better constrained. For transiting planets having Spitzer data, it should be possible to better measure absolute abundances by comparing Spitzer eclipse depths and the output of our cross-correlation analyses using various planetary atmosphere models.

With the further refinement of this technique and with the improved future implementation of next-generation spectrometers and coronagraphs, especially on the largest optical/infrared telescopes, we are optimistic that this method may be extended to the characterization of terrestrial atmospheres at Earth-like semi-major axes. This paper shows progress in that direction by presenting an algorithm capable of removing telluric lines whilst preserving planet lines even if the planet does not move significantly during the observations.

### Footnotes

- affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
- affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
- affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
- affiliation: King Abdullah University of Science and Technology, Thuwal Saudi Arabia
- affiliation: Lunar and Planetary Laboratory, University of Arizona, Tuscon, AZ 85721
- affiliation: Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802
- affiliation: Center for Exoplanets and Habitable Worlds, Pennsylvania State University, University Park, PA 16802
- affiliation: Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125
- affiliation: Naval Research Laboratory, Washington, DC 20375
- affiliation: Department of Astronomy, Yale University, 260 Whitney Avenue, New Haven, CT 0651
- affiliation: Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125
- affiliation: Department of Astronomy, University of California Berkeley, Berkeley, CA 94720
- affiliation: Harvard-Smithsonian Center for Astrophysics; Institute for Theory and Computation, Cambridge, MA 02138
- affiliationtext: Institute for Astronomy, University of Hawaii at Manoa, Honolulu, HI 96822
- The full set of RVs for this system is available as an electronic table online.
- PCA tutorial, http://stats.stackexchange.com/questions/

134282/relationship-between-svd-and-pca-how-to-use-svd-to-

perform-pca - Caltech Submillimeter Observatory, http://cso.caltech.edu/tau/
- Earth System Research Laboratory at Mauna Loa, ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_weekly_mlo.txt
- Earth System Research Laboratory at Mauna Loa, ftp://ftp.cmdl.noaa.gov/data/greenhouse_gases/ch4/in-situ/surface/mlo/ch4_mlo_surface-insitu_1_ccgg_DailyData.txt

### References

- Allard, N. F.,Allard, F., Hauschildt, P. H., Kielkopf, J. F., & Machin, L. 2003, A&A, 411, L473
- Asplund, M., Grevesse, N., Suval, J.A., Scott, P. 2009, ARAA, 47, 1
- Barman, T. S., Hauschildt, P. H., & Allard, F. 2001, ApJ, 556, 885
- Barman, T. S., Hauschildt, P. H., & Allard, F. 2005, ApJ, 632, 1132
- Benneke, B. 2015, arXiv:1504.07655
- Birkby, J. L., de Kok, R. J., & Brogi, M. et al. 2013, MNRAS, 436, L35
- Boogert, A. C. A., Blake, G. A., & Tielens, A. G. G. M. 2002, ApJ, 577, 271
- Borysow, A. 2002, A&A, 390, 4
- Brogi, M., Snellen, I.A.G., de Kok, R.J., et al. 2012, Nature, 486, 502
- Brogi, M., Snellen, I.A.G., de Kok, R.J., et al. 2013, ApJ, 767, 1
- Brogi M., de Kok R. J., Birkby J. L., Schwarz H. & Snellen I. A. G. 2014, A&A, 565A, 124B
- Brogi, M., de Kok, R.J., Albrecht, S. et al. 2016, ApJ, 817, 2
- Bryan, M.L., Knutson, H.A., Howard, A.W., et al. 2016, ApJ, 821, 2
- Burrows, A., & Volobuyev, M. 2003, ApJ, 583, 985
- Butler R. P., Marcy G. W., Williams E. et al 1996 PASP 108 500
- Butler, R.P., Wright, J.T., Marcy, G.W., et al. 2006, ApJ, 646, 505
- Chubak C., Marcy G., Fischer D. A. et al. 2012, arXiv:1207.6212
- de Kok R. J., Brogi M., Snellen I. A. G. et al. 2013, A&A, 554, A82
- Dulick, M., Bauschlicher, C. W., Jr., Burrows, A., et al. 2003, ApJ, 594, 1
- Fischer, D.W., Laughlin, G., Butler, P., et al. 2005, ApJ, 620, 481.
- Howard A. W., Johnson J. A., Marcy G. W. et al. 2010, ApJ, 721, 1467
- Howard A. W., Johnson J. A., Marcy G. W. et al. 2011, ApJ, 730, 10
- Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
- Hut, P. 1981, A&A, 99, 1
- Isaacson H. & Fischer D. 2010, ApJ, 725, 875
- Johnson J. A., Howard A. W., Marcy G. W. et al. 2010, PASP, 122, 888
- Kass, R.E., & Raftery, A.E. 1995, JASA, 90,773
- Lockwood, A.C., Johnson, J.A., Bender, C.F., et al. 2014, ApJ, 783, L29
- Madhusudhan, N., Knutson, H.A., Fortney, J.J., & Barman, T. 2012, Exoplanetary atmospheres. In Protostars and Planets VI (H. Beuther et al., eds.), pp. 739-762. Univ. of Arizona, Tuscon.
- Martins, J.H.C., Santos, N.C., Figueira, P., et al. 2015, A&A, 576, A134
- Mayor, M. & Queloz, D. 1995, Nature, 378, 23
- McLean, I. S., Becklin, E. E., Bendiksen, O., et al. 1998, Proc. SPIE 3354, 566
- Mortier, A., Santos, N.C., Souse, S.G., et al. 2013, A&A, 557, A70
- Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&A S, 112, 525
- Rodler, F., Lopez-Morales, M., & Ribas, I. 2012. ApJL, 753, L25
- Rothman, L.S., Gordon, I.E., Barbe, A., et al. 2009, J. Quant. Spectrosc. Radiat. Transfer, 110, 533
- Schwarz, H., Brogi, M., deKok, R.J., et al. 2015, A&A, 579, A111
- Snellen, I. A.G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
- Snellen, I.A.G., Brandl, B.R., deKok, R.J., et al. 2014, Nature, 503, 63
- Swain, M. R. Tinetti, G., Vashisht, G., et al. 2009, ApJ, 704, 1616.
- Tennyson, J., & Yurchenko, S. N. 2012, MRNAS, 425, 21
- Torres, G., Bakos, G. A., Hartman, J., et al. 2010, ApJ, 715, 458
- Valenti, J. A., Butler, R. P., & Marcy, G. W. 1995, PASP, 107, 966
- Vogt S. S., Allen S. L., Bigelow B. C. et al. 1994, Proc. SPIE, 2198, 362
- Wenger, M., Ochsenbein, F., Egret, D. et al. 2000, A&AS, 143, 1
- Wright J., Marcy G. W., Butler R. P. and Vogt S. S. 2004, ApJS, 152, 261
- Wright, J. T., G. W. Marcy, A. W. Howard, et al. 2012, ApJ, 753, 2
- Zucker, S. & Mazeh, T. 1994, ApJ, 420, 806