NEAR-IR DIRECT DETECTION OF WATER VAPOR IN TAU BOÖTIS b

Since the first detection of an exoplanet around a main-sequence star (Mayor & Queloz 1995), astronomers have discovered hundreds of exoplanets using the radial velocity (RV) technique (Wright et al. 2012). This powerful tool for discovering exoplanets only measures the minimum, or "indicative" mass, M_psin i, leaving a degeneracy between two interesting properties of the system, one physical and one orbital. While the indicative mass is useful for statistical studies, uniquely measuring the true mass of an exoplanet not only yields a key physical property, but also furthers our understanding of planetary formation and evolution via measurement of the true mass distribution of exoplanets (Zucker & Mazeh 2001; Weiss et al. 2013). For example, Ho & Turner (2011) demonstrate that knowledge of the true mass distribution is necessary to convert a minimum mass, M_psin i, into an estimate of a planet's true mass for RV-detected systems.

One class of objects for which the planet mass can be directly determined are those that transit their host star. Hundreds of transiting planets have been discovered and characterized, and the ongoing Kepler mission has found potential exoplanet candidates numbering in the thousands (Borucki et al. 2011; Batalha et al. 2013). Transit events also provide atmospheric information through transmission spectroscopy and secondary eclipses. Investigators have measured spectra of some of the larger transiting planets (Charbonneau et al. 2002), leading to the discovery of species such as water, methane, and carbon monoxide (Knutson et al. 2012; Berta et al. 2012; Crossfield et al. 2013; Baskin et al. 2013). A variety of spectral retrieval methods have been used to verify the results, confirming the hypothesis that O- and C-bearing gases are present in the atmospheres of hot Jupiters (Madhusudhan & Seager 2009; Line et al. 2012).

Recently, a technique previously used to detect low mass ratio spectroscopic binary stars has been applied to stars known to host exoplanet systems. The direct RV detection of an exoplanet involves separating the planetary and stellar components spectroscopically. The high flux ratio between the primary and companion makes these detections difficult, but not impossible, thanks to modern infrared echelle spectrographs. Snellen et al. (2010) detected CO on HD 209458 b with a precision of 2 km s⁻¹ using the R = 100,000 Very Large Telescope CRIRES instrument. This detection provided the spectroscopic orbit of the system and determined the true mass of the planet.

Indeed, this technique can be applied to non-transiting, RV-detected exoplanets to extract the unknown inclination and true mass. CRIRES was also used to detect CO on τ Boötis b (Brogi et al. 2012), the first ground-based detection of a short-period non-transiting exoplanet atmosphere, a result confirmed shortly thereafter by Rodler et al. (2012). These studies provide the true mass of the planet and probe the chemical composition of its atmosphere. A combination of high signal-to-noise, high spectral resolution, and coverage of multiple CO overtone lines was required to achieve the sensitivity required for these detections.

Despite the agreement between the two groups, the direct detection of exoplanets, especially τ Boo b, has a storied history. Collier Cameron et al. (1999) first reported a detection of reflected light from τ Boo b more than a decade ago and Wiedemann et al. (2001) reported a possible detection of CH₄ in the planet's atmosphere soon thereafter, using the same orbital solution. However, the planetary velocity from the earlier results disagrees with more recent findings discussed above.

Here we report a detection of water vapor in the atmosphere of τ Boo b, using spectroscopic observations centered around 3.3 μm. The τ Boötis system is a stellar binary comprised of an F-type "A" component and an M-dwarf "B" component. The hot Jupiter τ Boo b orbits the larger "A" component with a period of 3.312 days and a minimum (indicative) mass Msin i = 3.87 M_Jup (Butler et al. 1997). Our detection confirms the mass determined by Brogi et al. (2012) and serves to further characterize the atmospheric chemistry of this exoplanet. Data from five epochs reveal orbital motion of τ Boo b that is consistent with that reported by Brogi et al. (2012). The planetary template spectrum used in the cross-correlation is dominated by water vapor opacity, providing strong evidence of water in the atmosphere of a non-transiting hot Jupiter for the first time.

2.1. Observations and Data Reductions

We observed the τ Boo system on five separate nights in 2011 March, 2011 May, and 2012 April (Table 1), chosen to optimize phasing near orbital quadrature and allowing for maximum separation of stellar and planetary lines. We use the Near Infrared Echelle Spectrograph (NIRSPEC; McLean et al. 1995) at the W. M. Keck Observatory, which provides high resolution (R ∼ 25,000 for a three-pixel slit) in multiple orders at the wavelengths of interest to study near-infrared water emission. We obtained spectra covering 3.404–3.457 μm, 3.256–3.307 μm, 3.121–3.170 μm, and 2.997–3.043 μm. Each epoch covers a total elapsed time of approximately 1 hr and is comprised of a continuous sequence of hundreds of exposures.

Table 1. L-band Observations of Tau Boo b

Date	JD 2450000	Phase	Vbary	S/N3 μm
		(rad)	(km s−1)
2011 May 21	5702.8542	0.3199	−17.36	6525
2012 Apr 3	6021.0625	0.3847	2.26	4768
2012 Apr 1	6019.0625	0.7810	3.02	2871
2011 Mar 14	5634.9375	0.8164	11.75	5467
2011 Mar 24	5644.9375	0.8353	7.39	6575

Download table as:  ASCIITypeset image

We extract NIRSPEC spectra using a custom Interactive Data Language optimal extraction pipeline, similar to that described by Cushing et al. (2004). The spatial profile weighting intrinsic to optimal extraction provides a reliable method for detecting and removing bad pixels due to detector defects and cosmic rays. It adjusts for seeing variations that occur over the course of an observation, and can also minimize the contamination from nearby stars that happen to fall on the slit. We have determined that the contamination from the M-dwarf companion τ Boo B is negligible in all epochs of our extracted τ Boo A spectra. Arc lamps typically used for wavelength calibration do not provide useful reference lines in the thermal infrared. Instead, we wavelength calibrate our spectra using unblended telluric features with accurately known rest wavelengths. Between 15 and 30 individual telluric lines are used for each order.

We first correct the bulk telluric absorption with TERRASPEC (Bender et al. 2012), a synthetic forward-modeling algorithm that uses the line-by-line radiative transfer model LBLRTM (Clough et al. 2005) to generate a synthetic absorption function for the Earth's atmosphere (hereafter TAF). The TAF includes continuum absorption and line absorption for 28 molecular species, with line information provided by the HITRAN 2008 database (Rothman et al. 2009). Vertical mixing profiles for the seven most prominent species (H₂O, CO₂, O₃, N₂O, CO, CH₄, and O₂) can be scaled or adjusted from the LBLRTM default profiles, which include the US Standard 1976, tropical, midlatitude, and subartic models. These models represent the average atmosphere for their respective latitudes (in 1976), which is a good initial guess for the instantaneous atmosphere corresponding to a single observation. We use the tropical model to provide initial vertical profiles for spectroscopy obtained from Mauna Kea.

TERRASPEC convolves the TAF with an instrumental broadening profile (hereafter IP), which is measured from the data. The IP is parameterized as a central Gaussian surrounded by satellite Gaussians offset by a fixed percentage of the Gaussian width, and with adjustable amplitudes. This is similar to the IP parameterization described by Valenti et al. (1995) for use with I₂ absorption cells. Typically two to four satellite Gaussians serve to model the IP. The instrumental broadened TAF is then multiplied by a low-order wavelength dependent polynomial to correct for the combined effects of blaze sensitivity and stellar continuum. TERRASPEC uses the least-squares fitting algorithm MPFIT (Markwardt 2009) to optimize the parameters comprising the TAF, IP, and continuum. The FWHM of the profile is consistently 1.3 × 10⁻⁴ μm, which yields R = λ/Δλ = 24, 000. Spectral regions containing stellar absorption are excluded from the optimization by adaptive masks.

The initial telluric correction can be significantly affected by the instrumental fringing. We therefore mask strong stellar lines from the telluric corrected spectrum, and analyze the remaining spectral range with a Lomb–Scargle periodogram (Scargle 1982) to measure the frequency and power of the individual fringes present in our spectra. Two prominent fringes are seen at ∼1.75 cm⁻¹and 2.18 cm⁻¹. We then calculate the composite fringe function as the product of the individual fringes, and divide it into the non-telluric-corrected spectrum. The defringed spectrum is then reprocessed with TERRASPEC, yielding the telluric-free and continuum-normalized stellar spectrum, a measurement of the TAF, and a parameterization of the IP. This process was often iterated two to three times to ensure full removal of fringing and telluric absorption. Regions with atmospheric transmission ⩽40% of the continuum value, determined from the TAF, are masked and excluded from the final planet search. Figure 1 demonstrates the telluric removal process.

Zoom In Zoom Out Reset image size

2.2. Two-dimensional Cross Correlation

2.3. Maximum Likelihood Analysis

To find the most likely solution, each CCF must be converted into a probability, or likelihood, $\mathcal {L}$. To do this, we start with a relationship between $\mathcal {L}$ and the familiar χ² statistic

$$\begin{eqnarray} \mathcal {L} &=& \frac{1}{\sqrt{2 \pi } \sigma _{i}}\prod _{i} \exp { \left(- \frac{ \chi ^{2}_{i}}{2}\right)} \nonumber \\ \log \mathcal {L} &=& \kappa - \frac{\chi ^{2}}{2}, \end{eqnarray} \tag{ 1 }$$

where κ is a constant that does not matter when comparing relative likelihoods, assuming σ_i = σ = const. Denoting the observed spectra as S_i and the template spectra as f_i ≡ f(λ_i + Δλ) yields

$$\begin{eqnarray} \chi ^{2} &=& \sum _i \frac{(S_i - f_i)^2}{\sigma ^{2}} \nonumber \\ &=& \frac{1}{\sigma ^2} \sum _i \left(S_i^2 + f_i^2 - 2 S_i f_i\right). \end{eqnarray} \tag{ 2 }$$

Each part of the quotient can be summed individually. For a continuum-normalized spectrum with N points, this results in

$$\begin{equation} \sigma ^{2} = \frac{1}{N}\sum S_{i}^2 \end{equation} \tag{ 3 }$$

$$\begin{equation} \chi ^{2}_{i} = \mathrm{const} - \sum \frac{S_{i} f_i}{\sigma ^{2}}. \end{equation} \tag{ 4 }$$

The second term on the right-hand-side of Equation (4) is simply the CCF. Thus, the CCF and $\mathcal {L}$ are related by

$$\begin{equation} \log \mathcal {L} = \mathrm{const} + {\rm CCF.} \end{equation} \tag{ 5 }$$

The goal is to maximize ∑log L for both the stellar and planetary velocity shifts, where the sum is over all spectral orders for an individual epoch. Figure 2 (top panel) demonstrates that the stellar velocity clearly stands out as the most likely solution for one epoch, as is the case for all other epochs as well. Our NIRPSEC spectra have insufficient RV precision to be sensitive to the orbital motion of τ Boo A so the stellar RV is therefore fixed at the systemic value. Barycentric movement is accounted for.

Zoom In Zoom Out Reset image size

The ML solution of the planet's orbit is much more complex. The likelihood is proportional to the CCF, a two-dimensional surface that reflects coherence in features between template spectra and those of the target. Due to the multiplicity of rovibrational transitions in the asymmetric top spectrum of hot water, the correlation coefficient remains large (>0.9) over significant offsets. This results in multiple "peaks" at incorrect velocity lags, as can be seen in Figure 2. Thus, a single epoch would lead to a degeneracy of solutions. Fortunately, we do not worry about alignment between planet and telluric spectra because not only are the spectra very different, but the combined systemic and planetary RV shifts the planetary spectrum sufficiently to avoid collisions between spectral features.

2.4. Orbital Solution

To break the degeneracy several epochs of data must be brought to bear. Using the known period of the planet, we calculate the orbital phase of each epoch (Table 1) and use that to find the most likely planetary velocity consistent with a circular Keplerian orbit (the estimated planetary eccentricity is small, e ∼ 0.02; Butler et al. 1997; Brogi et al. 2012). Since the absolute orbital velocity of the planet is known from the period, we are actually interested in K_p, or semi-amplitude. A range of K_p is tested, each corresponding to a different inclination of the planet, as well as a unique mass. Furthermore, each K_p leads to a different velocity lag at a given phase, such that

$$\begin{equation} v_{{\rm pl}} = K_{p} \sin (\omega t + \phi) + \gamma, \end{equation} \tag{ 6 }$$

where ω = 2π/P, ϕ is a phase lag that can be set to zero by choosing the proper starting date, and γ is the combined stellar barycentric and systemic velocities (Vsin i ≈ 15 km s⁻¹ (Butler et al. 2006) and V_bary listed in Table 1). We then seek the best-fitting K_p by maximizing the sum of the likelihood of v_pl for all epochs, whose highest value indicates the most likely solution for K_p.

It is clear that the ML curve for each epoch in Figure 2 shows multiple peaks. Only one of these peaks per epoch can be the true peak that results from the planet signal in our spectra; the other peaks are artifacts caused by chance misalignments. The molecular lines in the model water spectrum can align by chance with both signal and noise features in the observed spectra and with other molecules in the planetary atmosphere that have not been included in the planetary template. However, the true peaks can be distinguished from artifacts by requiring a multi-epoch orbital solution that is consistent with a Keplerian orbit. The dashed, dotted, and solid lines in Figure 2 represent the v_pl for three Keplerian solutions with K_p = 70 km s⁻¹, 90 km s⁻¹, and 110 km s⁻¹, respectively. While none correspond to the highest peak for all epochs, certain orbital solutions find peaks more often than troughs.

To find the most likely orbit, the likelihoods at every epoch for a given K_p are combined. The second panel in Figure 3 presents this composite ML, along with the corresponding planetary mass. This function represents the sum of the log likelihoods of all five epochs. There are two peaks in the likelihood function, at K_p ∼ 70 km s⁻¹ and 111 km s⁻¹. We have demonstrated empirically the degenerate effects of cross-correlating a water spectrum. Now we explore the systematics of the analysis and consider how both the properties of the individual stellar and planetary template autocorrelation functions, as well as the correlation of the two template spectra, can affect the composite likelihood function.

Zoom In Zoom Out Reset image size

To do this, synthetic data sets are subjected to the same analysis as the original data. Synthetic target spectra of just the stellar and planetary templates are used, each shifted to the correct velocity for a given epoch. The results are given in Figure 3, whose top panel shows the normalized log likelihood from a planetary spectrum injected at 70 km s⁻¹, with a planet-to-star flux ratio of 10⁻⁵. The proper signal is clearly retrieved. However, the third panel shows that at an injected planetary RV of 110 km s⁻¹ and the same contrast ratio, a double-peaked log likelihood results. Thus, following exactly the same procedure as was used for the observations, a "perfect" target spectrum with no noise and no (terrestrial) atmosphere retrieves both the correct and an additional signal, with a structure that mimics the data, likely as a result of the complex hot water spectrum near 3 μm.

At a sufficiently high planet-to-star flux ratios, the true planetary signal should dominate the posterior likelihood. Indeed, using a flux ratio of 10⁻³, Birkby et al. (2013) detect water absorption from another (transiting) hot Jupiter, HD 189733 b. The right panel of Figure 4 shows that when the flux ratio of the planetary signal is increased to this level, the singular correct velocity is retrieved. We can constrain the spectroscopic contrast of τ Boo b relative to its host star by comparing the data to these synthetic fits. The simulations show that for α ⩾ 10⁻⁴, the correct 110 km s⁻¹ solution should present a larger ML. Since our data do not demonstrate this, we can put a 1σ upper limit on the high resolution spectroscopic flux ratio of 10⁻⁴ at ∼3.3μm. The left panel of Figure 4 shows that a planetary signal of 70 km s⁻¹ would be uniquely retrieved for all realistic flux ratios.

Zoom In Zoom Out Reset image size

The detection of τ Boo b has been a difficult quest. Contradictory results have been published over the past 15 years and the difficulty in both thoroughly removing the telluric absorption and also extracting the diminutive planetary signal is not to be underestimated. Multiple observations at different wavebands improve the validity of the detection, and facilitate a variety of statistical tests to ensure an accurate measurement.

Here, the validity of the most likely solutions for the RV of τ Boo b are explored. Of the two prominent velocities that fit the data, one is shown to result from the systematics of the cross-correlation analysis performed on this high flux ratio spectroscopic binary system. Interestingly, previous velocity studies of the planet in reflected light (Collier Cameron et al. 1999) and methane (Wiedemann et al. 2001) retrieved values close to this artifact. We suggest that residual water vapor in the atmosphere, after the telluric removal, might have been responsible for the false value previously reported from infrared observations. As methods for telluric corrections have improved over the past decade, this problem can be overcome and the correct RV retrieved.

Our analysis gives a 6σ detection of the planet at K_p = 111 ± 5 km s⁻¹ for τ Boo b, with a 1σ upper-limit on the 3.3 μm planet/star spectroscopic flux ratio of 10⁻⁴. To determine the significance of our detection, we injected synthetic signals at a variety of planetary velocities and planet-to-star spectroscopic flux ratios and constructed the chi-square surface of the ML fits. The orbital velocity is in good agreement with previous RV amplitude detections via CO by Brogi et al. (2012) and Rodler et al. (2012), and our analysis reveals the presence of water vapor in the planet's atmosphere. Furthermore, using a stellar mass of 1.341$^{+0.054}_{-0.039}$ M_☉ (Takeda et al. 2007) and a stellar velocity semi-amplitude of 0.4664 ± 0.0033 km s⁻¹ (Brogi et al. 2012), we derive a planetary mass of $M_{P} = 5.90^{+0.35}_{-0.20}\,M_{\rm Jup}$ with an orbital inclination of $i = 45^{+3}_{-4}$°.

The technique presented here is in its nascent stages, and the work is by no means complete. Additional quantitative characterization of the physical properties (temperature, opacity) and composition (and thus estimates of vertical mixing) of τ Boo b's atmosphere will require both significant simulations of the atmospheric radiative transfer and data analysis along with additional data at longer and shorter wavelengths. Such work is beyond the scope of this Letter, but expanded studies of τ Boo using the same methods described above, but also including molecules such as CH₄ and CO at longer wavelengths than examined here, are underway. Indeed, although the mole fractions of methane should be insignificant for τ Boo b, we have searched for the molecule using the methods described herein and have no detection to report in our data. With careful analysis in the future, the relative abundances of these molecules in the atmospheres of non-transiting exoplanets can be ascertained. Further applications of this technique, using water and methane as template molecules, to a variety of additional hot Jupiter exoplanets will be reported elsewhere.

The W. M. Keck Observatory is operated as a scientific partnership among the California Institute of Technology, the University of California and NASA, and was made possible by the financial support of the W. M. Keck Foundation. A.C.L. and G.A.B. gratefully acknowledge support from the NSF GRFP and AAG programs, JAJ the generous grants from the David and Lucile Packard and Alfred P. Sloan Foundations, and C.F.B. support from the Center for Exoplanets and Habitable Worlds, which is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. Basic research in infrared astronomy at the Naval Research Laboratory is supported by 6.1 base funding. We thank Jacques Sauval for kindly providing a copy of his solar linelist. Finally, the authors wish to acknowledge the significant cultural role of the summit of Mauna Kea.