Detection of Water Vapor in the Thermal Spectrum of the Non-transiting Hot Jupiter Upsilon Andromedae b

We used NIRSPEC (Near InfraRed SPECtrometer; McLean et al. 1998) at the Keck Observatory to observe ups And A and b on seven nights (2011 September 6, 7, and 9, 2013 October 27 and 29 and November 7, and 2014 October 7) in L, three nights (2016 September 19, November 12, and December 15) in K_r (the right, long-wavelength half of the dispersed, K band filtered light), and three nights (2014 October 5 and 2016 August 21 and September 19) in K_l (the left, short-wavelength half of the dispersed, K band filtered light). We obtained spectral resolutions of ∼25,000 in L and ∼30,000 in K using the 04 × 24'' slit setup and used an ABBA nodding pattern during data acquisition. In the L band, the echelle orders typically cover 3.4038–3.4565/3.2467–3.3069/3.1193–3.1698/2.995–3.044 μm. The echelle orders in the K_r band typically cover 2.38157–2.41566/2.31–2.34284/2.24245–2.27485/2.17894–2.20861/2.11878–2.14639/2.06170–2.08703 μm, while in the K_l band the echelle orders typically cover 2.34238–2.37535/2.27198–2.30374/2.20554–2.23653/2.14362–2.17298/2.08461–2.11312/2.02931–2.05634 μm. Altogether, the two K band setups provide near continuous wavelength coverage across the entire K band. Table 2 gives the details of these 13 observations.

Table 2.  NIRSPEC Observations of ups And b

Date	Modified Julian Datea	Mean Anomaly Mb	Barycentric Velocity vbary	Integration Time	S/NL,Kc
	(−2400,000.5 days)	(2π rad)	(km s−1)	(min)
L Band (3.0–3.4 μm)
2011 Sep 6	55810.639	0.25	21.07	60	5376
2011 Sep 7	55811.637	0.46	20.82	10d	2661
2011 Sep 9	55813.509	0.87	20.33	100	8265
2013 Oct 27	56592.526	0.59	1.89	140	9173
2013 Oct 29	56594.512	0.02	0.99	140	5937
2013 Nov 7	56603.609	0.99	−3.17	180	8686
2014 Oct 7	56937.553	0.32	10.64	50	5641
Kr Band (Long Wavelength Side of 2.0–2.4 μm)
2016 Sep 19	57650.361	0.73	17.14	100	11517
2016 Nov 12	57704.265	0.38	−5.37	230	12872
2016 Dec 15	57737.300	0.53	−18.63	70	7666
Kl Band (Short Wavelength Side of 2.0–2.4 μm)
2014 Oct 5	56935.579	0.87	11.47	70	7764
2016 Aug 21	57621.589	0.45	24.15	30	4369
2016 Sep 19	57650.501	0.73	17.14	120	10649

Notes.

^aJulian date refers to the middle of the observing sequence. ^bWe list only the mean anomalies (and no true anomalies) for our observations, since ups And b's orbit is nearly circular. ^cS/N_L, S/N${}_{{K}_{r}}$ and, S/N${}_{{K}_{l}}$ are calculated at 3.0, 2.1325, and 2.1515 μm, respectively. Each S/N calculation is for a single channel (i.e., resolution element) for the whole observation. ^dBecause the total integration time on ups And on 2011 September 7 is very short, we do not use principal component analysis to remove the telluric signals (see Section 2.2), and we exclude this epoch from the following analysis.

Download table as:  ASCIITypeset image

A top–down schematic of ups And b in orbit around ups And A is shown in Figure 1 with the expected orbital phase of each observational epoch marked. Figure 2 shows RV measurements of ups And A taken from Fischer et al. (2014) in comparison with expectations for the line-of-sight velocity of ups And b. We aim to take observations when the line-of-sight velocities of the star and planet are most distinct and when we expect to observe a decent amount of dayside radiation from the planet, thus maximizing the planet flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our data reduction and cleaning methods are parallel to those described in Piskorz et al. (2016) and are summarized here only briefly.

We use a Python pipeline in the style of Boogert et al. (2002) to flat field and dark subtract our data, remove bad pixels, and extract 1D spectra. For the wavelength calibration, we fit a fourth-order polynomial ($\lambda ={{ax}}^{3}+{{bx}}^{2}+{cx}+d$, where x is pixel number and a, b, c, and d are free parameters) that aligns our L band data to a telluric model or our K band data to a combined telluric and stellar model. Here, the difference in treatment of L and K band data stems from the fact that telluric lines are stronger in the L band than near 2 μm. Our stellar model is derived from the PHOENIX stellar library (Husser et al. 2013) and is described in more detail in Section 3.1. Finally, we fit an instrument profile to our data as in Valenti et al. (1995) and save it to apply to the models described in Sections 3.1 and 3.2.

We capitalize on the long time series of observations (roughly two minutes per nod, or four minutes per AB pair) taken at each epoch and perform a principal component analysis (PCA) to remove contributions to the spectra from the Earth's atmosphere. PCA rewrites a data set in terms of its principal components so that the variance of a data set with respect to a model or its mean is reduced. For our purposes, this means that PCA will identify the time-varying components of our time-series data, most notably, changes in the telluric spectrum over the course of a given epoch. The first principal component describes the most variance, the second, the second most, etc. We guide our PCA with a telluric model that best fits the data in terms of water, carbon dioxide, methane, and (where appropriate) ozone abundances, and determine the eigenvectors making up each observed spectrum. We calculate and remove the strongest principal components from our data, leave behind the parts of the spectra that are constant in time (the stellar and planet signals), combine every AB nod of data, and clip regions of substantial telluric absorption (>75%). More information on this technique is given in Piskorz et al. (2016). Figure 3 shows a raw spectrum of ups And taken on 2013 October 29, the first three principal components, and a cleaned spectrum of ups And.

Zoom In Zoom Out Reset image size

As in our analysis of HD 88133 data, we find that the telluric correction by PCA works well for all orders of L band data, but poorly for the K_r and K_l band orders spanning 2.06170–2.08703 μm and 2.02931–2.05634 μm, where there is a dense forest of telluric CO₂ lines. We also find that a few nights of K band observations were contaminated by significant issues with the read-out electronics. In these cases, we exclude the data on the "bad" side of the detector from our analysis; about 25% of the data is on the noisy side of the detector. Additionally, we remove the 2011 September 7 observations from our data set, since the 10 minute total integration time is not sufficient for PCA.

As in Piskorz et al. (2016), all but about 0.1% of the variance in each night's data set is encapsulated by the first principal component. The following results are roughly consistent for data sets with more than the first principal component removed. As discussed in Section 4.4 and shown in Figure 7, the expected photometric contrast α_phot at the observed wavelengths is ∼10⁻⁶. Based on the percent variance removed by each principal component, we determine that deletion of a signal of this magnitude requires the removal of upward of the first 15 principal components from our data. In the analysis that follows, our data set has the first five principal components removed, leaving the stellar and planetary signals intact.

Our PHOENIX stellar model is interpolated between the spectral grid points presented in Husser et al. (2013) for the effective temperature T_eff, surface gravity log g, and metallicity [Fe/H] values listed for ups And A in Table 1. We rotationally broaden this model assuming a stellar rotation rate of 9.62 km s⁻¹ (Valenti & Fischer 2005) and limb darkening coefficient of 0.29 (Claret 2000). For completeness, we instrumentally broaden the stellar model with the kernel determined in Section 2.2.

We compute a high-resolution (R = 250,000) thermal emission spectrum of ups And b according to the SCARLET framework (Benneke 2015). The thermal structure and equilibrium chemistry of the ups And b model spectrum are dependent upon the expected stellar flux at the location of the planet. The model assumes perfect heat redistribution (perhaps a flawed assumption, see Crossfield et al. 2010 and Section 4.3) and a solar elemental composition (Asplund et al. 2009). The temperature profiles are computed self-consistently for a 1 × solar, C/O = 0.54 atmosphere by iteratively recalculating the radiative-convective equilibrium and atmospheric equilibrium chemistry. We assume an internal heat flux of T_int = 75 K. Our default model in this paper has an inverted temperature structure due to the short-wavelength absorption of TiO and VO. The SCARLET framework includes molecular opacities of ${{\rm{H}}}_{2}{\rm{O}}$, CH₄, NH₃, HCN, CO, CO₂, and TiO (ExoMol database by Tennyson & Yurchenko 2012), molecular opacities of O₂, O₃, 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₂ (HITRAN database by Rothman et al. 2009), absorptions by alkali metals (VALD database by Piskunov et al. 1995), H₂-broadening (Burrows & Volobuyev 2003), and collision-induced broadening from H₂/H₂ and ${{\rm{H}}}_{2}/\mathrm{He}$ collisions (Borysow 2002).

Line positions and amplitudes are critical to obtaining the correct cross-correlation function. We use the line information from ExoMol for ${{\rm{H}}}_{2}{\rm{O}}$ and CH₄. The line lists were computed using ab initio calculations based on quantum mechanics. The line center wavelengths of these calculations are accurate. Line amplitudes are harder to compute in ab initio calculations, but we are using the best state-of-the-art line lists available, which is ExoMol for the temperatures encountered in hot Jupiters. Model spectra are convolved with the instrumental profile from Section 2.2 before the cross-correlation analysis.

We use the TODCOR algorithm (Zucker & Mazeh 1994) to cross-correlate each order of data for each epoch with the stellar and planet models, yielding a two-dimensional array of cross-correlation values for different stellar and planetary velocity shifts.

As in Piskorz et al. (2016), at this step, we eliminate the K_r and K_l band orders ranging from 2.3 to 2.4 μm from the analysis since there is high correlation between the stellar and planetary models themselves at these wavelengths. This means we remove any signal from carbon monoxide, and the dominant molecule in the planetary model in the remaining wavelengths is water vapor.

Following Lockwood et al. (2014), for each epoch of observations, we combine the correlation function for each order and produce nightly stellar and planetary maximum likelihood curves, a few of which are shown in Figure 4. For every epoch, we are able to confirm the expected velocity of the star

$$\begin{eqnarray}&&{v}_{\mathrm{pri}}={v}_{\mathrm{sys}}-{v}_{\mathrm{bary}}\end{eqnarray} \tag{ 1 }$$

(where v_sys is the systemic velocity of ups And A and v_bary is the barycentric velocity of the Earth at the time of observation) as is the shown by the strong peaks in the panels on the left-hand side of Figure 4. We suspect that the significant off-peak correlation signature in the primary velocity curve for the K_l band data implies that we were too aggressive in our clipping and that we have scratched the noise limit of our data (see Section 4.4).

Zoom In Zoom Out Reset image size

The right column of Figure 4 shows the maximum likelihood curves for shifts in the planetary velocity. Two aspects are notable. First, the likelihood variations of the K band data are an order of magnitude smaller than those of the L data, indicative of the small signals present in the K band data. Second, there are many peaks and troughs in the planetary maximum likelihood curves. Therefore, determining the line-of-sight velocity of the planet is not straightforward. Only one peak in each maximum likelihood curve represents the real planetary velocity for a given epoch; the other peaks are chance correlations with the repeating structure in the planetary model. The multi-epoch data are critical in breaking this degeneracy.

We use the cross-correlation functions for the planetary velocity shift v_sec at each epoch to determine the most likely value of the line-of-sight Keplerian velocity K_P. For the sake of completeness, we use the equation for orbital velocity, which considers eccentricity, even though the eccentricity of ups And b is very nearly zero. As a result of this near-zero eccentricity, the mean anomalies M of our observations are essentially the same as the true anomalies f. The velocity v_sec of the planet a function of its true anomaly f is

$$\begin{eqnarray}&&{v}_{\sec }(f)={K}_{{\rm{p}}}(\cos (f+\omega )+e\cos \omega )+{v}_{\mathrm{pri}},\end{eqnarray} \tag{ 2 }$$

where K_P is the planet's orbital velocity, ω is the longitude of periastron measured from the ascending node, and e is the eccentricity of the orbit. We test orbital velocities from −150 to 150 km s⁻¹ in steps of 1 km s⁻¹ and thus test a variety of planet masses and orbital inclinations. This results in a plot of maximum log likelihood versus the planet's orbital velocity (first column of Figure 5).

Zoom In Zoom Out Reset image size

Six L band cross-correlation functions similar to that in the upper right panel of Figure 4 are combined to produce the likelihood curve in the upper left panel of Figure 5 when combined with equal weighting. The single peak in K_P is at 55 ± 3 km s⁻¹. The error bars reported here are the 3σ error on the mean value of a Gaussian curve fit to the maximum likelihood peak with equal weighting to the points on the maximum likelihood curve. The error bars are not the full-width at half-maximum of the fitted Gaussian. We more robustly calculate the weighting of the points on the maximum likelihood curve and the error bars and significance of the K_P measurement based on the full 11 nights of data later in this section. Three K_r band cross-correlation functions similar to those in the middle right panel of Figure 4 produce the likelihood curve in the second row of the first column of Figure 5 and shows a peak at K_P = 53 ± 3 km s⁻¹. Finally, three K_l band cross-correlation functions similar to that in the bottom right panel of Figure 4 produce the likelihood curve in the third row of the first column of Figure 5 and shows a peak at K_P = 58 ± 3 km s⁻¹.

The combination of all nights of data is shown in the bottom left panel of Figure 5 and gives K_P = 55 km s⁻¹. We use this value of K_P to calculate the expected v_sec for each epoch of observation and note this value as a vertical line on the curves in the right column of Figure 4. For most cases (especially in L and K_r bands), the expected v_sec corresponds to a local maximum in likelihood. We also use K_P = 55 km s⁻¹ to calculate the secondary velocities plotted in Figure 2.

Given the full suite of data, we calculate the error bars of each point of the maximum likelihood curve using jackknife sampling. We remove one night of data from the sample at a time and recalculate the maximum likelihood curve. The error on each point is proportional to the standard deviation of the 12 resulting maximum likelihood curves. These errors are shown in the bottom left panel of Figure 5. These errors are an estimate only. For a Gaussian fit to the peak at 55 km s⁻¹, the reduced chi-squared value (chi-squared divided by the number of degrees of freedom) is 0.15, suggesting that the error bars are likely an overestimate. These large error bars are driven by a high variance in the jackknife samples. The Gaussian fit also gives error bars on the ultimate K_P measurement: 55 ± 9 km s⁻¹.

To determine the significance of this detection, we use the jackknifed error bars to fit a Gaussian (above) and a straight line and we compare the likelihoods of the fits with the Bayes factor B. Here, the Gaussian fit corresponds to the presence of a planetary signal and the linear fit corresponds to the lack thereof. The Bayes factor B is the ratio of the likelihood of two competing models (Kass & Raftery 1995). If 2lnB is greater than 10, then the model is very strongly preferred.

For the Gaussian fit compared to the linear fit, the value of 2lnB is 10.5, indicating that the signal at 55 km s⁻¹ is stronger than a straight line at about 3.7σ. Therefore, the line-of-sight orbital velocity of ups And b is 55 ± 9 km s⁻¹. Using the indicative mass of ups And b and the law of conservation of momentum, we calculate that the true mass of ups And b is ${1.7}_{-0.24}^{+0.33}$ M_J, and the orbital inclination of ups And b is 24° ± 4°.

With SCARLET, we can calculate the contributions of individual molecules (H₂O, CO, and CH₄) to the total spectrum to understand the dominant opacity structures. We cross-correlate these molecular planet models with our L, K_r, and K_l band data. Results of these single molecule cross-correlation calculations are shown in the middle column of Figure 5 for each band observed and indicate that the atmospheric opacity of ups And b is dominated by water vapor at the observed wavelengths. The likelihood curves for data correlated with CO- and CH₄-only planetary models show variations at least an order of magnitude smaller than the H₂O-only results. If carbon monoxide or methane are present at these wavelengths, they exist at levels below the detection limit of this data set. (See Section 4.3.) Note that we were forced to remove the CO band at 2.2935 μm from our data set because of the presence of CO features in the stellar spectrum.

Our initial test of the fidelity of the line-of-sight velocity detection at 55 km s⁻¹ is to vary the spectroscopic contrast α_spec. We test α_spec from 10⁻⁷ to 10⁻³ and find that the peak at 55 km s⁻¹ is robust down to 10^−6.5. α_spec is truly the ratio between the depths of the spectral lines, and so could be as low as zero for perfectly isothermal atmospheres.

Analagous to Piskorz et al. (2016), we produce a "shuffled" planetary model by randomly rearranging chunks of the planetary model. Cross-correlating our data with a shuffled model should show no peak near 55 km s⁻¹ if the planet truly exists with a line-of-sight orbital velocity of 55 km s⁻¹. For each band of data, we run this test three times and the results are shown in the right-hand column of Figure 5. The L, K_r, and K_l band detections show minima near 55 km s⁻¹, showing that the planet signal is successfully eliminated.

We use our inclination measurement of 24° ± 4° to compare our detection to the results presented in other works. The spectroscopic technique presented here would be unable to detect the motion of ups And b if i_b < 49 due to the size of a resolution element on NIRSPEC. Our inclination measurement is largely in agreement with previous works. Spitzer brightness measurements indicated i_b > 28° (Crossfield et al. 2010). Newtonian orbital simulations considering the orbital elements of ups And c and d suggested that orbits having i_b < 60° can be stable (McArthur et al. 2010). Analagous post-Newtonian orbital simulations prescribed a "region of stability" for i_b < 40°. Our measurement of i_b = 24° ± 4° lies within the error bars of these ranges.

Many previous works have characterized the ups And A system as on the edge of instability. Here we evaluate our calculation of the inclination of ups And b by running numerical simulations of the system with the Mercury software (Chambers 1999). Mercury is a hybrid-simplectic-Burlisch Stoer algorithm (Chambers 1999). We include the central star ups And A and the three planets ups And b, c, and d, set the time step to one-twentieth of the orbital period of ups And b, and consider general relativity.

Our method of calculating K_P and i_b provides no insight into the longitude of the ascending node of ups And b, Ω_b. As a result, we investigate values of i_b between 22° and 27° in steps of 1° and values of Ω_b between 0° and 360° in steps of 10°. We adjust M_b as is necessary given the value of i_b. All other orbital elements are taken from McArthur et al. (2010). Specifically, the orbital elements used for our simulations are listed in Table 1.

Of our 216 simulations, 122 were stable for more than 100,000 years. These simulations have Ω_b < 100° or Ω_b > 260°. Of these systems, 53 were stable for more than 1 Myr, having Ω_b < 40° and Ω_b > 320°. We extract the 24 simulations having 23^° < i < 25° and run them for 100 Myr. All but two are stable. It seems that for the successful simulations the orbital planes of planets b and d remain roughly aligned. For example, if i_b = 24° and Ω_b = 0°, then the mutual inclination of ups And b and c is about 29° and the mutual inclination of ups And b and d is about 2°. Recall that the mutual inclination of ups And c and d is 29°. Successful simulations tend to have mutual inclinations clustered about these values. In these simulations, the apsides of ups And c and d oscillate as in Chiang & Murray (2002), Barnes et al. (2011), and other works, and the orbital evolution of ups And b is secular (Figure 6). We stress that these simulations are stable not necessarily because of the value of ups And b's inclination, but because of the direction ups And b's inclination vector points over time. Our Mercury simulations provide evidence that stable ups And A systems do indeed exist for the inclination we have measured, and provide insight into the three-dimensional geometry of ups And b's orbit.

Zoom In Zoom Out Reset image size

In our planetary model, the L band opacity is dominated by water vapor. Therefore, our L band detection of ups And b's thermal emission spectrum suggests that radiative transfer in the planet's atmosphere is dominated by water vapor at these wavelengths. In fact, the source of the correlation signal for all wavelengths investigated is water vapor (see the middle column of Figure 5). Based on the analysis of α_spec presented in Section 4.1, the detection of H₂O suggests that its spectroscopic contrast α_spec > 10^−6.5.

We perform a comparison of the cross-correlation results given inverted and non-inverted model spectra. The main differences in the final maximum likelihood curves stem from the different line strengths at a given wavelength for each model. In other words, the differences stem from the optical depths as a function of wavelength. Therefore, the only conclusion we can draw at this time is the atmosphere of ups And b is dominated by water at the probed wavelengths.

Though the K band is typically dominated by CO absorption, the usable K band wavelengths in our data set do not include strong CO absorption. The non-detections of CO and CH₄ suggest that their spectroscopic contrasts are α_spec < 10^−6.5 at these wavelengths.

Our models do not account for cloud cover, atmospheric recirculation, or the differences between dayside and nightside spectra. Crossfield et al. (2010) reported a flux maximum in the Spitzer phase curve of ups And b at 80° before opposition, or at mean anomaly M = 0.4, in our formulation. M = 0.4 is almost directly between the phases of 2016 November 12 and 2016 August 21 observations as diagrammed in Figure 1. Fortuitously, this indicates that even if the planet's flux maximum is shifted from what would traditionally be expected, our measurements are still able to capture dayside emission.

From our raw data sets, we calculate the shot noise per resolution element for each observation (see Table 2). We compare the aggregate shot noise values to the expected photometric signal from the planet for each order observed, using the stellar and planet models described in Sections 3.1 and 3.2. As Figure 7 suggests, we easily achieve the required S/N to detect the planet with six nights of L band observations, but only marginally achieve that required with three nights of K_l and K_r observations. In fact, we achieve slightly better shot noise for K_r than for K_l, a possible reason for the stronger detection of the planet signal here (Figures 4 and 5).

Zoom In Zoom Out Reset image size

In this suite of observations, the K_l band data sets are equivalent to the K band data presented for HD 88133 in Piskorz et al. (2016). With four nights of K band data, Piskorz et al. (2016) were able to detect the signal from HD 88133 b, though not as clearly as in the six nights of L band data. This points to the general trend that, with NIRSPEC at Keck, L band observations may be more amenable to direct detection of exoplanet atmospheres than those in the K band. For hot Jupiters, the increase in the thermal background from K to L band is more than compensated for by the significant increase in planet flux relative to the star. In other words, though the increment of detection limit achieved per unit integration time is higher in the K band than in the L band, the star–planet contrast near 2 μm may be too small for a bona fide planet detection with our data. For this cross-correlation method, the superiority of L band observations over K band observations is a demonstration of the theoretical results presented in de Kok et al. (2014).