EVIDENCE FOR THE DIRECT DETECTION OF THE THERMAL SPECTRUM OF THE NON-TRANSITING HOT GAS GIANT HD 88133 b

Data were taken on six nights (2012 April 1 and 3, 2013 March 10 and 29, 2014 May 14, 2015 April 8) in the L band and three nights (2015 November 21, 2015 December 1, and 2016 April 15) in the K band in an ABBA nodding pattern with the NIRSPEC instrument at the W.M. Keck Observatory (McLean et al. 1998). With a 04 × 24'' slit NIRSPEC has an L-band resolution of 25,000 (30,000 in the 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.

Zoom In Zoom Out Reset image size

Table 3.  NIRSPEC Observations of HD 88133 b

Date	Julian Date	Mean Anomaly M	True Anomaly f	Barycentric Velocity vbary	Integration Time	S/NL,Ka
	(−2,440,000 days)	(2π rad)	(2π rad)	(km s−1)	(minutes)
L band (3.0–3.4 μm)
2012 Apr 1	16018.837	0.89	0.86	−20.96	140	1680
2012 Apr 3	16020.840	0.48	0.48	−21.66	140	2219
2013 Mar 10	16361.786	0.29	0.33	−11.59	180	2472
2013 Mar 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 Apr 8	17120.835	0.50	0.50	−23.01	160	2938
K band (2.0–2.4 μm)
2015 Nov 21	17348.129	0.06	0.68	29.95	60	2701
2015 Dec 1	17358.117	0.96	0.62	29.25	60	2823
2016 Apr 15	17166.300	0.54	0.53	−29.15	80	2466

Note.

^aS/N${}_{L}$ and S/N${}_{K}$ are calculated at 3.0 μm and 2.1515 μm, respectively. Each S/N calculation is for a single channel (i.e., resolution element) for the whole observation.

Download table as:  ASCIITypeset image

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 $\lambda ={{ax}}^{3}+{{bx}}^{3}+{cx}+d$, x is pixel number, and a, b, c, and d 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.

Zoom In Zoom Out Reset image size

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: 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 approximately constant. Note that we observe for only 2–3 hr 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 PCA,¹⁰ we first reduce our data set in AB pairs and construct a data matrix ${\boldsymbol{X}}$ having n rows and m columns, where n is the number of AB pairs and m 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 ${\boldsymbol{R}}$ according to

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{{ij}}=\displaystyle \frac{{{\boldsymbol{X}}}_{i}-M}{{\sigma }_{j}}\end{eqnarray} \tag{ 2 }$$

where i is the row number, j is the column number, M is either the mean of X_j or a telluric model, and ${\sigma }_{i}$ is the standard deviation of the values in column j. We guide our PCA with a telluric model for M (rather than the mean of X_j) that uses baseline values for water vapor,¹¹ carbon dioxide,¹² and methane¹³ 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 ${\boldsymbol{C}}$ of our mean-normalized data such that

$$\begin{eqnarray}&&{\boldsymbol{C}}=\displaystyle \frac{{{\boldsymbol{R}}}^{T}{\boldsymbol{R}}}{n-1}.\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&{\boldsymbol{C}}={{\boldsymbol{USV}}}^{T}\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{U}}$ contains the left singular vectors (or the eigenvectors), ${\boldsymbol{S}}$ is a diagonal matrix of the singular values (or the eigenvalues), and ${\boldsymbol{V}}$ 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 k principal components of the data set, given by ${{\boldsymbol{U}}}_{k}{{\boldsymbol{S}}}_{k}$, where ${{\boldsymbol{U}}}_{k}$ is the first k columns of ${\boldsymbol{U}}$ and ${{\boldsymbol{S}}}_{k}$ is the k × k upper-left part of ${\boldsymbol{S}}$. Rank-k data, ${{\boldsymbol{X}}}_{k}$, can be built as

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{k}={{\boldsymbol{U}}}_{k}{{\boldsymbol{S}}}_{k}{{\boldsymbol{U}}}_{k}^{T}.\end{eqnarray} \tag{ 5 }$$

To produce a telluric-corrected spectrum, ${{\boldsymbol{X}}}_{\mathrm{corr},i}$, each row ${{\boldsymbol{X}}}_{i}$ of ${\boldsymbol{X}}$ is divided by its corresponding un-mean normalized row in ${{\boldsymbol{X}}}_{k,i}$:

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{\mathrm{corr},i}=\displaystyle \frac{{{\boldsymbol{X}}}_{i}}{{{\boldsymbol{X}}}_{k,i}{\sigma }_{j}+M}.\end{eqnarray} \tag{ 6 }$$

Finally, we collapse the rows of data in ${{\boldsymbol{X}}}_{\mathrm{corr}}$ 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 principal 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 the 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 to 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 15 components to delete the planet signal. We have experimented with incrementally removing up to 10 principal components as input to the subsequent analysis; the results and conclusions presented below use data with the first five principal components removed.

Our synthetic stellar spectrum is a PHOENIX model (Husser et al. 2013) interpolated to match the published effective temperature ${T}_{\mathrm{eff}}$, surface gravity $\mathrm{log}g$, and metallicity $[{\rm{Fe}}/{\rm{H}}]$ of HD 88133 listed in Table 2. As HD 88133 has a $v\sin i$ < 5 km s⁻¹, instrumental broadening dominates, and we convolve the stellar model with the instrumental profile calculated in Section 3.2.

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 > 250 K, 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.

Zoom In Zoom Out Reset image size

The SCARLET model considers the molecular opacities of ${{\rm{H}}}_{2}{\rm{O}}$, CH₄, NH₃, HCN, $\mathrm{CO}$, CO₂, and $\mathrm{TiO}$ from the high-temperature ExoMol database (Tennyson & Yurchenko 2012), and O₂, O₃, $\mathrm{OH}$, ${{\rm{C}}}_{2}{{\rm{H}}}_{2}$, ${{\rm{C}}}_{2}{{\rm{H}}}_{4}$, ${{\rm{C}}}_{2}{{\rm{H}}}_{6}$, ${{\rm{H}}}_{2}{{\rm{O}}}_{2}$, and HO₂ 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 H₂-broadening prescription provided in Burrows & Volobuyev (2003). Collision-induced broadening from ${{\rm{H}}}_{2}/{{\rm{H}}}_{2}$ and ${{\rm{H}}}_{2}/\mathrm{He}$ 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 line lists 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 line list 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 ${\alpha }_{\mathrm{phot}}$ (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 ${\alpha }_{\mathrm{spec}}$. Since the cross-correlation analyses described in Section 4.3 are not very sensitive to contrast ratios, varying the value of ${\alpha }_{\mathrm{spec}}$ does not change our conclusions on the RV of the planet, and thus the system inclination (Lockwood et al. 2014). However, the specific value of ${\alpha }_{s}$ 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 ${\alpha }_{\mathrm{spec}}$, although the size and structure of the final maximum likelihood peak near the planet's signature will be.

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 to 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 to 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 CCF 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⁻¹.

Zoom In Zoom Out Reset image size

The retrieval of the planet velocity v_sec 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.

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

$$\begin{eqnarray}&&{v}_{\sec }(f)={K}_{p}(\cos (f+\omega )+e\cos \omega )+\gamma \end{eqnarray} \tag{ 7 }$$

where K_p is the planet's Keplerian orbital velocity, ω is the longitude of periastron measured from the ascending node, e 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.

Zoom In Zoom Out Reset image size

When the six maximum likelihood curves produced from L-band data (panels B–G in Figure 5) are combined with equal weighting according to Equation (7), we see that the most likely value for the radial projection of the planet's Keplerian orbital velocity K_P 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 Figure 5) suggests that K_P is 32 ± 12 km s⁻¹.

The combination of all nine nights of data yields K_P = 40 ± 15 km s⁻¹. It is this value of K_P that we use to calculate the secondary velocity curve shown in the bottom panel of Figure 2. The peak near K_P = 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 among 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 therefore 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 B 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. $2\mathrm{ln}B$ must be greater than 10 for a model to be very strongly preferred.

For the Gaussian distribution compared to the straight line, $2\mathrm{ln}B$ 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 K_P with the parameters given in Table 2, a K_P of 40 ± 15 km s⁻¹ implies that the true mass of HD 88133 b is ${1.02}_{-0.28}^{+0.61}{M}_{{\rm{J}}}$ and its orbital inclination is ${15}_{-5}^{+6^\circ }$.

Note that the values of v_sec implied by the most likely value of K_P often, but do not always, correspond with peaks in the maximum likelihood curves for each night, as indicated by the vertical gray 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.

We first check our detection of HD 88133 b's emission spectrum at a K_P RV projection of 36 km s⁻¹ by varying the spectroscopic contrast ${\alpha }_{\mathrm{spec}}$ uniformly with orbital phase. We tested nine values of ${\alpha }_{\mathrm{spec}}$ from 10⁻⁷ to 10⁻³, and find that the maximum likelihood peak near 40 km s⁻¹ is robust for ${\alpha }_{\mathrm{spec}}\geqslant {10}^{-5.5}$.

We create a "shuffled" planetary model by randomly rearranging chunks of each planetary atmosphere model. If the maximum likelihood peak at 36 km s⁻¹ in the L band data 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 K_P. 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 L-band 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 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 PCA 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}_{-0.24}^{+0.38}{M}_{{\rm{J}}}$, and an orbital inclination of ${50}_{-4}^{+3^\circ }$ for tau Boo b. This is in good agreement with Lockwood et al. (2014) (K_P = 111 ± 5 km s⁻¹) as well as Brogi et al. (2012) (K_P = 110 ± 3.2 km s⁻¹) and Rodler et al. (2012) (K_P = 115 ± 11 km s⁻¹), thereby validating our PCA-based methods.

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 K_P preferred by each data set are consistent. Although 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.

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 K_P 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.

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.