Thermal Spectrum of HD 88133 b

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


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

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
Table 1HD 88133 RV Measurements17
Figure 1.— RV data from the California Planet Survey with the best-fit stellar RV (primary velocity) curve overplotted in black. The colored points represent the NIRSPEC observations of this planet based on the observation phases and our qualitative expectations of their secondary velocities. The measure radial velocity of HD 88133 is 32.9 m/s

. In the course of this paper, we will show that the most likely value for the Keplerian orbital velocity of HD 88133 b is 40 15 km/s.

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


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
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)
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)

Table 2HD 88133 System Properties

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.

Figure 2.— Top-down schematic of the orbit of HD 88133 b around its star according to the orbital parameters derived by Fischer et al. (2005), Butler et al. (2006), and this work. Each point represents a single epoch’s worth of NIRSPEC observations of the system. Circles indicate L band observations and squares represent K band observations. The black arrow represents the line of sight to Earth.
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
Table 3NIRSPEC Observations of HD 88133 b

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.

Figure 3.— The data reduction and telluric correction process. (A): One order of a reduced AB pair of HD 881333 data taken on 2013 March 29 in L band, with a best-fit telluric spectrum overplotted with a green, dashed line. (B): The first principal component in arbitrary units of this time series of data which encapsulates changes in the stellar spectrum as the air mass varies during the observation. (C): The second principal component in arbitrary units which describes changes in abundances of telluric species. (D) The third principal component in arbitrary units which encompasses changes in plate scale. (E) Telluric-corrected data with the first five principal components removed shown in black. This is the data used for the cross-correlation analysis described in Section 4. Overplotted in orange is the stellar spectrum of HD 88133 adapted from the PHOENIX stellar library (Husser et al., 2013).

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 analysis18, we first reduce our data set in AB pairs and construct a data matrix having rows and columns, where is the number of AB pairs and is the number of pixels (1024 for individual NIRSPEC echelle orders). We linearly interpolate sub-pixel shifts between nods when aligning the AB pairs on the matrix grid, and then calculate the residual matrix according to


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 vapor19, carbon dioxide20, and methane21 abundances. For L band data, baseline values for ozone from a reference tropical model are also included for the orders ranging from 3.12-3.17 and 3.26-3.31 m. Next we calculate the covariance matrix of our mean-normalized data such that


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


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


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


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.

Figure 4.— Forward models for the planetary atmosphere of HD 88133 b produced by the PHOENIX and SCARLET models drawn at instrument resolution. Note that the flux calculated by the SCARLET model is shifted downward by 0.3 for clarity. Features shown here are principally due to water vapor. The correlation coefficient between these two models at zero-lag is 0.92.

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.

Figure 5.— Maximum likelihood functions for each observational epoch. (A): Maximum likelihood function of the stellar velocity shift of data taken on 2013 March 29. The black vertical dashed line indicates the expected stellar velocity shift. (B) - (G): Maximum likelihood function of for L band data from 3.0-3.4 m taken on 2012 April 1 and 3, 2013 March 10 and 29, 2014 May 14, and 2015 April 8, respectively. (H) - (J): Maximum likelihood function of for K band data from 2.10-2.20 m taken on 2012 November 21, 2015 December 1, and 2016 April 15, respectively. Note in (B) - (J), the blue dashed curve shows the maximum likelihood function for the PHOENIX model, the red curve shows the maximum likelihood function for the SCARLET model, and the grey vertical dashed lines indicate the planetary velocity shift on that date given an orbital solution having = 40 km/s. Based on , the error on the calculated planetary velocity shift is about 1.2 km/s.

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


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.

Figure 6.— Normalized log likelihood as a function of Keplerian orbital velocity . Note that the vertical axes cannot be directly compared, and the color scheme is the same as Figure 5. (A): Normalized log likelihood curve for six nights of L band data from 3.0-3.4 m. (B): Normalized log likelihood curve for three nights of K band data from 2.10-2.20 m. (C): Normalized log likelihood for all epochs and orders of NIRSPEC data used in panels (A) and (B). The grey region represents the one-sigma error bars determined by jackknife sampling for data cross-correlated with the SCARLET model. (D): Normalized log likelihood curve for six nights of L band data cross-correlated with shuffled planetary spectra.

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.

The authors would like to thank Heather Knutson for helpful discussions throughout the preparation of this manuscript. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This work was partially supported by funding from the NSF Astronomy & Astrophysics and NASA Exoplanets Research Programs (grants AST-1109857 and NNX16AI14G, G.A. Blake P.I.), and the Center for Exoplanets and Habitable Worlds, which is supported by the Pennsylvania State Unviersity, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. Basic research in infrared astrophysics at the Naval Research Laboratory is supported by 6.1 base funding. Finally, we thank an anonymous reviewer for helpful insights which improved the content of this paper.


  1. affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
  2. affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
  3. affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
  4. affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
  5. affiliation: King Abdullah University of Science and Technology, Thuwal Saudi Arabia
  6. affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125
  7. affiliation: Lunar and Planetary Laboratory, University of Arizona, Tuscon, AZ 85721
  8. affiliation: Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802
  9. affiliation: Center for Exoplanets and Habitable Worlds, Pennsylvania State University, University Park, PA 16802
  10. affiliation: Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125
  11. affiliation: Naval Research Laboratory, Washington, DC 20375
  12. affiliation: Department of Astronomy, Yale University, 260 Whitney Avenue, New Haven, CT 0651
  13. affiliation: Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125
  14. affiliation: Department of Astronomy, University of California Berkeley, Berkeley, CA 94720
  15. affiliation: Harvard-Smithsonian Center for Astrophysics; Institute for Theory and Computation, Cambridge, MA 02138
  16. affiliationtext: Institute for Astronomy, University of Hawaii at Manoa, Honolulu, HI 96822
  17. The full set of RVs for this system is available as an electronic table online.
  18. PCA tutorial,
  19. Caltech Submillimeter Observatory,
  20. Earth System Research Laboratory at Mauna Loa,
  21. Earth System Research Laboratory at Mauna Loa,


  1. Allard, N. F.,Allard, F., Hauschildt, P. H., Kielkopf, J. F., & Machin, L. 2003, A&A, 411, L473
  2. Asplund, M., Grevesse, N., Suval, J.A., Scott, P. 2009, ARAA, 47, 1
  3. Barman, T. S., Hauschildt, P. H., & Allard, F. 2001, ApJ, 556, 885
  4. Barman, T. S., Hauschildt, P. H., & Allard, F. 2005, ApJ, 632, 1132
  5. Benneke, B. 2015, arXiv:1504.07655
  6. Birkby, J. L., de Kok, R. J., & Brogi, M. et al. 2013, MNRAS, 436, L35
  7. Boogert, A. C. A., Blake, G. A., & Tielens, A. G. G. M. 2002, ApJ, 577, 271
  8. Borysow, A. 2002, A&A, 390, 4
  9. Brogi, M., Snellen, I.A.G., de Kok, R.J., et al. 2012, Nature, 486, 502
  10. Brogi, M., Snellen, I.A.G., de Kok, R.J., et al. 2013, ApJ, 767, 1
  11. Brogi M., de Kok R. J., Birkby J. L., Schwarz H. & Snellen I. A. G. 2014, A&A, 565A, 124B
  12. Brogi, M., de Kok, R.J., Albrecht, S. et al. 2016, ApJ, 817, 2
  13. Bryan, M.L., Knutson, H.A., Howard, A.W., et al. 2016, ApJ, 821, 2
  14. Burrows, A., & Volobuyev, M. 2003, ApJ, 583, 985
  15. Butler R. P., Marcy G. W., Williams E. et al 1996 PASP 108 500
  16. Butler, R.P., Wright, J.T., Marcy, G.W., et al. 2006, ApJ, 646, 505
  17. Chubak C., Marcy G., Fischer D. A. et al. 2012, arXiv:1207.6212
  18. de Kok R. J., Brogi M., Snellen I. A. G. et al. 2013, A&A, 554, A82
  19. Dulick, M., Bauschlicher, C. W., Jr., Burrows, A., et al. 2003, ApJ, 594, 1
  20. Fischer, D.W., Laughlin, G., Butler, P., et al. 2005, ApJ, 620, 481.
  21. Howard A. W., Johnson J. A., Marcy G. W. et al. 2010, ApJ, 721, 1467
  22. Howard A. W., Johnson J. A., Marcy G. W. et al. 2011, ApJ, 730, 10
  23. Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  24. Hut, P. 1981, A&A, 99, 1
  25. Isaacson H. & Fischer D. 2010, ApJ, 725, 875
  26. Johnson J. A., Howard A. W., Marcy G. W. et al. 2010, PASP, 122, 888
  27. Kass, R.E., & Raftery, A.E. 1995, JASA, 90,773
  28. Lockwood, A.C., Johnson, J.A., Bender, C.F., et al. 2014, ApJ, 783, L29
  29. 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.
  30. Martins, J.H.C., Santos, N.C., Figueira, P., et al. 2015, A&A, 576, A134
  31. Mayor, M. & Queloz, D. 1995, Nature, 378, 23
  32. McLean, I. S., Becklin, E. E., Bendiksen, O., et al. 1998, Proc. SPIE 3354, 566
  33. Mortier, A., Santos, N.C., Souse, S.G., et al. 2013, A&A, 557, A70
  34. Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&A S, 112, 525
  35. Rodler, F., Lopez-Morales, M., & Ribas, I. 2012. ApJL, 753, L25
  36. Rothman, L.S., Gordon, I.E., Barbe, A., et al. 2009, J. Quant. Spectrosc. Radiat. Transfer, 110, 533
  37. Schwarz, H., Brogi, M., deKok, R.J., et al. 2015, A&A, 579, A111
  38. Snellen, I. A.G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  39. Snellen, I.A.G., Brandl, B.R., deKok, R.J., et al. 2014, Nature, 503, 63
  40. Swain, M. R. Tinetti, G., Vashisht, G., et al. 2009, ApJ, 704, 1616.
  41. Tennyson, J., & Yurchenko, S. N. 2012, MRNAS, 425, 21
  42. Torres, G., Bakos, G. A., Hartman, J., et al. 2010, ApJ, 715, 458
  43. Valenti, J. A., Butler, R. P., & Marcy, G. W. 1995, PASP, 107, 966
  44. Vogt S. S., Allen S. L., Bigelow B. C. et al. 1994, Proc. SPIE, 2198, 362
  45. Wenger, M., Ochsenbein, F., Egret, D. et al. 2000, A&AS, 143, 1
  46. Wright J., Marcy G. W., Butler R. P. and Vogt S. S. 2004, ApJS, 152, 261
  47. Wright, J. T., G. W. Marcy, A. W. Howard, et al. 2012, ApJ, 753, 2
  48. Zucker, S. & Mazeh, T. 1994, ApJ, 420, 806
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
Request answer
The feedback must be of minumum 40 characters
Add comment
Loading ...