Discovery of Water at High Spectral Resolution in the Atmosphere of 51 Peg b

We observed the bright star 51 Peg (G2.5V, V = 5.46 mag, K = 3.91 mag) for 3.7 hr during the night beginning 2010 October 21, using CRIRES⁸ (Kaeufl et al. 2004) mounted at Nasmyth A at the VLT (8.2 m UT1/Antu), Cerro Paranal, Chile. The observations were collected as part of the ESO large program 186.C-0289. The instrument setup consisted of a 0.2 arcsec slit centred on 3236 nm (order 17), in combination with the Multi-Application Curvature Adaptive Optic system (MACAO; Arsenault et al. 2003). The CRIRES infrared detector is comprised of four Aladdin III InSb-arrays, each with 1024 × 512 pixels, and separated by a gap of 280 pixels. The resulting wavelength coverage of the observations was $3.1806\lt \lambda \,(\mu {\rm{m}})\lt 3.2659$ with a resolution of $R\approx$ 100,000 per resolution element (see Figure 1).

Zoom In Zoom Out Reset image size

We observed 51 Peg continuously while its hot Jupiter companion passed through orbital phases $0.55\lesssim \phi \lesssim 0.58$, corresponding to an expected change in the planet's radial velocity of ${\rm{\Delta }}{\mathrm{RV}}_{{\rm{P}}}=-23$ km s⁻¹ (15 pixels on the CRIRES detectors). In total, we obtained 42 spectra, with the first 20 spectra each consisting of two sets of 5 × 20 s exposures, and the remainder each consisting of two sets of 5 × 30 s exposures. The increase in the exposure time was aimed at maintaining a constant signal-to-noise ratio (S/N) in the continuum of the observed stellar spectra after a sudden and significant deterioration of the seeing (increasing from 0.75 to 1.4 arcsec between one set of frames, see Section 2.3). To enable accurate sky-background subtraction, the telescope was nodded along the slit by 10 arcsec between each set of exposures in a classic ABBA sequence, with each of the final 42 extracted spectra consisting of an AB or BA pair. A standard set of CRIRES calibration frames was taken the following morning. Later in this paper, we will compare our results to those of Brogi et al. (2013), who observed 51 Peg b at $2.3\,\mu {\rm{m}}$ at the same spectral resolution with CRIRES/VLT on dates either side of these observations, including 2010 October 16, 17, and 25.

Throughout the data reduction, the four CRIRES detectors were treated independently and separately. We used the CRIRES esorex pipeline (v2.2.1) to first process the observed 2D images, including nonlinearity and bad pixel corrections, flat-fielding, background subtraction and combination of the nodded exposures, and finally the optimal extraction of the 1D spectra. The 42 extracted spectra were stored as four 1024 × 42 matrices. The matrix for detector 1 is shown in the upper panel of Figure 2. The x-axis corresponds to pixel number (i.e., wavelength channel) and the y-axis denotes the frame number (i.e., orbital phase or time). Remaining singular bad pixels, bad regions, and bad columns in these matrices were identified iteratively by eye and replaced by spline interpolation values from their horizontal neighboring pixels. There was a total of 0.4–1.0 per cent bad pixels in each matrix, with detector 4 requiring the most corrections.

Zoom In Zoom Out Reset image size

A gradual drift occurs in the position of the spectrum in the dispersion direction on the detectors over the course of the observations. To correct this, we apply a global shift to each spectrum on each detector using a spline interpolation to align it to the telluric features of the spectrum with the highest S/N. The shifts were determined by cross-correlating the spectrum in question with the highest S/N spectrum using iraf.fxcor. Detector 1 required the largest corrections, with the last spectrum deviating from the first by 0.7 pixels (the equivalent of 1 km s⁻¹ or 0.01 Å). Brogi et al. (2013) note that their $2.3\,\mu {\rm{m}}$ observations of 51 Peg b also experience a similar drift, which correlates with the temperatures of the instrument pre-optics system, grating, and its stabilizer. For small fluctuations ($\lt 0.05$ K), the drift did not exceed 0.5 pixels, but a 1.5 K change in these temperatures resulted in much larger drifts (1.5 pixels) for their observations on 2010 October 25, which resulted in a non-detection of the 51 Peg b signal. The drift in our $3.5\,\mu {\rm{m}}$ observations does not correlate with these instrumental temperatures, which remained stable throughout the night. Instead, the drift correlates with the ambient temperature of the telescope dome and the primary mirror temperature, which both cooled by 2 K over the course of the observations. However, the drift of our $3.5\,\mu {\rm{m}}$ spectra is comparatively small, thus we do not expect the alignment correction to significantly affect our subsequent analysis.

Finally, a common wavelength solution per detector was calculated using a synthetic telluric transmission spectrum (see the second panel of Figure 1) from ATRAN⁹ (Lord 1992) to identify the wavelengths of the telluric features in the highest S/N spectrum. Line positions were identified using iraf.identify, and fitted with a third-order Chebyshev polynomial to obtain the wavelength solution. This replaced the default solution from the CRIRES pipeline from which it differed by up to 1.9 Å. The wavelength solutions were not linearized and thus retained the pixel-spacing information. Detector 2 contains significant telluric contamination (see Figures 1 and 2) such that no useful planetary signal can be extracted. Following Birkby et al. (2013), we therefore discarded detector 2 and exclude it from all further analyses in this paper.

Telluric absorption from the Earth's atmosphere is the dominant spectral feature in our observed spectra (see Figure 1), while the Doppler-shifted features of 51 Peg b are expected at the ${10}^{-3}\mbox{--}{10}^{-4}$ level with respect to the stellar continuum. Thus, we needed to remove the telluric contamination. In a previous analysis of high-resolution spectra of 51 Peg b, Brogi et al. (2013) included an additional step to remove stellar lines before removing the telluric features. However, they studied the $2.3\,\mu {\rm{m}}$ region, which contains multiple strong CO lines from the Sun-like host star. In the $3.2\,\mu {\rm{m}}$ region under consideration here, a comparison with a proxy solar model spectrum for 51 Peg from Robert Kurucz's stellar model database¹⁰ at $R\,=$ 100,000 indicates that there are few strong absorption lines from the host star, and that they mostly fall on gaps between detectors, or on the discarded detector 2 (see the third panel of Figure 1). Consequently, we do not perform any pre-removal of the stellar lines in the $3.2\,\mu {\rm{m}}$ data set.

The removal of the telluric features in our spectra was achieved using our implementation of Sysrem, which is an algorithm based on principle component analysis but also allows for unequal error bars per data point (Tamuz et al. 2005; Mazeh et al. 2007). It is commonly used in ground-based wide-field transit surveys to correct systematic effects common to all light curves (e.g., SuperWASP; Collier Cameron et al. 2006). Each wavelength channel (i.e., pixel column) in the matrix shown in panel B of Figure 2 was treated as a "light curve", where the errors per data point are the quadrature sum of the Poisson noise and the error from the optimal extraction of the spectrum at a given pixel. Each column had its mean subtracted before being passed through Sysrem, and regions of saturated telluric features, which contain essentially no flux, were also masked. These masked regions are marked by the vertical solid gray regions in Figure 2. Sysrem then searched for common modes between the 1024 light curves per matrix, such as variation with airmass, and subtracted them, resulting in the removal of the quasi-static telluric and stellar lines, leaving only the Doppler-shifting planet spectrum in each spectrum plus noise. However, in practice, once the dominant telluric and stellar spectral features are removed, Sysrem will begin to remove the planet features too. This is because the sub-pixel shift of the planet spectrum between frames creates a small but detectable common mode between adjacent columns. Thus, we must determine when to halt the Sysrem algorithm before it removes the planet signal. To do this, we injected the best-matching model planet spectrum proposed by Brogi et al. (2013) based on observations at $2.3\,\mu {\rm{m}}$ at their measured planet velocity and ephemeris (${V}_{\mathrm{sys}}=-33$ km s⁻¹, ${K}_{{\rm{P}}}\ =134$ km s⁻¹, with a phase shift of ${\rm{\Delta }}\phi =0.0095$), and iterated Sysrem ten times. The model was injected at a nominal strength of 1. At each iteration, we used the cross-correlation method described in Sections 3.2 and 3.3 to determine the significance of the detection at the injected velocity. The results of this analysis are shown in Figure 3 for each detector. We find that Sysrem begins to remove the planet signal after only one iteration on detector 3, but that this does not occur until after two and three iterations on detectors 1 and 4, respectively. The difference in iterations required per chip may be partly due to a known "odd–even" effect,¹¹ which only affects detectors 1 and 4. It is caused by variations in the gain between neighboring columns along the spectral direction. This odd–even effect has been seen previously in similar analyses of high-resolution spectroscopy of exoplanets from CRIRES/VLT (Brogi et al. 2012, 2013, 2014, 2016a; Birkby et al. 2013).

Zoom In Zoom Out Reset image size

Guided by the results of executing Sysrem on the injected signal as described above, we adopted Sysrem iterations of two, one, and three for detectors 1, 3, and 4, respectively, for the analysis of the observed data (see Section 3.2). The standard deviation of the final residuals in panels D of Figure 2 are 0.0050, 0.038 0.0082, 0.0040, for detectors 1, 2, 3, and 4, respectively. The trends removed from each detector during each Sysrem iteration are shown in Figure 4. Possible physical causes of these trends are shown in Figure 5, and these are discussed in further detail in Section 3.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final two steps in the telluric removal process were: (i) the application of a high-pass filter with a 64 pixel width smoothing function, which removes a heavily smoothed version of each residual spectrum from itself to filter out low-order variation across the matrix, and then (ii) each column is divided by its variance to account for variation in S/N as a function of wavelength. The final product of this process for detector 1 is shown in panel (D) of Figure 2. For illustrative purposes, the bottom panel of Figure 2 shows how the final matrix would appear if a model planet was injected at the planet velocity but at $100\times$ nominal value before running Sysrem. Many individual lines from the planet spectrum are clearly visible as they blueshift across the matrix. With the tellurics and stellar continuum removed from each detector, we could proceed to extract the observed planetary signal contained within the noise of the residual spectra.

The molecular templates are parameterized by a grid of atmospheric temperature–pressure (T − P) profiles and trace gas abundances (i.e., volume mixing ratios; VMR). To generate the model spectra, we employed the same radiative transfer code of de Kok et al. (2013), performing line-by-line calculations, including H₂–H₂ collision-induced absorption along with absorption from the trace gases, which is assumed to follow a Voigt line shape. We used line data from HITEMP 2010 (Rothman et al. 2010) to create the H₂O and CO₂ models, and HITRAN 2008 (Rothman et al. 2009) for CH₄. The model atmospheres are clear (i.e., cloud-free) with uniformly mixed gases. The T − P structure follows a relatively simple profile. Deep in the atmosphere, at pressures p₁ and higher, we assume a uniform T − P profile at a fixed temperature t₁. Between pressures p₁ and p₂, we assume a constant lapse rate (i.e., a constant rate of change of temperature with log pressure). At altitudes higher (and pressures lower) than p₂ we again assume a uniform T − P profile at fixed temperature t₂. The pressure p₁ took values of (1, 0.1, 0.01) bar, and p₂ was varied using values of ($1\times {10}^{-3},1\times {10}^{-4},1\times {10}^{-5}$) bar. The basal temperatures were guided by the effective temperature of the planet assuming external heating (see Equation (1) of López-Morales 2007), a low Bond albedo (${A}_{B}\lt 0.5$), and considering the full range of heat circulation from instantaneous reradiation to full advection. Thus, t₁ took values of (1000, 1250, 1500) K, while t₂ was varied using values of (500, 1500) K. Note that certain combinations result in inverted T − P profiles, where the temperature increases with increasing altitude (decreasing pressure). These spectra have features in emission, rather than absorption. The gas abundance volume mixing ratios took values appropriate for hot Jupiters over our considered temperature range (e.g., Madhusudhan 2012), including ${10}^{-4.5}$, 10⁻⁴, or ${10}^{-3.5}$ for water, and 10⁻⁷, 10⁻⁵, or 10⁻³ for CO₂ and CH₄.

Before cross-correlating the residuals with the molecular template grid, we convolved the models to the spectral resolution of CRIRES, and subtracted their baseline level. Note that the telluric removal process in Section 2.3 has also removed the continuum information in the observed planet spectrum, such that our analysis is only sensitive to the relative, not absolute, depth of the spectral features with respect to the stellar continuum.

The cross-correlation analysis was performed for planet radial velocities in the range $-249\leqslant {\mathrm{RV}}_{{\rm{P}}}(\mathrm{km}\,{{\rm{s}}}^{-1})\leqslant 249$ in intervals of 1.5 km s⁻¹, interpolating the convolved grid of molecular templates onto the Doppler-shifted wavelengths. The interval size is set by the velocity resolution of the CRIRES pixels. The cross-correlation functions (CCFs) were determined separately for each residual spectrum on each detector, and then summed equally with their corresponding CCF on the other detectors, resulting in a single summed CCF matrix of dimension 333 × 42. The matrix created by the best-matching template is shown in upper left panel of Figure 6, where the template is a water-only model with the following parameters: ${t}_{1}=1500$ K, ${t}_{2}=500$ K, ${p}_{1}=0.1$ bar, ${p}_{2}=1\times {10}^{-5}$ bar, and a water volume mixing ratio of ${\mathrm{VMR}}_{{\rm{H}}2{\rm{O}}}={10}^{-4}$. These parameters differ to the best-matching model reported by Brogi et al. (2013) for molecular absoprtion at $2.3\,\mu {\rm{m}}$; however, see Sections 3.3 and 4 for further discussion on other models in the grid that produce signals within $1\sigma$ of this result. We refer the reader to Section 3.4 for a more detailed discussion of Figure 6.

Zoom In Zoom Out Reset image size

We note that in our first attempt to perform the cross-correlation analysis of these data, we used the ephemeris and orbital solution for 51 Peg b in Butler et al. (2006), and refrained from using the phase shift invoked by Brogi et al. (2013) to match the planet signal to the systemic velocity of the host star, even though their ${{\rm{\Delta }}}_{\phi }=0.0095$ phase shift was within the $1\sigma$ uncertainty range ${\rm{\Delta }}\phi =\pm 0.012$ of the Butler et al. (2006) ephemeris. However, we also found that this resulted in the strongest cross-correlation signal being offset from the known systemic velocity of the host star by −9 ± 2 km s⁻¹, corresponding to a phase shift of ${\rm{\Delta }}\phi =0.011$. This is still within the uncertainty of the original Butler et al. (2006) orbital solution. However, since the discovery of 51 Peg b, the RV of its host star has been monitored sporadically throughout the decades, hence we endeavoured to measure the most up-to-date ephemeris for the planets orbital solution in the hope that this would negate the need for the phase shift in our analysis.

3.2.1. A Refined Orbital Solution for 51 Peg b

To ensure we had its most up-to-date orbital solution, we compiled an extensive repository of literature and archival radial velocity measurements of the star 51 Peg from multiple observatories. These data are given in Table 1 and span observing dates from 1994 September 15 to 2014 July 9, resulting in 639 RV measurements over 20 years. The table includes the discovery measurements from the ELODIE spectrograph at Observatoire Haute Provence (Mayor & Queloz 1995) and subsequent additional monitoring. We took these RV measurements from the Naef et al. (2004) compilation. We also included the legacy data set from Lick Observatory observed with the Hamilton spectrograph, taking measurements from the self-consistent reprocessing of all the Lick spectra presented by Fischer et al. (2014). Finally, we included more recent additional monitoring from HIRES at the Keck Observatory (Howard & Fulton 2016), and extracted RVs from observations with HARPS at the ESO-3.6 m telescope in 2013 (ESO program ID 091.C-0271, PI: Santos). The reduced HARPS spectra were obtained from the ESO Science Archive,¹² and the RVs, their errors, and timestamps were obtained from the headers of the CCF data product files, using keywords drs ccf rvc, drs dvrms, and drs bjd, respectively. We were careful to note the format of the timestamps reported for all our data sets, which vary between JD, HJD, and BJD, due to different conventions being adopted over time. The timestamps in Table 1 were all converted to ${\mathrm{BJD}}_{\mathrm{TDB}}$ (Barycentric Dynamical Time).

Table 1.  Stellar Radial Velocity Measurements of 51 Peg from Multiple Observatories

${\mathrm{BJD}}_{\mathrm{TDB}}$	RV (m s−1)	${\sigma }_{\mathrm{RV}}$ (m s−1)	Data set
2449610.532755	−33258.0	9.0	ELODIE
2449612.471656	−33225.0	9.0	ELODIE
2449655.311263	−33272.0	7.0	ELODIE

Note. The Lick8 RVs reported in this table include a $+13.1$ m s⁻¹ velocity offset correction to the RVs extracted from Vizier (http://vizier.cfa.harvard.edu/viz-bin/VizieR-3?-source=J/ApJS/210/5/table2), to account for the instrumental systematic reported in Fischer et al. (2014). The following additional offsets, determined from a circular orbit fit to each data set using Exofast, can be applied to place all of the RV measurements onto the same zero-point: ${\gamma }_{\mathrm{Lick}13}=-21.70$ m s⁻¹, ${\gamma }_{\mathrm{Lick}8}=-4.52$ m s⁻¹, ${\gamma }_{\mathrm{Lick}6}\ =-14.64$ m s⁻¹, ${\gamma }_{\mathrm{ELODIE}}=+33251.59$ m s⁻¹, ${\gamma }_{\mathrm{HIRES}}=+2.24$ m s⁻¹, ${\gamma }_{\mathrm{HARPS}}=\,+33152.54$ m s⁻¹.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We used Exofast ¹³ (Eastman et al. 2013) to model the orbital components constrained by each data set. To find the best-fitting model to the radial velocity data, Exofast employs a nonlinear solver (amoeba), which uses a downhill simplex to explore the parameter space that minimizes the ${\chi }^{2}$ of the orbital solution. In order to negate any systematic underestimate of the uncertainties on the RV data, Exofast rescales the RV uncertainties by a constant multiplicative factor, such that the reduced ${\chi }^{2}$ of the best fit is unity (${\chi }_{\nu }^{2}=1$). Consequently, a poor fit would be reflected by larger uncertainties on the derived parameters. Exofast does not include an additive jitter term to the RV solution, as previous work with Exofast found no statistically significant difference between the two approaches to uncertainty scaling (e.g., Lee et al. 2011). Once the best fit is found, Exofast then executes a differential evolution Markov chain Monte Carlo method to obtain the uncertainties on the derived orbital elements. These algorithms are explained in detail in Eastman et al. (2013). The code requires priors for the stellar effective temperature (${T}_{\mathrm{eff}}=5787\pm 233$ K), metallicity ($[\mathrm{Fe}/{\rm{H}}]=0.200\pm 0.030$), and surface gravity ($\mathrm{log}({g}_{* })=4.449\pm 0.060$), which we supplied based on Valenti & Fischer (2005). We also gave the logarithm of the period (${\mathrm{log}}_{10}(P)={\mathrm{log}}_{10}(4.230785\pm 0.000036)$) as a prior based on Butler et al. (2006), and restricted the period range to 4–5 days.

Combining relative RV measurements from different observatories is hindered by velocity offsets in part due to instrumental systematics. We consider our repository to consist of six separate data sets: HARPS, HIRES, ELODIE, plus three Lick data sets (Lick13, Lick8, and Lick6). The numbers in the Lick data set names denote the dewar associated with each upgrade to the Hamilton spectrograph, which introduced different velocity offsets (see Fischer et al. 2014 for further information). Our first step was to model the data sets independently to determine if they supported an eccentric orbit. We ran two models, one with and one without a free variable for a long-term linear trend in the RVs. We excluded the HARPS data set in this assessment, as it has poor phase coverage, being clustered within $\phi =\pm 0.1$ of superior and inferior conjunction and thus lacking in strong constraint on the point of maximum absolute radial velocity. We ran the Lucy–Sweeney test to determine if the derived small eccentricities in our orbital solutions of the remaining data sets were significant (Lucy & Sweeney 1971). The probability of small eccentricity values arising by chance were $\gt 5 \%$ in all cases. Thus, we fixed the eccentricity to zero (circular orbits) in our subsequent modeling with Exofast, which allowed us to determine the velocity offset values for each data set. These offsets are listed in the notes of Table 1, and were subsequently used to place all the RVs onto the same zero-point velocity.

Prompted by the report of a 1.64 m s⁻¹ yr⁻¹ trend in the RVs of 51 Peg by Butler et al. (2006) and scatter in the discovery RVs reported by Mayor & Queloz (1995), we assessed the significance of long-term linear trends in our circular orbit solutions. We found that the earliest Lick RVs (Lick13), spanning 791 days, supported a $-{1.64}_{-1.10}^{+1.17}$ m s⁻¹ yr⁻¹, in agreement with Butler et al. (2006). However, the remaining Lick data sets, spanning 1175 days and 3354 days, respectively, returned non-significant linear trends of $-{0.58}_{-0.88}^{+0.84}$ m s⁻¹ yr⁻¹ and ${0.029}_{-0.29}^{+0.28}$ m s⁻¹ yr⁻¹, respectively. The ELODIE observations, including concurrent RVs with the Lick13 measurements, however, suggested a much smaller trend ($-{0.15}_{-0.40}^{+0.37}$ m s⁻¹ yr⁻¹), in better agreement with the more recent HIRES RVs ($-{0.33}_{-0.19}+0.19$ m s⁻¹ yr⁻¹). It is possible that the ensemble RV data for 51 Peg probe the turnover point of an RV curve for a long-period companion. However, it is beyond the scope of this paper to search for such a companion. To obtain our final orbital solution for 51 Peg b, we place all the RVs onto the zero-point described above, and model the system as a circular orbit, with no long-term linear trend. Our orbital solution is given in Table 2 and shown in Figure 7. The correlation between the parameters in the fit are shown in the Appendix along with their covariance values. We note that the additional RV scatter from any long-term companion does not cause prohibitively large error bars in the ephemeris for our later analysis.

Zoom In Zoom Out Reset image size

Table 2.  Updated Orbital Solution and Planet Properties for 51 Peg b

Parameter	Units	Value
Stellar Parameters:
M⋆	Mass (${M}_{\odot }$)	1.100 ± 0.066
R⋆	Radius (${R}_{\odot }$)	${1.020}_{-0.079}^{+0.084}$
L⋆	Luminosity (${L}_{\odot }$)	${1.05}_{-0.24}^{+0.32}$
${\rho }_{\star }$	Density (cgs)	${1.46}_{-0.27}^{+0.34}$
$\mathrm{log}({g}_{\star })$	Surface gravity (cgs)	${4.452}_{-0.059}^{+0.061}$
${T}_{\mathrm{eff}}$	Effective temperature (K)	5790 ± 230
$[\mathrm{Fe}/{\rm{H}}]$	Metallicity	0.198 ± 0.029
Planetary Parameters:
P	Period (days)	${4.2307869}_{-0.0000046}^{+0.0000045}$
a	Semimajor axis (au)	${0.0528}_{-0.0011}^{+0.0010}$
Teqa	Equilibrium Temperature (K)	${1226}_{-69}^{+72}$
$\langle F\rangle$	Incident flux (109 erg s−1 cm−2)	${0.51}_{-0.11}^{+0.13}$
RV Parameters:
K⋆	RV semi-amplitude (m s−1)	54.93 ± 0.18
${M}_{{\rm{P}}}\sin i$	Minimum mass ($\,{M}_{{\rm{J}}}$)	0.466 ± 0.019
${M}_{{\rm{P}}}/{M}_{\star }$	Mass ratio	${0.0004043}_{-0.0000080}^{+0.0000085}$
${\gamma }_{\mathrm{zp}}$	Residual zero-point offsetb (m s−1)	${0.032}_{-0.077}^{+0.076}$
TC	Time of conjunction (${\mathrm{BJD}}_{\mathrm{TDB}}$)	2456326.9314 ± 0.0010
${P}_{T,G}$	A priori transit probability	${0.092}_{-0.013}^{+0.014}$
Mp	Planet mass ($\,{M}_{{\rm{J}}}$)	${0.476}_{-0.031}^{+0.032}$
i	Inclination (°)	70–90

Notes. The stellar RV parameters were derived with Exofast, based on priors for $\mathrm{log}(g)$, ${T}_{\mathrm{eff}}$, [Fe/H], and ${\mathrm{log}}_{10}(P)$ whose values are noted in the main text of Section 3.2.1. The eccentricity was fixed to e = 0 and we did not allow for a long-term linear trend. The true planet mass and its range of inclination were derived from the CRIRES spectra as detailed in Section 4.1.

^aT_eq was calculated assuming zero Bond albedo, ${A}_{B}=0,$ and perfect redistribution of incident flux following Hansen & Barman (2007). ^bThis is the residual scatter around the zero-point determined from the independent orbital fits. It is consistent with zero within the $1\sigma$ error bars, and we note that the velocity resolution of the CRIRES pixels (1.5 km s⁻¹) is insensitive to this small discrepancy. We adopt a literature systemic velocity of −33.25 km s⁻¹ (Brogi et al. 2013).

Download table as:  ASCIITypeset image

For reference, we also show in the top panel of Figure 7 photometric monitoring data obtained with Hipparcos (ESA 1997) which were also extracted from the NASA Exoplanet Archive. The data were obtained between 1989 November and 1992 December. The plotted data only include measurements with a quality flag of zero, and are zoomed such that variations at the few percent level could be clearly seen. There are no transit events within the scatter of the data when folded on the orbital period of 51 Peg b. This is entirely consistent with more recent ground- and space-based photometric monitoring of 51 Peg, which also report no transit events corresponding to Earth-size planets or larger at the orbital period of 51 Peg b (Guinan et al. 1995; Mayor et al. 1995; Henry et al. 1997; Walker et al. 2006). This places an upper limit on the orbital inclination angle, which is discussed in Section 4.1.

With the orbital solution updated, we re-ran the cross-correlation analysis using our new ephemeris for 51 Peg b. This resulted in the strongest cross-correlation signal appearing $-13.3\pm 2$ km s⁻¹ from the known systemic velocity of the host star i.e., even further offset than when using the Butler et al. (2006) orbital calculation. In order to match the signal to the known systemic velocity, we must invoke a phase shift of ${\rm{\Delta }}\phi =0.1$, or equivalently ${\rm{\Delta }}{T}_{C}=0.07$, which is significantly larger than our $1\sigma$ error on T_C. One therefore might conclude that the strongest cross-correlation signal is not associated with 51 Peg b. However, given that we see a strong signal at almost identical offsets across multiple data sets (i.e., those reported here and at $2.3\,\mu {\rm{m}}$ by Brogi et al. 2013), which targeted different wavelength regimes and different molecules, and underwent different data processing techniques, the case of association with the star remains plausible. While it is possible that additional long-period companions in the 51 Peg system could be affecting the orbit of 51 Peg b, it seems unlikely that such a companion could induce the magnitude of shift in phase we have measured. We instead note that constraints on certain parameters in the orbital solution using solely the stellar RVs are not strong, especially for the argument of periastron (${\omega }_{\star }$). Using our orbital solution in Table 2, but allowing a small eccentricity (e = 0.001), which cannot be constrained by the existing stellar RV data, we find that an offset in ${\omega }_{\star }$ of just 6^◦ from the standard definition (${\omega }_{\star }=90^\circ$) aligns the strongest cross-correlation signal with the known systemic velocity of the host star. This is considerably smaller than the typical ${\omega }_{\star }$ uncertainties reported in the literature from stellar RV measurements. For example, using the Exoplanet Data Explorer¹⁴ (Han et al. 2014), we find that for similar, non-transiting hot Jupiter systems ($P\lt 10$ days) with small orbital eccentricity ($e\lt 0.1$), the smallest reported error on ${\omega }_{\star }$ is $\pm 11^\circ$, and the rest are all larger than 36^◦. Given that such small changes in the orbital solution can result in alignment of the star and planet systemic velocities, we conclude that the most likely scenario arising from our data is that strongest cross-correlation signal is associated with the 51 Peg system, and thus of planetary nature. However, for simplicity, we adopt a phase shift of ${\rm{\Delta }}\phi =0.1$ in the circular orbital solution, as noted above, to align the signal, rather than modify the eccentricity and ${\omega }_{\star }$.

We further note that recent studies have highlighted that the cross-correlation of water models from different molecular line list databases can result in a velocity offset due to mismatching lines (Brogi et al. 2016b). In the case of 51 Peg b, the velocity shift is seen with both water and CO models, where the latter molecule has a very robust line list. Consequently, we think it unlikely that water line lists are causing the velocity shift we observe in 51 Peg b, but we highlight it as a potential issue for other future studies of exoplanet atmospheres at high spectral resolution.

We further conducted our cross-correlation analysis using other molecular templates containing signatures of CH₄, CO₂ for similar grid parameters, but we found no significant ($\gt 3\sigma$) signal from these molecules that would indicate their presence at the abundances probed by our model grid in the atmosphere of 51 Peg b (see Section 4.2 for a discussion on the possible causes for these non-detections).

In the final CCF matrix the planet signal appears as a dark diagonal trail that, in this case, is blueshifting across the matrix. The trail is a small section of the planet RV curve and its slope determines the planet velocity. For visual purposes only, the bottom left panel of Figure 6 shows the CCF matrix after the best-matching template was injected into the observed spectra, before the removal of telluric contamination, at $7\times$ its nominal value at the detected velocity of the planet. When shifted into the planet's rest frame velocity, the trail becomes vertically aligned, as shown in the middle panels of Figure 6. The size of the shift to vertically align the signal in each CCF is related to the RV semi-amplitude of the planet (K_P). Although Brogi et al. (2013) reported a tentative value of ${K}_{{\rm{P}}}=134.1\pm 1.8$ km s⁻¹ for 51 Peg b, we opted to search for the planet signal over a wide range of K_P values to act as a blind search that would independently confirm this value. We therefore aligned CCFs for a range of K_P values, from $20\leqslant {K}_{{\rm{P}}}(\mathrm{km}\,{{\rm{s}}}^{-1})\leqslant 180$ in steps of 1 km s⁻¹. To determine when the trail was vertically aligned, i.e., in the planet rest frame, thus yielding the value of K_P, we performed a Welch T-test on each aligned CCF matrix. This statistic compared the distribution of a 3 pixel wide sliding column of "in-trail" CCFs values in the aligned matrix to the distribution of those outside it ("out-of-trail"), and determined the probability that they were drawn from the same parent distribution. The sliding of the column allows for the location of ${V}_{\mathrm{sys}}$, as well as K_P. These probabilities were converted into σ-values, using the erf IDL function, and are displayed as a matrix in Figure 8, where ${V}_{\mathrm{sys}}$ denotes the systemic velocity of the central column of the in-trail data. The most discrepant in- and out-of-trail distributions deviated by $5.6\sigma$ and corresponded to ${V}_{\mathrm{sys}}=-33\pm 2$ km s⁻¹, which is coincident with the known systemic velocity of the host star, and at ${K}_{{\rm{P}}}={133}_{-3.5}^{+4.3}$ km s⁻¹, as marked by the black cross in Figure 8. The error bars are the $1\sigma$ marginalized uncertainties, although it is clear they are correlated based on the ellipsoidal shape of the $1\sigma$ error contour shown in white in Figure 8. This is expected, because as we approach the correct K_P for the planet when aligning the CCFs, the detection strength will increase. If the observations had occurred before superior conjunction, the contour would be slanted in the opposite direction. Combining multiple nights of similar data, but spanning different phase ranges, would significantly constrain the $1\sigma$ contour (see e.g., Brogi et al. 2012). Following Brogi et al. (2014), in Figure 8 we also explore negative orbital velocities for the planet. Negative velocity implies a retrograde orbit which we know is not true for 51 Peg b. However, any residual correlated noise that could produce a false positive would interact with the model spectrum in a similar way in both positive and negative velocity space, with a false positive appearing as a mirror image of the positive signal in the ${K}_{{\rm{p}}}\lt 0$ km s⁻¹ plane. The lack of significant negative velocity signals in Figure 8 therefore serves to highlight the robustness of the positive velocity detection.

Zoom In Zoom Out Reset image size

A final step in our analysis was to determine which other model water templates also provided a detection within $1\sigma$ of the peak detection significance of $5.6\sigma$ and thus also adequately describe the observations. Non-inverted models with shallow temperature gradients (i.e., ${t}_{1}=1000$ K and ${t}_{2}=500$ K) returned comparatively lower significance values between $2.4\mbox{--}4.5\sigma$. All remaining non-inverted models in the grid, which span the full range of p₁, p₂, and ${\mathrm{VMR}}_{{\rm{H}}2{\rm{O}}}$ that we explored, were within $1\sigma$ of our peak detection significance. However, we found that models in our grid with a temperature inversion (i.e., ${t}_{2}=1500$ K) were all negatively correlated with the observed planet spectrum. This means that the emission lines in the model temperature inversion spectra correlated negatively with the absorption lines in the observed planetary spectrum. Consequently, we confidently rule out the inverted models in our grid as good descriptors of the observations.

A histogram of the corresponding in- and out-of-trail distributions are shown in Figure 9. Note that there are two in-trail distributions shown. The gray one includes all of the observed spectra, while the other (red) one contains only a subset as described below. The latter was used in our analysis described in Section 3.3.

Zoom In Zoom Out Reset image size

We note that in the model-injected CCFs, the trail fades between frames 18–24 (see bottom-right panel of Figure 6), suggesting a fundamental degradation of the observations, and corresponds to the period of unstable seeing noted in Section 2.3. These CCFs act to scramble the distribution of the in-trail signal. Their inclusion in the T-test results in a $\gt 1\sigma$ decrease in the detection significance ($4.3\sigma$) compared to when they are excluded (see the gray histogram in Figure 9 and note how it shifts back toward the out-of-trail distribution). We have explored multiple reasons for why these spectra are degraded, as shown in Section 2.3. We initially attempted to weight the CCFs by the physical parameters displayed in Figure 5. However, only excessively large weightings gave an improvement, resulting in only a few CCFs dominating the final signal. Given that the injected trail also fades, we conclude that these spectra are not useful. We suspect a combination of poor atmospheric conditions and thus possible slit losses are responsible for the degradation. Indeed, Figure 5 shows that frames with the worse seeing required drift corrections on top of the temperature-induced trend, and that the flux received by the AO system is considerably lower during these times. Consequently, we chose to exclude three frames either side of the seeing spike, with the central frame being the first one with the increased exposure time. We note that this may exclude some planetary signal, thus our detection significance is conservative. The results of Section 3.3 therefore include only 35 of the observed 42 spectra.

The detection of molecular features in the atmosphere of 51 Peg b indicates that the 51 Peg Ab system is a double-lined spectroscopic binary. Our measured K_P of 51 Peg b can be combined with K_⋆ and the system mass ratio determined from precision stellar RV measurements to determine the mass of 51 Peg b, independent of its inclination, via $\tfrac{{K}_{\star }}{{K}_{{\rm{p}}}}=\tfrac{{M}_{{\rm{p}}}}{{M}_{\star }}$. Using the stellar properties given in Table 2, the measured true mass of 51 Peg b from our observations is ${M}_{{\rm{p}}}={0.476}_{-0.031}^{+0.032}{M}_{{\rm{J}}}$, placing it firmly in the planetary mass range, and laying to rest any lingering doubts on the true nature of the first reported exoplanet orbiting a solar-like star. The planet's mass places it at the boundary between Jovian and Neptunian worlds, according to recent work by Chen & Kipping (2016). They found data-driven evidence for a break in the power-law relationship between mass and radius for gaseous worlds at ${M}_{{\rm{P}}}=0.41\pm 0.07{M}_{{\rm{J}}}$. This tipping point can be physically interpreted as the mass at which any further accretion of gas into the outer layers of a Neptunian atmosphere overcomes the barrier for self-compression by gravity, leading to Jovians with slightly smaller radii. Using the mass–radius relations calculated by Chen & Kipping (2016), we predict that the radius of 51 Peg b should be between ${R}_{{\rm{p}}}=1.01-1.32{R}_{{\rm{J}}}$, which includes a $14.6 \%$ dispersion in the mass–radius relation. We note that this is considerably smaller than the ${R}_{{\rm{p}}}=1.9{R}_{{\rm{J}}}$ suggested by recent attempts to detect the reflected light from 51 Peg b at optical wavelengths (Martins et al. 2015), and we discuss this further in Section 4.2.2.

Now that we have measured K_P for 51 Peg b, we can also solve for its orbital inclination using Kepler's third law, such that ${{PK}}_{{\rm{P}}}/2\pi a={K}_{{\rm{P}}}/{V}_{{\rm{P}}}=\sin (i)={0.977}_{0.032}^{0.038}$ (assuming ${K}_{\star }\ll {K}_{{\rm{P}}}$, where V_P is the orbital velocity of the planet). This yields an inclination of $i=78^\circ$, but the $1\sigma$ uncertainties permit values within the range $70\lt i(^\circ )\lt 90$. The geometric probability of 51 Peg b transiting its host star is $P(\mathrm{transit})={R}_{\star }/a=9 \%$; however, the lack of transits in photometric monitoring of 51 Peg b (see Figure 7) place an upper limit on the inclination angle at $i\lt 82\buildrel{\circ}\over{.} 2$, as described by Brogi et al. (2013). This is also consistent with estimates of the inclination of the stellar rotation axis at $70{^\circ }_{-30}^{+11}$ (Simpson et al. 2010), and suggests spin–orbit alignment.

4.1.1. Improving Ephemerides for Non-transiting Planets

The best-matching model to the planet spectra from our grid contained molecular features from water only, with a volume mixing ratio of ${\mathrm{VMR}}_{{{\rm{H}}}_{2}{\rm{O}}}={10}^{-4}$, and a T − P profile that steadily reduced in temperature as altitude increased, from ${t}_{1}=1500$ K at ${p}_{1}=0.1$ bar, to ${t}_{2}=500{\rm{K}}$ at ${p}_{2}=1\times {10}^{-5}$ bar. This model, convolved to the $R=$ 100,000 resolution of CRIRES, is shown in Figure 10, along with a simple schematic of its T − P structure. However, it is important to note that a number of other models give detection strengths within $1\sigma$ of the best-matching model at velocities within the $1\sigma$ uncertainties on K_P and ${V}_{\mathrm{sys}}$. This includes models that span the full pressure range and abundances tested by our grid. The only models in the grid we can rule out are those that include a temperature inversion within our probed pressure range, and those with ${T}_{1}\leqslant 1000$ K. Inversion layers, in which the best-matching templates have emission lines, are yet to be identified with this method of atmospheric characterization (Schwarz et al. 2015).

Zoom In Zoom Out Reset image size

The deepest water lines in the best-matching template have a relative depth of $1.0\times {10}^{-3}$. This reduces to $0.9\times {10}^{-3}$when convolved to the resolution of CRIRES. These contrast ratios were determined by scaling and injecting the negative of the best-matching template into the observed spectra before the telluric removal process, at the detected planet velocity parameters, until the peak detection significance was reduced to zero. The best-matching model shown in Figure 10 required a scaling factor of 3.0 times its nominal value to completely cancel the planet signal. The detection of 51 Peg b via ground-based, high-resolution spectroscopy highlights the method's effectiveness in detecting faint companions among the glare of their much brighter host stars, routinely reaching contrast ratios between ${10}^{-4}\,\mathrm{and}\,{10}^{-3}$ at angular separations of 1–3 ms.

No significant detections ($\gt 3\sigma$) were seen when cross-correlating with models containing carbon dioxide or methane. While it is likely that the abundances of these molecules are too low for these data to constrain (see e.g., Madhusudhan 2012; Moses et al. 2013 for typical abundances of hot Jupiter atmospheres as a function of temperature and C/O ratio), we do not rule out that incorrect line positions in the line list used to create the models are responsible for their non-detection. Hot methane and other line lists are known to have inaccuracies (Hargreaves et al. 2015), particularly at high resolution (Hoeijmakers et al. 2015).

4.2.1. Confirmation of Molecular Absorption at $2.3\,\mu {\rm{m}}$

4.2.2. Other Atmospheric Observations of 51 Peg b

Martins et al. (2015) reported a tentative ($3\sigma$) detection of reflected light from 51 Peg b using high-resolution (R = 115,000) spectra observed with HARPS/ESO-3.6 in 2013, implying a planet-to-star flux ratio on the order of ${F}_{{\rm{p}}}/{F}_{{\rm{s}}}={10}^{-4}$ at optical wavelengths. They calculated a geometric albedo of ${A}_{g}=0.5$ under the assumption of a highly inflated planet at ${R}_{{\rm{P}}}=1.9{R}_{{\rm{J}}}$, while smaller planetary radii (e.g., ${R}_{{\rm{P}}}=1.2{R}_{{\rm{J}}}$) would require very high albedos (${A}_{g}\gt 1$), inconsistent with a Lambertian sphere (i.e., isotropic reflection in all directions) implying a strongly backscattering atmosphere. The CCF derived by Martins et al. (2015) for 51 Peg b was also significantly broadened to a $\mathrm{FWHM}=22.6\pm 3.6$ km s⁻¹, which at face value implies rapid rotation of the planet's atmosphere at the altitudes (pressures) probed by the optical spectra. Assuming that 51 Peg b is tidally locked to its host star, the expected rotation speed of the planet is ${V}_{\mathrm{rot}}=2\pi {R}_{{\rm{P}}}/P=1.2-2.3$ km s⁻¹ for ${R}_{{\rm{P}}}=1-1.9{R}_{{\rm{J}}}$. This is close to or below the resolution of the instrument profile for both HARPS and CRIRES. We show the CCF from the CRIRES observations at the detected planet velocity (${K}_{{\rm{P}}}=133$ km s⁻¹) in Figure 11. The CRIRES CCF matches the cross-correlation of the best-matching template cross-correlated with a version of itself, after being broadened to match the R = 100,000 resolution of CRIRES. Using the measured $\mathrm{FWHM}=5.6$ km s⁻¹ of the CRIRES CCF, combined with the lower limit on the orbital inclination of the system and assuming the spin axes are aligned, we derive an upper limit to the rotational broadening of 51 Peg b of 5.8 km s⁻¹ (${P}_{\mathrm{rot}}\gt 0.9$ days). This result applies to the atmospheric rotation at the altitude and pressure probed by our infrared observations. Rather than invoke dynamical arguments e.g., wind sheer, or an exo-ring system (Santos et al. 2015) for the discrepancy with the HARPS optical measurement of the FWHM, we instead highlight that Martins et al. (2015) urged caution regarding their FWHM measurement, stressing that injected planet spectra and the host star spectrum both resulted in over-broadened CCFs in their data analysis, possibly due to non-Gaussian noise in the data, indicating that the parameters of their CCFs were strongly affected by the noise for signals at the $3\sigma$ level.

Zoom In Zoom Out Reset image size

We also note that the orbital solution used by Martins et al. (2015) assumed a fixed, zero eccentricity, but did not explore the systemic velocity parameter space, which may have led to a stronger signal being detected at the offset ${V}_{\mathrm{sys}}$ we found in Section 3.2.1. If the optical data can be analyzed such that they do not induce the excess broadening seen in Martins et al. (2015), we can compare the detected planet ${V}_{\mathrm{sys}}$ at optical and infrared wavelengths. If the reflected light signal is not offset in ${V}_{\mathrm{sys}}$, while the infrared remains discrepant, then assuming it is astrophysical it potentially indicates an offset hot spot in the non-transiting planet's highly irradiated atmosphere, akin to that seen in the transiting hot Jupiters HD 189733 b and WASP-43 b (Knutson et al. 2007; Stevenson et al. 2014, respectively), enabling a means of mapping the atmospheres of close-in non-transiting planets. This would require small errors on the measured ${V}_{\mathrm{sys}}$, which can be significantly reduced by observing the planet over a range of orbital phase, especially before and after superior conjunction. To date, the best constrained ${V}_{\mathrm{sys}}$ of a hot Jupiter is for τ Boo b (Brogi et al. 2012) with an uncertainty of ±0.1 km s⁻¹ (an order of magnitude better than reported here for 51 Peg b). τ Boo is of similar brightness to 51 Peg b and required three half-nights of CRIRES observations to achieve this precision on ${V}_{\mathrm{sys}}$, a factor of three times longer than our 51 Peg b observations. Similar observing durations would be required in the optical. This suggests that mapping of non-transiting hot Jupiters could be achieved with current instrumentation on 4–8 m class telescopes.

New and future instruments, such as ESPRESSO/VLT, and G-CLEF/GMT, will enable more sensitive probes of the reflected light and the albedo of even smaller (non-)transiting planets, while newly commissioned, improved, and upcoming high-resolution infrared spectrographs (e.g., IGRINS, GIANO, iSHELL, CARMENES, NIRSPEC, and CRIRES+, to name a few; see Crossfield 2014 for an extensive list) will independently confirm the CRIRES infrared detections to date, and probe smaller, cooler planets.

Finally, we note that dedicated, stable, uninterrupted, and consistent observing sequences under stable atmospheric conditions best serve the analysis of exoplanet atmospheres at high spectral resolution.