TWO EXOPLANETS DISCOVERED AT KECK OBSERVATORY^*

We used the HIRES spectrometer (Vogt et al. 1994) at Keck Observatory for 3–5 years to obtain a temporal sequence of R = 65,000 spectra for each star. An iodine cell in the beam imprinted a rich set of molecular absorption lines on the stellar spectrum. These iodine lines constrain the wavelength scale, point-spread functions, and Doppler shift for each individual observation (Marcy & Butler 1992; Butler et al. 1996). An exposure meter was used to adjust each exposure time to achieve a consistent signal-to-noise ratio of 200 per extracted pixel, which alleviates some types of systematic errors. The exposure meter was also used to determine the photon-weighted midpoint of each exposure, which improves the precision of the barycentric velocity correction.

Spectra obtained after 2004 August with the new HIRES detector mosaic include the Ca ii H & K lines, which provide a diagnostic of chromospheric activity. We characterized line core emission in terms of the S index (Vaughan et al. 1978; Duncan et al. 1991). To improve measurement precision, we matched the line wings and neighboring continua for each observation of a given star to the mean value for all observations of that star (Isaacson 2009).

We obtained photometric observations of HD 179079 and HD 73534 with the T12 0.8 m automated photometric telescope (APT) at Fairborn Observatory in southern Arizona. The T12 APT and its precision photometer are very similar to the T8 APT described in Henry (1999). The precision photometer uses two temperature-stabilized EMI 9124QB photomultiplier tubes to measure photon count rates simultaneously through Strömgren b and y filters.

The telescope was programmed to measure each target star with respect to three nearby comparison stars in the following sequence: DARK, A, B, C, D, A, SKY_A, B, SKY_B, C, SKY_C, D, SKY_D, A, B, C, D, where A, B, and C are comparison stars and D is the program star. Each complete sequence, referred to as a group observation, was reduced to form three independent measures of each of the six differential magnitudes D−A, D−B, D−C, C−A, C−B, and B−A. The differential magnitudes were corrected for differential extinction with nightly extinction coefficients and transformed to the standard Strömgren system with yearly mean transformation coefficients. To filter out observations taken under nonphotometric conditions, an entire group of observation was discarded if the standard deviation of any of the six mean differential magnitudes exceeded 0.01 mag. We also combined the Strömgren b and y differential magnitudes into a single (b + y)/2 passband to improve the precision.

We determined stellar properties by analyzing one spectrum of each star obtained without the iodine cell. We used Spectroscopy Made Easy (SME) to model our spectra (Valenti & Piskunov 1996) with the procedure of Valenti & Fischer (2005) to fit each observed spectrum with a synthetic spectrum, obtaining the stellar effective temperature T_eff, surface gravity log g, metallicity [M/H], projected rotational velocity vsin i, and elemental abundances of Na, Si, Ti, Fe, and Ni. Our [M/H] parameter here (as in Valenti & Fischer 2005) scales solar abundances for elements other than Na, Si, Ti, Fe, and Ni that have significant spectral lines in our fitted wavelength intervals. We report iron abundance [Fe/H] rather than [M/H] because [Fe/H] is better defined. We use the abundance of Si relative to iron [Si/Fe] as a proxy for alpha-element enrichment [α/Fe].

We obtained a bolometric correction (BC) for each star by interpolating the "high-temperature" grid of VandenBerg & Clem (2003) to our spectroscopically determined T_eff, log g, and [Fe/H]. Using formulae and constants in Valenti & Fischer (2005), we calculated stellar luminosity from apparent visual magnitude V, BC, and distance d of the star, and then stellar radius R_⋆ from L_⋆ and T_eff.

Following Valenti & Fischer (2005), we determined stellar mass M_⋆, age, and surface gravity log g_iso by interpolating tabulated Yonsei–Yale isochrones (Demarque et al. 2004), first to our measured [α/Fe], then to our measured [Fe/H], then to our measured T_eff, and finally to the observed L_⋆. This yields zero, one, or multiple possible stellar models. We repeat the interpolation for a grid of 17⁴ values of [α/Fe], [Fe/H], T_eff, and L_⋆ that span the range of measurement uncertainties. We weight each outcome by the likelihood of the interpolation parameter and by the lifetime of the output evolutionary state. This favors relatively stable main-sequence states over rapidly evolving states.

The Valenti & Fischer (2005) procedure described above does not require agreement between log g from spectroscopy and log g_iso from isochrones, though the two are usually close. Here we introduce an outer iteration loop in which the SME analysis is repeated with log g fixed at the value of log g_iso from the preceding iteration. After a few iterations, the output log g from SME agrees with the input log g_iso. The price for this self-consistency is a greater dependence on models and worse χ²_ν for the spectrum fits. However, we are already completely dependent on stellar evolution models for M_⋆, and systematic errors in spectral line data dominate χ²_ν. Using evolutionary tracks as an additional constraint on log g while fitting spectra should improve the accuracy of the other derived parameters. We give examples later in this paper. Figure 1 shows the analysis procedure graphically.

Zoom In Zoom Out Reset image size

Using the periodogram procedure described in Marcy et al. (2005), we analyzed the velocity time series for each star to identify planet candidates and prospective periods. Once a sufficient number of observations were obtained, we fitted the velocities with Keplerian orbits, using a Levenberg–Marquardt algorithm (Wright & Howard 2009). The free parameters are orbital period P, velocity semiamplitude K, eccentricity e, time of periastron passage T_p, argument of periastron referenced to the line of nodes ω, and velocity offset of the center of mass γ.

When fitting Keplerian models, we adopted a total variance equal to the quadrature sum of our velocity measurement precision (1–2 m s⁻¹) and a velocity jitter term of 3.5 m s⁻¹. Jitter includes both systematic measurement errors and astrophysical velocity perturbations caused by photospheric flows and inhomogeneities. Our adopted jitter of 3.5 m s⁻¹ for the two subgiants in this paper yields χ²_ν near unity and is consistent with the range of values in Wright (2005). The exact choice of jitter has an insignificant effect on our derived values of P and K.

We used 10⁵ Monte Carlo trials to estimate uncertainties in our derived orbital parameters. In each trial, we constructed a simulated observation by adding errors to the Keplerian model that best matches our observed velocities. These errors were selected randomly from the distribution of observed residuals. A particular residual could be used multiple times or not at all in any given trial. We fitted each simulated observation using exactly the same procedure that we used to analyze the actual observation. Each trial included a periodogram analysis, which for stars with few observations can occasionally yield a distinct period. After 10⁵ trials, the width of the distribution function for a particular orbital parameter yields an estimate of the uncertainty in our observed value of that parameter.

HD 179079 (HIP 94256, V = 7.95, B − V = 0.744 ± 0.013) is a G5 subgiant at a distance of 65.5 ± 3.3 pc (ESA 1997; van Leeuwen 2008). High-resolution spectroscopic analysis (see Section 2.3) yields T_eff = 5684 ± 44 K, log g = 4.06 ± 0.06, [Fe/H] = +0.25 ± 0.03 dex, and v sin i < 1.0 km s⁻¹. These spectroscopic results imply a bolometric correction of BC = −0.094 and hence a stellar luminosity of L_⋆ = 2.41 ± 0.27 L_☉, where the uncertainty in luminosity is dominated by the uncertainty in distance. Using L_⋆ and T_eff, we obtain a stellar radius of R_⋆ = 1.60 ± 0.09 R_☉. Stellar evolutionary tracks imply a stellar mass of M_⋆ = 1.15 ± 0.03 M_☉.

HD 179079 is chromospherically inactive, based on weak emission in the cores of the Ca ii H & K lines. We measure S_HK= 0.153, which yields log R'_HK = −5.06 and implies a rotation period of roughly 38 d (Noyes et al. 1984). The level of chromospheric activity implies an approximate age of 7 Gyr, which is consistent with the 68% credible interval of 6.1–7.5 Gyr from the isochrone analysis. We estimate that surface motion in the photosphere of this early subgiant contributes 3.5 m s⁻¹ of velocity jitter that we add in quadrature to velocity measurement uncertainties, when modeling the data. Table 1 summarizes the stellar parameters. Table 2 shows how key stellar parameters changed during the iteration procedure introduced in Section 2.3.

We obtained 74 observations of HD 179079 beginning in 2004 July. The observation dates, radial velocities, and measurement uncertainties are listed in Table 3. Typical exposure times were about 2 minutes. The median velocity measurement uncertainty is 1.2 m s⁻¹, which is small compared to the velocity jitter of 3.5 m s⁻¹.

Figure 2 shows a periodogram of the 74 measured radial velocities, with an unambiguous peak in power at 14.46 days. The false alarm probability (FAP) is less than 0.0001, i.e., the probability that a random set of data would produce a peak with the observed power is less than 0.01%. The FAP test checks for spurious peaks that can arise because of window functions in the data. To calculate the FAP, simulated velocities for actual observation times were drawn randomly from the distribution of observed velocities, allowing reuse of any value. After 10⁴ trials, we found no peaks for simulated data with as much power as the observed peak at 14.48 days. Thus, the FAP is less than the reciprocal of the number of trials, i.e., less than 10⁻⁴.

Zoom In Zoom Out Reset image size

Our best-fitting Keplerian model has a period P = 14.476 ± 0.011 days, velocity semiamplitude K = 6.64 ± 0.60 m s⁻¹, and eccentricity e = 0.12 ± 0.09. Uncertainties for the orbital parameters were derived from Monte Carlo trials, as described in Section 2.4. The rms residual about the model fit is 3.9 m s⁻¹. Including jitter of 3.5 m s⁻¹, we find χ²_ν = 1.04. The eccentricity is poorly constrained, as indicated by the large uncertainty found in our Monte Carlo trials. A circular orbit yields χ²_ν = 1.03, which is as probable as the solution obtained with eccentricity and T_p as free parameters. The challenge of measuring low eccentricities in low amplitude systems is discussed in Section 5.1.

Figure 3 shows the phase-folded RV data together with the Keplerian (solid line) and circular (dotted line) models that best fit the data. Observations are plotted using phases based on the Keplerian fit, rather than the circular fit. In this and subsequent velocity plots, error bars are dominated by velocity jitter, rather than velocity measurement precision. Adopting a stellar mass of 1.146 ± 0.028 M_☉, we derive a minimum mass of M sin i = 27.5 ± 2.5 M_⊕ and a semimajor axis of 0.121 ± 0.001 AU. The complete Keplerian orbital solution is given in Table 4.

Zoom In Zoom Out Reset image size

Because the orbital period for this planet is relatively short, we have carried out an extensive search for transits, both with photometry (described in the following section) and by phase-folding the radial velocities to search for velocities that might have been obtained serendipitously during the ∼3 hr transit window. As the planet traverses the stellar disk, it blocks light from first one side (say the approaching, blueshifted edge of the stellar disk) and then the other side (the receding, redshifted edge) of the star. Our Doppler analysis interprets the spectral line asymmetry as excess velocity shifts during ingress and egress. This phenomenon is known as the Rossiter–McLaughlin effect.

The amplitude of the Rossiter–McLaughlin effect increases with both projected stellar rotation velocity and planet size relative to the star. HD 179079 has a low rotational velocity (v sin i = 0.5 km s⁻¹) and a small radius (expected to be similar to Neptune for an M sin i = 25.4 M_Jup planet). Therefore, the amplitude of the Rossiter–McLaughlin effect is expected to be no more than our measurement errors of about 2 m s⁻¹ for HD 179079. Nevertheless, we phase-folded the radial velocities and calculated the ingress and egress times for eccentricities ranging from 0 to 0.1. Uncertainties in the orbital eccentricity lead to shifts of up to 16 hr in the prospective transit time at the epochs of our radial velocity measurements. Velocity residuals during these broad transits windows were no larger than velocity residuals at other orbital phases.

Our 243 good brightness measurements of HD 179079 were made between 2007 June and 2008 June and cover parts of the 2007 and 2008 observing seasons. The comparison stars A, B, and C were HD 177552 (V = 6.51, B − V = 0.36, F1 V), HD 181420 (V = 6.57, B − V = 0.44, F2), and HD 180086 (V = 6.63, B − V = 0.35, F0), respectively. The differential magnitudes C−A, C−B, and B−A demonstrated that all three comparison stars were constant to 0.002 mag or better. To minimize the effect of any low-level intrinsic variation in the three comparison stars, we averaged the three D−A, D−B, and D−C differential magnitudes of HD 179079 within each group into a single value, representing the difference in brightness between HD 179079 and the mean brightness of the three comparison stars: D − (A + B + C)/3. The standard deviation of these ensemble differential magnitudes for the complete data set is 0.00190 mag. This is comparable to the typical precision of a single observation with this telescope, indicating there is little or no photometric variability in HD 179079.

Solar-type stars often exhibit brightness variations caused by cool, dark photospheric spots as they are carried into and out of view by stellar rotation (e.g., Gaidos et al. 2000). Periodogram analyses of the D−A, D−B, and D−C differential magnitudes for HD 179079 yield no significant periodicity between 1 and 100 days, consistent with the star's low level of chromospheric activity and its low vsin i (Table 1). We see no significant power at any period within a factor of 2 of 38 days, which is the rough rotation period implied by the chromospheric activity level.

The 243 ensemble (b + y)/2 differential magnitudes of HD 179079 are plotted in the top panel of Figure 4. Phases are computed from the orbital period given in Table 4 and the epoch JD 2, 454, 678.8 ± 0.67, a recent time of mid-transit derived from the orbital elements. A least-squares sine fit on the orbital period yields a semiamplitude of only 0.00007 ± 0.00016 mag. This very low limit to photometric variability on the radial velocity period is strong evidence that the low-amplitude radial velocity variations observed in the star are, in fact, due to reflex motion induced by a low-mass companion and not to activity-induced intrinsic variations in the star itself (e.g., Paulson et al. 2004).

Zoom In Zoom Out Reset image size

The photometric observations of HD 179079 near the predicted time of transit are replotted with an expanded horizontal scale in the bottom panel of Figure 4. The solid curve shows the predicted time (0.00 phase units) and duration (±0.009 phase units) of transits with a depth of 0.08% computed from estimated stellar and planetary radii. The error bar in the upper right of both panels represents the mean precision of a single observation (0.0019 mag). The horizontal error bar immediately below the transit in both panels represents the uncertainty in the predicted time of mid-transit (±0.67 days or ±0.046 phase units). It is clear from the data in the bottom panel that we cannot rule out the possibility of shallow transits of HD 179079b.

HD 73534 (HIP 42446, V = 8.23, B − V = 0.962 ± 0.021) is a G5 subgiant at a distance of 81.0 ± 4.9 pc (ESA 1997; van Leeuwen 2008). High-resolution spectrum synthesis modeling yields T_eff = 5041 ± 44 K, log g = 3.78 ± 0.06, [Fe/H] = +0.23 ± 0.03 dex, and v sin i < 1.0km s⁻¹. The resulting bolometric correction of BC = −0.266 yields a stellar luminosity of L_⋆ = 3.33 ± 0.43 L_☉. Combining T_eff and L_⋆, we obtain a stellar radius of R_⋆ = 2.39 ± 0.16 R_☉. Stellar evolutionary tracks yield a stellar mass of M_☉ = 1.23 ± 0.06 M_☉.

HD 73534 has minimal emission in the Ca ii H & K line cores, implying chromospheric inactivity, which is typical for subgiants (Wright 2004). Applying Noyes et al. (1984) relationships that were calibrated using main-sequence stars, our measured S_HK = 0.155 yields log R'_HK = −5.13 and a crude rotation period of 53 days. The level of chromospheric activity implies an approximate age of 9 Gyr, which differs significantly from the 68% credible interval of 5.2–7.2 Gyr from the isochrone analysis.

Subgiants have slightly more velocity jitter than inactive main-sequence stars, as evidenced by the greater rms scatter seen in subgiants without detected planets. We estimate that the intrinsic stellar jitter of HD 73534 is 3.5 m s⁻¹. Table 1 lists our derived stellar parameters. Table 5 shows how key parameters for HD 73435 changed during the iteration procedure introduced in Section 2.3.

We began observing HD 73534 in 2004 July as part of the N2K program (Fischer et al. 2005). No short-period velocity variations were detected, but we continued to obtain a few velocity measurements each year to map out an emerging low-amplitude, long-period planet. We now have a total of 30 observations that span five years. The raw velocity measurements have a precision of 1.1 m s⁻¹, but the (unmodeled) rms is 13.8 m s⁻¹. Exposure times were two to five minutes. The observation dates, radial velocities and associated uncertainties are listed in Table 6.

The best-fit Keplerian model has a period, P = 1770 ± 40 days, velocity semiamplitude, and K = 16.2 ± 1.1 m s⁻¹. The orbital eccentricity, e = 0.07 ± 0.07, is not significantly different from zero. The rms residual of the fit is 3.36 m s⁻¹ with χ²_ν= 0.91 after including stellar jitter of 3.5 m s⁻¹. Adopting a stellar mass of M_☉ = 1.228 ± 0.060 M_☉, we derive M sin i = 1.103 ± 0.087 M_Jup and a semimajor axis of a = 3.067 ± 0.068 AU. The complete set of orbital parameters are listed in Table 4. Figure 5 shows the phased radial velocity data with the best-fit Keplerian (solid line) and circular (dotted line) models overplotted.

Zoom In Zoom Out Reset image size

From 2004 November to 2008 December, we collected 521 good photometric observations of HD 73534 during five consecutive observing seasons. We do not see a correlation between radial velocities and activity or photometric measurements, however, we present the photometric data for posterity.

The comparison stars A, B, and C were HD 72943 (V = 6.33, B − V = 0.34, F0 IV), HD 73347 (V = 8.00, B − V = 0.40, F0), and HD 73821 (V = 7.82, B − V = 0.30, F0), respectively. Comparison star A (HD 72943) was found to be variable with an amplitude of 0.010 mag and a period of 0.0919 d; it is probably a δ Scuti variable. Comparison stars B and C were constant to 0.002 mag or better, so we created and analyzed ensemble differential magnitudes using only those two comparison stars: D − (B + C)/2.

The standard deviation of all 521 observations is 0.0018 mag, which is the typical precision of our photometry. To search for low-amplitude periodic brightness variations, we calculated power spectra for each observing season and combined the results. We did not detect any significant power for periods between 1 and 100 days, which more than spans the range of plausible rotation periods. The absence of significant optical variability on rotational time scales is consistent with the low level of chromospheric activity detected in HD 73534.

To search for low-amplitude, brightness variations on longer time scales, we computed the mean brightness for each of our five observing seasons. For HD 73534, the seasonal means have a full range of 0.0029 mag and a standard deviation of 0.0012 mag. For the two comparison stars, the seasonal mean values of C − B have a full range of only 0.0014 mag and a standard deviation of only 0.0006 mag, a factor of 2 smaller. Wright et al. (2008) demonstrated that for solar-type stars, we can measure seasonal means with a precision (standard deviation) of 0.0002 mag. HD 73534 varies by about 0.001 mag from year to year, but the variations are not systematic.

Finally, we compared our seasonal mean brightness measurements with our seasonal mean radial velocities, obtaining a linear correlation coefficient of 0.3345. With only five data points, the correlation coefficient would be at least this large 42% of the time, if brightness and velocity are uncorrelated random variables. A significant correlation between stellar brightness and radial velocity would raise doubts about the planetary origin of the stellar velocity variations, but in this case we do not detect a significant correlation.

Eccentricity cannot be negative. For a circular orbit, errors in individual radial velocity measurements can only drive the measured eccentricity away from the true value of zero (e.g., Shen & Turner 2008). As radial velocity planet searches push to lower amplitude systems, this bias becomes significant for low-eccentricity planets.

When observed velocity constraints for a circular orbit are uniformly distributed in orbital phase (approximately true for most radial velocity planet detections), Lucy (1971) derives the probability distribution for measured eccentricity,

$$\begin{equation} p(e)de={e\over \sigma _e^2}\exp \left(-{e^2\over 2\sigma _e^2}\right)de, \end{equation} \tag{ 1 }$$

where the eccentricity uncertainty, σ_e, is given by Luyten (1936) as

$$\begin{equation} \sigma _e={\sigma \over K}\left({2\over N}\right)^{0.5}, \end{equation} \tag{ 2 }$$

where σ is the typical uncertainty in velocity and N is the number of velocity measurements. The detection of nonzero eccentricity with better than 95% confidence requires approximately e/σ_e > 2.45.

To demonstrate the challenge in recovering zero eccentricity for HD 179079, we created 10⁵ synthetic data sets based on an e = 0 fit of the observed velocities for HD 179079. We added errors by drawing randomly (with replacement) from the observed distribution of residuals about the e = 0 fit. Finally, we fitted each simulated data set with a Keplerian, leaving eccentricity as a free parameter. Figure 6 shows the resulting distribution of imprecisely measured eccentricities.

Zoom In Zoom Out Reset image size

Even though we simulated a circular orbit, the measured eccentricities in Figure 6 are significantly biased toward positive values. The eccentricity distribution has a median of 0.164 and a standard deviation of 0.075. Clearly, our observed eccentricity of 0.115 ± 0.087 does not rule out a circular orbit.

Our N = 74 observations of HD 179079 are spread nearly uniformly in orbital phase and have an rms residual of σ = 3.9 m s⁻¹ (dominated by stellar jitter with an amplitude of about 3.5 m s⁻¹). With a measured eccentricity of e = 0.11 and a velocity semiamplitude of K = 6.6 m s⁻¹, we obtain e/σ_e = 1.1, which is well below the approximate threshold for a significant detection.

The dotted line in Figure 6 shows the predicted distribution according to Equation (1), but with σ = 4.2 m s⁻¹, rather than the observed rms of 3.9 m s⁻¹, to better match the observed distribution. This slight excess in σ_e relative to Equation (2) may be due to velocity constraints that are not uniformly distributed in phase or nonoptimal behavior of the Levenberg–Marquardt algorithm used to fit Keplerian orbits. The Keplerian fitting routine also returns the limiting value of e = 0 more often than predicted by the analytic approximation.

Figure 7 shows observed eccentricity distributions for 204 well-characterized planets, a subset of 163 planets with K > 20 m s⁻¹, and a further subset of 70 planets with K > 80 m s⁻¹. In all three cases, we exclude 53 planets with periods less than 10 days that may have experienced orbital evolution due to tidal interaction with the star.

Zoom In Zoom Out Reset image size

The observed eccentricity distributions become flatter as lower K planets are excluded. This flattening cannot be due to eccentricity measurement bias (Section 5.1), which would create the opposite slope, as positively biased measurements for lower K planets are included in the distribution. The sequence of observed eccentricity distributions in Figure 7 may reflect dynamical processes that link planet mass and eccentricity (e.g., Jurić & Tremaine 2008; Ford & Rasio 2008, for high-mass planets), but below we discuss briefly a possible origin based on observational bias.

At fixed K, planets in more eccentric orbits induce large velocity shifts for a smaller fraction of the orbit, making them harder to detect. For this reason, the dependence of the observed eccentricity distribution on K threshold may reflect an observational bias, rather than measurement bias. If true, the flat eccentricity distribution for K > 80 m s⁻¹ in Figure 7 may be the true distribution for giant planets. For the K > 20 m s⁻¹ threshold, detection of highly eccentric planets may become more difficult, causing the eccentricity distribution to drop above e = 0.4. Finally, the lowest K planets may be very difficult to detect reliably except in nearly circular orbits. If this hypothesis is true, a large population of low-K, high-e planets have yet to be detected. It would be interesting to quantify the expected detection bias as a function of eccentricity in actual radial velocity planet search programs.