RETIRED A STARS AND THEIR COMPANIONS. VII. 18 NEW JOVIAN PLANETS^*

The details of the target selection of our Doppler survey of evolved stars at Keck Observatory have been described in detail by, e.g., Johnson et al. (2006, 2010c) and Peek et al. (2009). In summary, we have selected subgiants from the Hipparcos catalog (van Leeuwen 2007) based on B − V colors and absolute magnitudes M_V so as to avoid K-type giants that are observed as part of other Doppler surveys (e.g., Hatzes et al. 2003; Sato et al. 2005; Reffert et al. 2006) and exhibit jitter levels in excess of 10 m s⁻¹ (Hekker et al. 2006). We also selected stars in a region of the temperature–luminosity plane in which stellar model grids of various masses are well separated and correspond to masses M_⋆ > 1.3 M_☉ at solar metallicity according to the Girardi et al. (2002) model grids. However, some of our stars have subsolar metallicities ([Fe/H] <0) and correspondingly lower masses down to ≈1 M_☉. Our sample of 240 subgiants monitored at Keck Observatory (excluding the Lick Observatory sample described by Johnson et al. 2006) is shown in Figure 1 and compared with the full target sample of the California Planet Survey (CPS).

Zoom In Zoom Out Reset image size

We obtained spectroscopic observations of our sample of subgiants at Keck Observatory using the High Resolution Echelle Spectrometer (HIRES) with a resolution of R ≈ 55,000 with the B5 decker (086 width) and red cross-disperser (Vogt et al. 1994). We use the HIRES exposure meter to ensure that all observations receive uniform flux levels independent of atmospheric transparency variations and to provide the photon-weighted exposure midpoint, which is used for the barycentric correction (BC). Under nominal atmospheric conditions, a V = 8 target requires an exposure time of 90 s and results in a signal-to-noise ratio of 190 at 5800 Å for our sample comprising mostly early K-type stars.

Normal program observations are made through a temperature-controlled Pyrex cell containing gaseous iodine, which is placed just in front of the entrance slit of the spectrometer. The dense set of narrow molecular lines imprinted on each stellar spectrum from 5000 Å to 6200 Å provides a robust, simultaneous wavelength calibration for each observation, as well as information about the shape of the spectrometer's instrumental response (Marcy & Butler 1992). RVs are measured with respect to an iodine-free "template" observation that has had the HIRES instrumental profile removed through deconvolution. Differential Doppler shifts are measured from each spectrum using the forward-modeling procedure described by Butler et al. (1996), with subsequent improvements over the years by the CPS team (e.g., Howard et al. 2011a). The instrumental uncertainty of each measurement is estimated based on the weighted standard deviation of the mean Doppler shift measured from each of ≈700 independent 2 Å spectral regions. In a few instances we made two or more successive observations of the same star and averaged the velocities in 2 hr time intervals, thereby reducing the associated measurement uncertainty.

We have also obtained additional spectra for HD 1502 in collaboration with the McDonald Observatory planet search team. A total of 54 RV measurements were collected for HD 1502: 32 with the 2.7 m Harlan J. Smith Telescope and its Tull Coude Spectrograph (Tull et al. 1995), and 22 with the High Resolution Spectrograph (HRS; Tull 1998) at the Hobby–Eberly Telescope (Ramsey et al. 1998). On each spectrometer we use a sealed and temperature-controlled iodine cell as velocity metric and to allow line-spread function reconstruction. The spectral resolving power for the HRS and Tull spectrograph is set to R = 60,000. Precise differential RVs are computed using the Austral I₂-data modeling algorithm (Endl et al. 2000).

The RV measurements are listed in Tables 1–18 together with the Heliocentric Julian Date (HJD) of observation and internal measurement uncertainties, excluding the jitter contribution described in Section 2.5.

We use the iodine-free template spectra to estimate atmospheric parameters of the target stars with the LTE spectroscopic analysis package Spectroscopy Made Easy (SME; Valenti & Piskunov 1996), as described by Valenti & Fischer (2005) and Fischer & Valenti (2005). Subgiants have lower surface gravities than dwarfs, and the damping wings of the Mg i b triplet lines therefore provide weaker constraints on the surface gravity, which is in turn degenerate with effective temperature and metallicity. To constrain log g, we use the iterative scheme of Valenti et al. (2009), which ties the SME-derived value of log g to the gravity inferred from interpolating the stellar luminosity, temperature, and metallicity onto the Yonsei–Yale (Y²; Yi et al. 2004) stellar model grids, which also give the stellar age and mass. The model-based log g is held fixed in a second SME analysis, and the process is iterated until convergence is met between the model-based and spectroscopically measured surface gravity, which results in best-fitting estimates of T_eff, log g, [Fe/H], and V_rsin i.

We perform our model-grid interpolations using a Bayesian framework similar to that described by Takeda et al. (2008). We incorporate prior constraints on the stellar mass based on the stellar initial mass function and the differential evolutionary timescales of stars in various regions of the theoretical H-R diagram. These priors tend to decrease the stellar mass inferred for a star of a given effective temperature, luminosity, and metallicity compared with a naive interpolation onto the stellar model grids (Lloyd 2011).

We determine the luminosity of each star from the apparent V-band magnitude and parallax from Hipparcos (van Leeuwen 2007) and the bolometric correction based on the effective temperature relationship given by VandenBerg & Clem (2003).¹⁰ Stellar radii are estimated using the Stefan–Boltzmann relationship and the measured L_⋆ and T_eff. We also measure the chromospheric emission in the Ca ii line cores (Wright et al. 2004; Isaacson & Fischer 2010), providing an S_HK value on the Mt. Wilson system.

The stellar properties of the 18 stars presented herein are summarized in Table 19.

We acquired photometric observations of 17 of the 18 planetary candidate host stars with the T3 0.4 m automatic photometric telescope (APT) at Fairborn Observatory. T3 observed each program star differentially with respect to two comparison stars in the following sequence, termed a group observation: K,S,C,V,C,V,C,V,C,S,K, where K is a check (or secondary comparison) star, C is the primary comparison star, V is the target star, and S is a sky reading. Three V − C and two K − C differential magnitudes are computed from each sequence and averaged to create group means. Group mean differential magnitudes with internal standard deviations greater than 0.01 mag were rejected to eliminate the observations taken under nonphotometric conditions. The surviving group means were corrected for differential extinction with nightly extinction coefficients, transformed to the Johnson system with yearly mean transformation coefficients, and treated as single observations thereafter. The precision of a single group-mean observation is usually in the range ∼0.003–0.006 mag (e.g., Henry et al. 2000), depending on the brightness of the stars within the group, the quality of the night, and the air mass of the observation. Further information on the operation of the T3 APT can be found in Henry et al. (1995a, 1995b) and Eaton et al. (2003).

Our photometric observations are useful for eliminating potential false positives from our sample of new planets. For example, Queloz et al. (2001) and Paulson et al. (2004) have demonstrated how rotational modulation in the visibility of starspots on active stars can result in periodic RV variations and, therefore, the potential for erroneous planetary detections. Photometric results for the 17 stars in the present sample are given in Table 20. Columns 7–10 give the standard deviations of the V − C and K − C differential magnitudes in the V and B passbands with 3σ outliers removed. All of the standard deviations are small and consistent with the measurement precision of the telescope. Periodogram analysis of each data set found no significant periodicity between 1 and 100 days.

We conclude that all 17 planetary candidate stars in Table 20, as well as all of their comparison and check stars, are constant to the limit of our photometric precision. The lack of evidence for photometric variability provides support for the planetary interpretation of the RV variations.

Although we do not have photometric measurements of HD 142245, we note from Table 19 that HD 142245 has one of the lowest values for S_HK in the sample. Therefore, like the rest of the sample, HD 142245 should be photometrically stable.

As in Johnson et al. (2010c), we perform a thorough search of the RV time series of each star for the best-fitting Keplerian orbital model using the partially linearized, least-squares fitting procedure described in Wright & Howard (2009) and implemented in the IDL package RVLIN.¹¹ The free parameters in our model are the velocity semiamplitude K, period P, argument of periastron ω, time of periastron passage T_p, and the systemic velocity offset γ. When fitting RVs from separate observatories, we include additional offsets γ_i for the different data sets. As described in Section 2.6, we also explore the existence of a constant acceleration $\dot{\gamma }$ in each RV time series.

In addition to the parameters describing the orbit, we also include an additional error contribution to our RV measurements due to stellar "jitter," which we denote by s. The jitter accounts for any unmodeled noise sources intrinsic to the star such as rotational modulation of surface inhomogeneities and pulsation (Saar et al. 1998; Wright 2005; Makarov et al. 2009; Lagrange et al. 2010) and is added in quadrature to the internal uncertainty of each RV measurement.

Properly estimating the jitter contribution to the uncertainty of each measurement is key to accurately estimating the confidence intervals for each fitted parameter. Ignoring jitter will lead to underestimated parameter uncertainties, rendering them less useful in future statistical investigations of exoplanet properties. Similarly, the equally common practice relying on a single value of the jitter based on stars with properties similar to the target of interest ignores variability in the jitter observed from star to star and can potentially overestimate the parameter uncertainties. For these reasons we take the approach of allowing the jitter term to vary in our orbit analyses, as described by, e.g., Ford & Gregory (2007).

We estimate parameter uncertainties using a Markov Chain Monte Carlo (MCMC) algorithm (see, e.g., Ford 2005; Winn et al. 2007). MCMC is a Bayesian inference technique that uses the data together with prior knowledge to explore the shape of the posterior probability density function (pdf) for each parameter of an input model. MCMC with the Metropolis–Hastings algorithm in particular provides an efficient means of exploring high-dimensional parameter space and mapping out the posterior pdf for each model parameter.

At each chain link in our MCMC analysis, one parameter is selected at random and is altered by drawing a random variate from a transition probability distribution. If the resulting value of the likelihood $\mathcal {L}$ for the trial orbit is greater than the previous value, then the trial orbital parameters are accepted and added to the chain. If not, then the probability of adopting the new value is set by the ratio of the probabilities from the previous and current trial steps. If the current trial is rejected, then the parameters from the previous step are adopted. The size of the transition function determines the efficiency of convergence. If it is too narrow, then the full exploration of parameter space is slow and the chain is susceptible to local minima; if it is too broad, then the chain exhibits large jumps and the acceptance rates are low.

Rather than minimizing χ²_ν, we maximize the logarithm of the likelihood of the data, given by

$$\begin{equation} \ln {\mathcal {L}} = -\sum _{i=1}^{N_{\rm obs}} \ln {\sqrt{2\pi (\sigma _i + s)^2}} - \frac{1}{2} \sum _{i=1}^{N_{\rm obs}} \left[\frac{v_i - v_{m}(t_i)}{\sigma _i + s}\right]^2, \end{equation} \tag{ 1 }$$

where v_i and σ_i are the ith velocity measurement and its associated measurement error, v_m(t_i) is the Keplerian model at time t_i, s is the jitter, and the sum is performed over all N_obs measurements. If s = 0, then the first term on the right side—the normalization of the probability—is a constant, and the second term becomes $\frac{1}{2}\chi ^2$. Thus, maximizing $\ln \mathcal {L}$ is equivalent to minimizing χ². Larger jitter values more easily accommodate large deviations of the observed RV from the model prediction, but only under the penalty of a decreasing (more negative) normalization term, which makes the overall likelihood smaller.

We impose uninformative priors for most of the free parameters (either uniform or modified Jeffreys; e.g., Gregory & Fischer 2010). The notable exception is jitter, for which we use a Gaussian prior with a mean of 5.1 m s⁻¹ and a standard deviation of 1.5 m s⁻¹ based on the distribution of jitter values for a similar sample of intermediate-mass subgiants from Johnson et al. (2010c).

We use the best-fitting parameter values from RVLIN as initial guesses for our MCMC analysis. We choose normal transition probability functions with constant (rather than adaptive) widths. The standard deviations are iteratively chosen from a series of smaller chains so that the acceptance rates for each parameter are between 20% and 30%; each main chain is then run for 10⁷ steps. The initial 10% of the chains are excluded from the final estimation of parameter uncertainties to ensure uniform convergence. We select the 15.9 and 84.1 percentile levels in the posterior distributions as the "1σ" confidence limits. In most cases the posterior probability distributions were approximately Gaussian.

To determine whether there is evidence for a linear velocity trend, we use two separate methods: the Bayesian Information Criterion (BIC; Schwarz 1978; Liddle 2004) and inspection of the MCMC posterior pdfs, as described by Bowler et al. (2010). The BIC rewards better-fitting models but penalizes overly complex models, and is given by

$$\begin{equation} \mathrm{BIC} \equiv -2 \ln \mathcal {L_\mathrm{max}} + k \ln N, \end{equation} \tag{ 2 }$$

where $\mathcal {L_\mathrm{max}}$ is the maximum likelihood for a particular model with k free parameters and N data points. The relationship between $\mathcal {L_\mathrm{max}}$ and χ²_min is only valid under the assumption that the RVs are normally distributed, which is approximately valid for our analyses. A difference of ≳2 between BIC values with and without a trend indicates that there is sufficient evidence for a more complex model (Kuha 2004).

We also use the MCMC-derived pdf for the velocity trend parameter to estimate the probability that a trend is actually present in the data. We only adopt the model with the trend if the 99.7 percentile of the pdf lies above or below 0 m s⁻¹ yr⁻¹. The BIC and MCMC methods yield consistent results for the planet candidates presented in Section 3, and in many cases the RV trend is evident by visual inspection of Figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size