THE DISCOVERY OF HD 37605c AND A DISPOSITIVE NULL DETECTION OF TRANSITS OF HD 37605b^*

Jupiter analogs orbiting other stars represent the first signposts of true solar system analogs, and the eccentricity distribution of these planets with a > 3 AU will reveal how rare or frequent true Jupiter analogs are. To date, only nine "Jupiter analogs" have been well characterized in the peer reviewed literature¹² (defined here as P > 8 years, 4 > Msin i > 0.5 M_Jup, and e < 0.3; Wright et al. 2011, exoplanets.org). As the duration of existing planet searches approaches 10–20 years, more Jupiter analogs will emerge from their longest-observed targets (Wittenmyer et al. 2012; Boisse et al. 2012).

Of the over 700 exoplanets discovered to date, nearly 200 are known to transit their host star (Wright et al. 2011, exoplanets.org; Schneider et al. 2011, exoplanet.eu), and thousands more candidates have been discovered by the Kepler telescope. Of all of these planets, only three orbit stars with V < 8 and all have P < 4 days¹³. Long-period planets are less likely than close-in planets to transit unless their orbits are highly eccentric and favorably oriented, and indeed only two transiting planets with P > 20 days have been discovered around stars with V < 10, and both have e > 0.65 (HD 80606, Laughlin et al. 2009; Fossey et al. 2009; HD 17156, Fischer et al. 2007; Barbieri et al. 2007; both highly eccentric systems were discovered first with radial velocities).

Long-period planets not known to transit can have long transit windows due to both the large duration of any edge-on transit and higher phase uncertainties (since such uncertainties scale with the period of the orbit). Long-term radial velocity (RV) monitoring of stars, for instance, for the discovery of low-amplitude signals, can produce collateral benefits in the form of orbit refinement for a transit search and the identification of Jupiter analogs (e.g., Wright et al. 2009). Herein, we describe an example of both.

The inner planet in the system, HD 37605b, was the first planet discovered with the Hobby–Eberly Telescope (HET) at McDonald Observatory (Cochran et al. 2004). It is a super-Jupiter (Msin i = 2.41 M_Jup) on an eccentric orbit e = 0.67 with an orbital period in the "period valley" (P = 55 days; Wright et al. 2009).

W.C., M.E., and P.J.M. of the University of Texas at Austin, continued observations in order to get a much better orbit determination and to begin searching for transits. With the first new data in the fall of 2004, it became obvious that another perturber was present in the system, first from a trend in the RV residuals (i.e., a non-zero dv/dt; Wittenmyer et al. 2007), and later from curvature in the residuals. By 2009, the residuals to a one-planet fit were giving reasonable constraints on the orbit of a second planet, HD 37605c, and by early 2011 the orbital parameters of the c component were clear, and the Texas team was preparing the system for publication.

The Transit Ephemeris Refinement and Monitoring Survey (TERMS; Kane et al. 2009) seeks to refine the ephemerides of the known exoplanets orbiting bright, nearby stars with sufficient precision to efficiently search for the planetary transits of planets with periastron distances greater than a few hundredths of an AU (Kane et al. 2011b; Pilyavsky et al. 2011; Dragomir et al. 2011). This will provide the radii of planets not experiencing continuous high levels of insolation around nearby, easily studied stars.

In 2010, S.M. and J.T.W. began RV observations of HD 37605b at HET from Penn State University for TERMS, to refine the orbit of that planet for a future transit search. These observations, combined with Keck radial velocities from the California Planet Survey (CPS) consortium from 2006 onward, revealed that there was substantial curvature to the RV residuals to the original Cochran et al. (2004) solution. In 2010 October monitoring was intensified at HET and at Keck Observatory by A.W.H., G.W.M., J.T.W., and H.I., and with these new RV data and the previously published measurements from Wittenmyer et al. (2007) they obtained a preliminary solution for the outer planet. The discrepancy between the original orbital fit and the new fit (assuming one planet) was presented at the 2011 January meeting of the American Astronomical Society (Kane et al. 2011a).

In early 2011, the Texas and TERMS teams combined efforts and began joint RV analysis, dynamical modeling, spectroscopic analysis, and photometric observations (Kane et al. 2012). The resulting complete two-planet orbital solution allows for a sufficiently precise transit ephemeris for the b component to be calculated for a thorough transit search. We herein report the transit exclusion of HD 37605b and a stable dynamical solution to the system.

In Section 2, we describe our spectroscopic observations and analysis, which provided the radial velocities and the stellar properties of HD 37605. Section 3 details the orbital solution for the HD 37605 system, including a comparison with MCMC Keplerian fits, and our dynamical analysis. We report our photometric observations on HD 37605 and the dispositive null detection¹⁴ of non-grazing transits of HD 37605b in Section 4. After Section 5, Summary and Conclusion, we present updates on Msin i of two previously published systems (HD 114762 and HD 168443) in Section 6. In the Appendix we describe the algorithm used in the package BOOTTRAN (for calculating orbital parameter error bars, see Section 3.2).

Observations on HD 37605 at HET started 2003 December. In total, 101 RV observations took place over the course of almost 8 years, taking advantage of the queue scheduling capabilities of HET. The queue scheduling of HET allows for small amounts of telescope time to be optimally used throughout the year, and for new observing priorities to be implemented immediately, rather than on the next allocated night or after the time allocation and scheduling process (Shetrone et al. 2007). The observations were taken through the High Resolution Spectrograph (HRS; Tull 1998) situated at the basement of the HET building. This fiber-fed spectrograph has a typical long-term Doppler error of 3–5 m s⁻¹ (Baluev 2009). The observations were taken with the spectrograph configured at a resolving power of R = 60,000. For more details, see Cochran et al. (2004).

Observations at Keck were taken starting 2006 August. A set of 33 observations spanning over 5 years were made through the High Resolution Echelle Spectrometer (HIRES) (Vogt et al. 1994) on the Keck I telescope, which has a long-term Doppler error of 0.9–1.5 m s⁻¹ (e.g., Howard et al. 2009). The observations were taken at a resolving power of R = 55,000. For more details, see Howard et al. (2009) and Valenti et al. (2009).

Both our HET and Keck spectroscopic observations were taken with an iodine cell placed in the light path to provide wavelength standard and information on the instrument response function¹⁵ (IRF) for RV extraction (Marcy & Butler 1992; Butler et al. 1996). In addition, we also have observations taken without an iodine cell to produce stellar spectrum templates—on HET and Keck, respectively. The stellar spectrum templates, after being deconvolved with the IRF, are necessary for both RV extraction and stellar property analysis. The typical working wavelength range for this technique is roughly 5000–6000 Å.

In this section, we describe our data reduction and Doppler analysis of the HET observations. We reduced the Keck data with the standard CPS pipeline, as described in, for example, Howard et al. (2011) and Johnson et al. (2011a).

We have constructed a complete pipeline for analyzing HET data—from raw data reduction to RV extraction. The raw reduction is done using the REDUCE package by Piskunov & Valenti (2002). This package is designed to optimally extract echelle spectra from two-dimensional (2D) images (Horne 1986). Our pipeline corrects for cosmic rays and scattered light. In order to make the data reduction process completely automatic, we have developed our own algorithm for tracing the echelle orders of HRS and replaced the original semi-automatic algorithm from the REDUCE package.

After the raw data reduction, the stellar spectrum template is deconvolved using IRF derived from an iodine flat on the night of observation. There were two deconvolved stellar spectrum templates (DSST) derived from HET/HRS observations and one from Keck/HIRES. Throughout this work, we use the Keck DSST, which is of better quality thanks to a better known IRF of HIRES and a superior deconvolution algorithm in the CPS pipeline (Howard et al. 2009, 2011).

Then the pipeline proceeds with barycentric correction and RV extraction for each observation. We have adopted the Doppler code from CPS (e.g., Howard et al. 2009, 2011; Johnson et al. 2011a). The code is tailored to be fully functional with HET/HRS-formatted spectra, and it is capable of working with either an HET DSST or a Keck one.

The 101 HET RV observations include 44 observations which produced the published velocities in Cochran et al. (2004) and Wittenmyer et al. (2007), 34 observations also done by the Texas team in follow-up work after 2007, and 23 observations taken as part of the TERMS program. We have performed re-reduction on these 44 observations together with the rest of the 57 HET observations through our pipeline. This has the advantage of eliminating one free parameter in the Keplerian fit—the offset between two Doppler pipelines.

Two out of the 101 HET observations were excluded due to very low average signal-to-noise ratio per pixel (<20), and one observation taken at twilight was also rejected as such observation normally results in low accuracy due to the significant contamination by the residual solar spectrum (indeed this velocity has a residual of over 100 m s⁻¹ against best Keplerian fit, much larger than the ∼8  m s⁻¹ RV error).

All the HET and Keck radial velocities used in this work (98 from HET and 33 from Keck) are listed in Table 1.

HD 37605 is a K0 V star (V ∼ 8.7) with high proper motion at a distance of 44.0 ± 2.1 pc (ESA 1997; van Leeuwen 2008). We derived its stellar properties based on analysis on a high-resolution spectrum taken with Keck HIRES (without iodine cell in the light path). Table 2 lists the results of our analysis,¹⁶ including the effective temperature T_eff, surface gravity log g, iron abundance [Fe/H], projected rotational velocity vsin i, bolometric correction BC, bolometric magnitude M_bol, stellar luminosity L_⋆, stellar radius R_⋆, stellar mass M_⋆, and age. HD 37605 is found to be a metal-rich star ([Fe/H] ∼ 0.34) with M_⋆ ∼ 1.0 M_☉ and R_⋆ ∼ 0.9 R_☉.

We followed the procedure described in Valenti & Fischer (2005) and also in Valenti et al. (2009) with improvements. Briefly, the observed spectrum is fitted with a synthetic spectrum using Spectroscopy Made Easy (SME; Valenti & Piskunov 1996) to derive T_eff, log g, [Fe/H], vsin i, and so on, which are used to derive the bolometric correction BC and L_⋆ consequently. Then an isochrone fit by interpolating tabulated Yonsei-Yale isochrones (Demarque et al. 2004) using derived stellar parameters from SME is performed to calculate M_⋆ and log g_iso values (along with age and stellar radius). Next, Valenti et al. (2009) introduced an outside loop which re-runs SME with log g fixed at log g_iso, followed by another isochrone fit deriving a new log log g_iso using the updated SME results. The loop continues until log g values converge. This additional iterative procedure to enforce self-consistency on log g is shown to improve the accuracy of other derived stellar parameters (Valenti et al. 2009). The stellar radius and log g reported here in Table 2 are derived from the final isochrone fit, which are consistent with the purely spectroscopic results. The gravity (log g = 4.51) is also consistent with the purely spectroscopic gravity (4.44) based on strong Mg b damping wings, so for HD 37605 the iteration process is optional.

Cochran et al. (2004) reported the values of T_eff, log g, and [Fe/H] for HD 37605, and their estimates agree with ours within 1σ uncertainty. Santos et al. (2005) also estimated T_eff, log g, [Fe/H], and M_⋆, all of which agree with our values within 1σ. Our stellar mass and radius estimates are also consistent with the ones derived from the empirical method by Torres et al. (2010).

Our SME analysis indicates that the rotation of the star (vsin i) is likely <1 km s⁻¹ (corresponding to rotation period ≳ 46 days). We have used various methods to estimate stellar parameters from the spectrum, including the incorporation of color and absolute magnitude information and the Mg b triplet to constrain log g, and various macroturbulent velocity prescriptions. All of these approaches yield results consistent with an undetectable level of rotational broadening, with an upper limit of 1–2 km s⁻¹, consistent with the tentative photometric period 57.67 days derived from the Automatic Photoelectric Telescope (APT) data (see Section 4.1).

The traditional parameters for reporting the ephemerides of spectroscopic binaries are P, K, e, ω, and T_p, the last being the time of periastron passage (Wright & Howard 2009). This information is sufficient to predict the phase of a planet at any point in the future in principle, but the uncertainties in those parameters alone are insufficient to compute the uncertainty in orbital phase without detailed knowledge of the covariances among the parameters.

This problem is particularly acute when determining transit or secondary eclipse times for planets with near-circular orbits, where $\sigma _{T_p}$ and σ_ω can be highly covariant. In such cases the circular case is often not excluded by the data, and so the estimation of e includes the case e = 0, where ω is undefined. If the best or most likely value of e in this case is small but not zero, then it is associated with some nominal value of ω, but σ_ω will be very large (approaching π). Since T_p represents the epoch at which the true anomaly equals 0, T_p will have a similarly large uncertainty (approaching P), despite the fact that the phase of the system may actually be quite precisely known!

In practice even the ephemerides of planets with well-measured eccentricities suffer from a lack of knowledge of the covariance in parameters, in particular T_p and P (whose covariance is sensitive to the approximate epoch chosen for T_p). To make matters worse, the nature of "1σ" uncertainties in the literature is inconsistent. Some authors may report uncertainties generated while holding all or some other parameters constant (for instance, by seeing at what excursion from the nominal value χ² is reduced by 1), while others using bootstrapping or MCMC techniques may report the variance in a parameter over the full distribution of trials. In any case, covariances are rarely reported, and in some cases authors even report the most likely values on a parameter-by-parameter basis rather than a representative "best fit," resulting in a set of parameters that is not self-consistent.

The TERMS strategy for refining ephemerides therefore begins with the recalculation of transit time uncertainties directly from the archival RV data. We used bootstrapping (see the Appendix) with the time of conjunction, T_c (equivalent to transit center, in the case of transiting planets) computed independently for each trial. For systems whose transit time uncertainty makes definitive observations implausible or impossible due to the accumulation of errors in phase with time, we sought additional RV measurements to "lock down" the phase of the planet.

There are in total 137 radial velocities used in the Keplerian fit for the HD 37605 system. In addition to the 98 HET velocities and 33 Keck ones (see Section 2.2), we also included six velocities from Cochran et al. (2004)¹⁷ which were derived from observations taken with the McDonald Observatory 2.1 m Telescope (hereafter the 2.1 m telescope).

We used the RVLIN package by Wright & Howard (2009) to perform the Keplerian fit. This package is based on the Levenberg–Marquardt algorithm and is made efficient in searching parameter space by exploiting the linear parameters. The uncertainties of the parameters are calculated through bootstrapping (with 1000 bootstrap replicates) using the BOOTTRAN package, which is described in detail in the Appendix.¹⁸

The best-fit Keplerian parameters are listed in Table 3. The joint Keplerian fit for HD 37605b and HD 37605c has 13 free parameters: the orbital period P, time of periastron passage T_p, velocity semi-amplitude K, eccentricity e, and the argument of periastron referenced to the line of nodes ω for each planet; and for the system, the velocity offset between the center of the mass and barycenter of solar system γ and two velocity offsets between the three telescopes (Δ_Keck and Δ_HET, with respect to the velocities from the 2.1 m telescope as published in Cochran et al. 2004). We did not include any stellar jitter or RV trend in the fit (i.e., fixed to zero). The RV signals and the best Keplerian fits for the system, HD 37605b only, and HD 37605c only are plotted in the three panels of Figure 1, respectively.

Zoom In Zoom Out Reset image size

Adopting a stellar mass of M_⋆ = 1.000 ± 0.017 M_☉ (as in Table 2), we estimated the minimum mass (Msin i) for HD 37605b to be 2.802 ± 0.011 M_Jup and 3.366 ± 0.072 M_Jup for HD 37605c. While HD 37605b is on a close-in orbit at a = 0.2831 ± 0.0016 AU that is highly eccentric (e = 0.6767 ± 0.0019), HD 37605c is found to be on a nearly circular orbit (e = 0.013 ± 0.015) out at a = 3.814 ± 0.058 AU, which qualifies it as one of the "Jupiter analogs."

In order to discover whether the period and mass of the outer planet, HD 37605c, are well constrained, we mapped out the χ²_ν values for the best Keplerian fit in the P_c−M_csin i space (subscript "c" denoting parameters for the outer planet, HD 37605c). Each χ²_ν value on the P_c−M_csin i grid was obtained by searching for the best-fit model while fixing the period P_c for the outer planet and requiring constraints on K_c and e_c to maintain Msin i fixed. As shown in Figure 2, our data are sufficient to have both P_c and M_csin i well constrained. This is also consistent with the tight sampling distributions for P_c and M_csin i found in our bootstrapping results.

Zoom In Zoom Out Reset image size

The rms values against the best Keplerian fit are 7.86 m s⁻¹ for HET, 2.08 m s⁻¹ for Keck, and 12.85 m s⁻¹ for the 2.1 m telescope. In the case of HET and Keck, their rms values are slightly larger than their typical reported RV errors (∼5 m s⁻¹ and ∼1 m s⁻¹, respectively). This might be due to stellar jitter or underestimated systematic errors in the velocities. We note that the χ²_ν is reduced to 1.0 if we introduce a stellar jitter of 3.6 m s⁻¹ (added in quadrature to all the RV errors).

We compared our best Keplerian fit from RVLIN and uncertainties derived from BOOTTRAN (abbreviated as RVLIN+BOOTTRAN hereafter) with that from a Bayesian framework following Ford (2005) and Ford (2006) (referred to as the MCMC analysis hereafter). Table 4 lists the major orbital parameters from both methods for a direct comparison. Figure 3 illustrates this comparison, but with the MCMC results presented in terms of 2D confidence contours for P, e, K, Msin i, and ω of both planets, as well as for T_c of HD 37605b.

Zoom In Zoom Out Reset image size

For the Bayesian analysis, we assumed priors that are uniform in log of orbital period, eccentricity, argument of pericenter, mean anomaly at epoch, and the velocity zero-point. For the velocity amplitude (K) and jitter (σ_j), we adopted a prior of the form p(x) = (x + x_o)⁻¹[log(1 + x/x_o)]⁻¹, with K_o = σ_{j, o} = 1 m s⁻¹, i.e., high values are penalized. For a detailed discussion of priors, strategies to deal with correlated parameters, the choice of the proposal transition probability distribution function, and other details of the algorithm, we refer the reader to the original papers: Ford (2005, 2006) and Ford & Gregory (2007). The likelihood for RV terms assumes that each RV observation (v_i) is independent and normally distributed about the true RV with a variance of σ²_i + σ_j², where σ_i is the published measurement uncertainty. σ_j is a jitter parameter that accounts for additional scatter due to stellar variability, instrumental errors, and/or inaccuracies in the model (i.e., neglecting planet–planet interactions or additional, low-amplitude planet signals).

We used an MCMC method based upon Keplerian orbits to calculate a sample from the posterior distribution (Ford 2006). We calculated five Markov chains, each with ∼2 × 10⁸ states. We discarded the first half of the chains and calculated Gelman–Rubin test statistics for each model parameter and several ancillary variables. We found no indications of non-convergence among the individual chains. We randomly drew 3 × 10⁴ solutions from the second half of the Markov chains, creating a sample set of the converged overall posterior distribution of solutions. We then interrogated this sample on a parameter-by-parameter basis to find the median and 68.27% (1σ) values reported in Table 4. We refer to this solution set below as the "best-fit" MCMC solutions.

We note that the periods of the two planets found in this system are very widely separated (P_c/P_b ∼ 50), so we do not expect planet–planet interactions to be strong, hence we have chosen to forgo a numerically intensive N-body DEMCMC fitting procedure (see e.g., Johnson et al. 2011b; Payne & Ford 2011) as the non-Keplerian perturbations should be tiny (detail on the magnitude of the perturbations is provided in Section 3.4). However, to ensure that the Keplerian fits generated are stable, we took the results of the Keplerian MCMC fits and injected those systems into the Mercury n-body package (Chambers 1999) and integrated them forward for ∼10⁸ years. This allows us to verify that all of the selected best-fit systems from the Keplerian MCMC analysis are indeed long-term stable. Further details on the dynamical analysis of the system can be found in Section 3.4.

We assumed that all systems are coplanar and edge-on for the sake of this analysis, hence all of the masses used in our n-body analyses are minimum masses.

As shown in Table 4 and Figure 3, the parameter estimates from RVLIN+BOOTTRAN and MCMC methods agree with each other very well (all within 1σ error bar). In some cases, the MCMC analysis reports error bars slightly larger than bootstrapping method (∼20% at most). We note that the relatively large MCMC confidence intervals are not significantly reduced if one conducts an analysis at a fixed jitter level (e.g., σ_J = 3.5 m s⁻¹) unless one goes to an extremely low jitter value (e.g., ∼1.5 m s⁻¹). That is, the larger MCMC error bars do not simply result from treating the jitter as a free parameter. For the uncertainties on minimum planet mass Msin i and semi-major axes a, the MCMC analysis does not incorporate the errors on the stellar mass estimate. Note here, as previously mentioned in Section 3.1, that the "best-fit" parameters reported by the MCMC analysis here listed in Table 4 are not a consistent set, as the best estimates were evaluated on a parameter-by-parameter basis, taking the median from marginalized posterior distribution of each. Assuming no jitter, the best Keplerian fit from RVLIN has a reduced chi-square value χ²_ν = 2.28, while the MCMC parameters listed in Table 4 give a higher χ²_ν value of 2.91.

We used the best-fit Keplerian MCMC parameters as the basis for a set of long-term numerical (n-body) integrations of the HD 37605 system using the Mercury integration package (Chambers 1999). We used these integrations to verify that the best-fit systems (1) are long-term stable, (2) do not exhibit significant variations in their orbital elements on the timescale of the observations (justifying the assumption that the planet–planet interactions are negligible), (3) do not exhibit any other unusual features. We emphasize again that the planets in this system are well separated and we do not expect any instability to occur: for the masses and eccentricities in question, a planet at a_b ∼ 0.28 AU will have companion orbits which are Hill stable for a ≳ 0.83 AU (Gladman 1993), so while Hill stability does not preclude outward scatter of the outer planet, the fact that a_c ∼ 3.8 ≫ 0.83 AU suggests that the system will be far from any such instability.

We integrated the systems for >10⁸ years (∼10⁷ × the orbital period of the outer planet and >10² × the secular period of the system), and plot in Figure 4 the evolution of the orbital elements a, e, and, ω. On the left-hand side of the plot we provide short-term detail, illustrating that over the ∼10 year time period of our observations, the change in orbital elements will be very small. On the right-hand side we provide a much longer-term view, plotting 10⁷ out of >10⁸ years of system evolution, demonstrating that (1) the secular variation in some of the elements (particularly the eccentricity of the outer planet; see e_c in red) over a time span of ∼4 × 10⁵ years can be significant: in this case we see 0.03 < e_c < 0.11; but (2) the system appears completely stable, as one would expect for planets with a period ratio P_c/P_b ∼ 50. Finally, at the bottom of the figure we display the range of parameter space covered by the e_icos ω_i, e_isin ω_i parameters (i = b in blue for inner planet and i = c in red for outer planet), demonstrating that the orbital alignments circulate, i.e., they do not show any signs of resonant confinement, which confirms our expectation of minimal planet–planet interaction as mentioned before.

Zoom In Zoom Out Reset image size

As noted above, our analysis assumed coplanar planets. As such the planetary masses used in these dynamical simulations are minimum masses. We note that for inclined systems, the larger planetary masses will cause increased planet–planet perturbations. To demonstrate this is still likely to be unimportant, we performed a 10⁸ year simulation of a system in which 1/sin i = 10, pushing the planetary masses to ∼30 M_Jup. Even in such a pathological system the eccentricity oscillations are only increased by a factor of ∼2 and the system remains completely stable for the duration of the simulation.

We also performed a separate transit timing variation (TTV) analysis, using the best-fit MCMC systems as the basis for a set of highly detailed short-term integrations. From these we extracted the times of transit and found a TTV signal ∼100 s, or ∼0.001 day, which is much smaller than the error bar on T_c (∼0.07 days). Therefore we did not take into account the effect of TTV when performing our transit analysis in the next section.

The coincidence of the solar activity cycle period of 11 years and Jupiter's orbital period of nearly 12 years illustrates how activity cycles could, if they induced apparent line shifts in disk-integrated stellar spectra, confound attempts to detect Jupiter analogs around Sun-like stars. Indeed, Dravins (1985) predicted apparent RV variations of up to 30 m s⁻¹ in solar lines due to the solar cycle, and Deming et al. (1987) reported a tentative detection of such a signal in NIR CO lines of 30 m s⁻¹ in just 2 years, and noted that such an effect would severely hamper searches for Jupiter analogs. That concern was further amplified when Campbell et al. (1991) reported a positive correlation between RV and chromospheric activity in the active star κ¹ Cet, with variations of order 50–100 m s⁻¹.

Wright et al. (2008) found that the star HD 154345 has an apparent Jupiter analog (HD 154345 b), but that this star also shows activity variations in phase with the RV variations. They noted that many Sun-like stars, including the precise RV standard star HD 185144 (σ Dra) show similar activity variations and that rarely, if ever, are these signals well correlated with signals similar in strength to that seen in HD 154345 (∼15 m s⁻¹), and concluded that the similarity was therefore likely just an inevitable coincidence. Put succinctly, activity cycles in Sun-like stars are common (Baliunas et al. 1995), but few Jupiter analogs have been discovered, meaning that the early concern that activity cycles would mimic giant planets is not a severe problem.

Nonetheless, there is growing evidence that activity cycles can, in some stars, induce RV variations (Santos et al. 2010; Dumusque et al. 2012), and the example of HD 154345 still warrants care and concern. Most significantly, Dumusque et al. (2011) found a positive correlation between chromospheric activity and precise RV in the average measurements of a sample of HARPS stars, and Lovis et al. (2011) provided a formula for predicting the correlation strength as a function of the metallicity and effective temperature of the star. Their formulae predict a value of 2 m s⁻¹ for the most suspicious case in the literature, HD 154345 (compared to an actual semi-amplitude of ∼15 m s⁻¹), but are rather uncertain. It is possible that in a few, rare cases, the formula might significantly underestimate the amplitude of the effect.

The top panel of Figure 5 plots the T12 APT observations from all five observing seasons (data provided in Table 5; see the details on APT photometry in Section 4.1). The dashed line marks the mean relative magnitude (Δ(b + y)/2) of the first season. The seasonal mean brightness of the star increases gradually from year to year by a total of ∼0.002 mag, which may be due to a weak long-term magnetic cycle. However, no evidence is found in support of such a cycle in the Mount Wilson chromospheric Ca ii H and K indices (Isaacson & Fischer 2010), although the S values vary by approximately 0.1 over the span of a few years. The formulae of Lovis et al. (2011) predict a corresponding RV variation of less than 2 m s⁻¹ due to activity, far too small to confound our planet detection with K = 49 m s⁻¹.

Zoom In Zoom Out Reset image size

Since we do not have activity measurements for this target over the span of the outer planet's orbit in HD 37605, we cannot definitively rule out activity cycles as the origin of the effect, but the strength of the outer planetary signal and the lack of such signals in other stars known to cycle strongly dispels concerns that the longer signal is not planetary in origin.

The T12 0.8 m APT, located at Fairborn Observatory in southern Arizona, acquired 696 photometric observations of HD 37605 between 2008 January 16 and 2012 April 7. Henry (1999) provides detailed descriptions of observing and data reduction procedures with the APTs at Fairborn. The measurements reported here are differential magnitudes in Δ(b + y)/2, the mean of the differential magnitudes acquired simultaneously in the Strömgren b and y bands with two separate EMI 9124QB bi-alkali photomultiplier tubes. The differential magnitudes are computed from the mean of three comparison stars: HD 39374 (V = 6.90, B − V = 0.996, K0 III), HD 38145 (V = 7.89, B − V = 0.326, F0 V), and HD 38779 (V = 7.08, B − V = 0.413, F4 IV). This improves the precision of each individual measurement and helps to compensate for any real microvariability in the comp stars. Intercomparison of the differential magnitudes of these three comp stars demonstrates that all three are constant to 0.002 mag or better from night to night, consistent with typical single-measurement precision of the APT (0.0015–0.002 mag; Henry 1999).

Figure 5 illustrates the APT photometric data and our transit search. As mentioned in Section 3.5, the top panel shows all of our APT photometry covering five observing seasons, which exhibits a small increasing trend in the stellar brightness. To search for the transit signal of HD 37605b, the photometric data were normalized so that all five seasons had the same mean (referred to as the "normalized photometry" hereafter). The data were then phased at the orbital period of HD 37605b, 55.01307 days, and the predicted time of mid-transit, T_c, defined as Phase 0. The normalized and phased data are plotted in the middle panel of Figure 5. The solid line is the predicted transit light curve, with the predicted transit duration (0.352 days or 0.0064 phase unit) and transit depth (1.877% or ∼0.020 mag) as estimated above. The scatter of the phased data from their mean is 0.00197 mag, consistent with APT's single-measurement precision, and thus demonstrates that the combination of our photometric precision and the stability of HD 37605 is easily sufficient to detect the transits of HD 37605b in our phased data set covering 5 years. A least-squares sine fit of the phased data gives a very small semi-amplitude of 0.00031 ± 0.00011 mag (consistent with zero) and so provides strong evidence that the observed radial velocity variations are not produced by rotational modulation of surface activity on the star.

The bottom panel of Figure 5 plots the phased data around the predicted time of mid-transit, T_c, at an expanded scale on the abscissa. The horizontal bar below the transit window represents the ±1σ uncertainty on T_c (0.138 days or 0.0025 phase unit for T_c's near BJD 2, 455, 901.361, see Section 3.2). The light curve appears to be highly clustered, or binned, due to the near integral orbital period (P ∼ 55.01 days) and consequent incomplete sampling from a single observing site. Unfortunately, none of the data clusters chance to fall within the predicted transit window, so we are unable to rule out transits of HD 37605b with the APT observations.

Periodogram analysis of the five individual observing seasons revealed no significant periodicity between 1 and 100 days. This suggests that the star is inactive and the observed K ∼ 200 m s⁻¹ RV signal (for HD 37605b) is unlikely to be the result of stellar activity.

Analysis of the complete, normalized data set, however, suggests a weak periodicity of 57.67 ± 0.30 days with a peak-to-peak amplitude of just 0.0012 ± 0.0002 mag (see Figure 6). We tentatively identify this as the stellar rotation period. This period is consistent with the projected rotational velocity of vsin i < 1 km s⁻¹ derived from our stellar analysis described in Section 2.3. It is also consistent with the analysis of Isaacson & Fischer (2010), who derived a Mount Wilson chromospheric Ca ii H and K index of S = 0.165, corresponding to log R'_HK = −5.03. Together, these results imply a rotation period ≳ 46 days and an age of ∼7 Gyr (see Table 2). Similarly, Ibukiyama & Arimoto (2002) find an age of >10 Gyr using isochrones along with the Hipparcos parallax and space motion, supporting HD 37605's low activity and long rotation period.

Zoom In Zoom Out Reset image size

As noted earlier, the near-integer period of HD 37605b makes it difficult to observe from a single longitude. The brightness of the target and the relatively long predicted transit duration creates additional challenges for ground-based observations. We thus observed HD 37605 during 2011 December 5–6 (around the predicted T_c at BJD 2, 455, 901.361 as listed in Table 3) with the MOST satellite launched in 2003 (Walker et al. 2003; Matthews et al. 2004) in the direct imaging mode. This observing technique is similar to ground-based CCD photometry, allowing to apply traditional aperture and point-spread function procedures for data extraction (see, e.g., Rowe et al. 2006 for details). Outlying data points caused by, e.g., cosmic rays were removed.

MOST is orbiting with a period of ∼101 minutes (14.19 cycles day⁻¹, cd⁻¹), which leads to a periodic artifact induced by the scattered light from the earthshine. This signal and its harmonics are further modulated with a frequency of 1 cd⁻¹ originating from the changing albedo of the earth. To correct for this phenomenon, we constructed a cubic fit between the mean background and the stellar flux, which was then subtracted from the data. The reduced and calibrated MOST photometric data are listed in Table 6.

The MOST photometry is shown in Figure 7 for the transit window observations. The vertical dashed lines indicate the beginning and end of the 1σ transit window defined by adding $\sigma _{T_c}$ (0.069 days) on both sides of the predicted transit duration of 0.352 days. The solid line shows the predicted transit model for the previously described planetary parameters. The mean of the relative photometry is 0.00% (or 0.00 mag), with an rms scatter of 0.17%, and within the predicted transit window there are 58 MOST observations. Therefore, the standard error on the mean relative photometry is 0.17%$/\sqrt{58} =$ 0.022%. This means that, for the predicted transit window and a predicted depth of 1.877%, we can conclude a null detection of HD 37605b's transit with extremely high confidence (149σ).

Zoom In Zoom Out Reset image size

Note that the above significance is for an edge-on transit with an impact parameter of b = 0.0. A planetary trajectory across the stellar disk with a higher impact parameter will produce a shorter transit duration. However, the gap between each cluster of MOST measurements is 0.06 days which is 17% of the edge-on transit duration. In order for the duration to be fit within the data gaps, the impact parameter would need to be b > 0.996. To estimate a more conservative lower limit for b, we now assume the most unfortunate case where the transit center falls exactly in the middle of one of the measurement gaps, and also consider the effect of limb darkening by using the nonlinear limb darkening model by Mandel & Agol (2002) with their fitted coefficients for HD 209458. Even under this scenario, we can still conclude the null detection for any transit with b < 0.951 at ≳ 5σ (taking into account that there are at least ∼20 observations will fall within the transit window in this case, though only catching the shallower parts of the transit light curve).

All of the above is based on the assumption that the planet has the predicted radius of 1.1 R_Jup. If in reality the planet is so small that even a b = 0 transit would fall below our detection threshold, it would mean that the planet has a radius of <0.36 R_Jup (a density of >74.50 g cm⁻³), which seems unlikely. It is also very unlikely that our MOST photometry has missed the transit window completely due to an ill-predicted T_c. In the sampling distribution of T_c from BOOTTRAN (with 1000 replicates; see Section 3.2 and the Appendix), there is no T_c that would put the transit window completely off the MOST coverage. In the marginalized posterior distribution of T_c calculated via MCMC (see Section 3.3 and Figure 3), there is only 1 such T_c out of 3 × 10⁴ (0.003%).

The uncertainties listed for the orbital parameter estimates¹⁹ and transit mid-time T_c are calculated via bootstrapping (Freedman 1981; Davison & Hinkley 1997) using the package BOOTTRAN, which we have made publicly available (see Section 3.2). It is designed to calculate error bars for transit ephemerides and the Keplerian orbital fit parameters output by the RVLIN package (Wright & Howard 2009), but can also be a stand-alone package. Thanks to the simple concept of bootstrapping, it is computationally very time efficient and easy to use.

The basic idea of bootstrap is to resample based on original data to create bootstrap samples (multiple data replicates); then for each bootstrap sample, derive orbital parameters or transit parameters through orbital fitting and calculation. The ensemble of parameters obtained in this way yields the approximate sampling distribution for each estimated parameter. The standard deviation of this sampling distribution is the standard error for the estimate.

We caution the readers here that there are regimes in which the "approximate sampling distribution" (a frequentist's concept) is not an estimate of the posterior probability distribution (a Bayesian concept), and there are regimes (e.g., when limited sampling affects the shape of the χ² surface) where there are qualitative differences and the bootstrap method dramatically underestimates uncertainties (e.g., long-period planets when the observations are not yet sufficient to pin down the orbital period; Ford 2005; Bender et al. 2012). In situations with sufficient RV data, good phase coverage, a sufficient time span of observations and a good orbital fit, bootstrap often gives a useful estimate of the parameter uncertainties. For the data considered in this paper, it was not obvious that the bootstrap uncertainty estimate would be accurate, as the time span of observations is only slightly longer than the orbital period of planet c. Nevertheless, we find good agreement between the uncertainty estimates derived from bootstrap and MCMC calculations.

The RV data are denoted as $\lbrace {\bm t}, {\bm v},{\bm \sigma } \rbrace$, where each t_i, v_i, σ_i represents RV v_i observed at time (BJD) t_i with velocity uncertainty σ_i. Extreme outliers should be rejected in order to preserve the validity of our bootstrap algorithm. We first derive our estimates for the true orbital parameters from the original RV data via orbital fitting, using the RVLIN package (Wright & Howard 2009):

$$\begin{equation} {{\bm \beta }} = \mu ( {\bm t}, {\bm v},{\bm \sigma }), \end{equation} \tag{ A1 }$$

where ${\bm \beta }$ is the best fitted orbital parameters.²⁰ From ${\bm \beta }$, we derive $\lbrace {\bm t}, {\bm v}_{{\rm best}}({\bm \beta }) \rbrace$, the best-fit model (here ${\bm t}$ are treated as predictors and thus fixed). Then we can begin resampling to create bootstrap samples.

Our resampling plan is model-based resampling, where we draw from the residuals against the best-fit model. For data that come from the same instrument or telescope, in which case no instrumental offset needs to be taken into account, we simply draw from all residuals, $\lbrace {\bm v}-{\bm v}_{{\rm best}} \rbrace$, with equal probability for each (v_i − v_{best, i}). This new ensemble of residuals, denoted as ${{\bm r}^*}$, is then added to the best-fit model ${\bm v}_{{\rm best}}$ to create one bootstrap sample, ${\bm v}^*$.²¹ Associated with ${{\bm r}^*}$, the uncertainties ${\bm \sigma }$ are also re-assigned to ${\bm v}^*$—that is, if v_j − v_{best, j} is drawn as r_k and added to v_k to generate v*_k, then the uncertainty for v*_k is set to be σ_j.

For data that come from multiple instruments or multiple telescopes, we incorporate our model-based resampling plan to include stratified sampling. In this case, although data from each instrument or telescope are close to homoscedastic, the entire set of data are usually highly heteroscedastic due to stratification in instrument/telescope RV precision. Therefore, the resampling process is done by breaking down the data into different groups, $\lbrace {\bm v}_1,\ {\bm v}_2,\ \ldots \rbrace$, according to instrument and/or telescope, and then resample within each subgroup of data with the algorithm described in last paragraph. The bootstrap sample is then ${\bm v}^* = \lbrace {\bm v}^*_1,\ {\bm v}^*_2,\ \ldots \rbrace$.

To construct the approximate sampling distribution of the orbital parameter estimates ${\bm \beta }$, we compute

$$\begin{equation} {\bm \beta }^* = \mu ({\bm t},\ {\bm v}^*,\ {\bm \sigma }^*) \end{equation} \tag{ A2 }$$

for each bootstrap sample, $\lbrace {\bm t},\ {\bm v}^*,\ {\bm \sigma }^* \rbrace$. The sampling distribution for each orbital parameter estimate β_i can be constructed from the multiple sets of ${\bm \beta }^*$ calculated from multiple bootstrap samples (${\bm \beta }^{*(1)},\ {\bm \beta }^{*(2)},\ \ldots$ from ${\bm v}^{*(1)},\ {\bm v}^{*(2)},\ \ldots$). The standard errors for ${\bm \beta }$ are simply the standard deviations of the sampling distributions.²²

The sampling distribution of the estimated transit mid-time, T_c, is calculated likewise. Here T_c is the transit time for a certain planet of interest in the system, and is usually specified to be the first transit after a designated time T. However, the situation is complicated by the periodic nature of T_c. Our approach is to first calculate, based on the original RV data, T_c0, the estimated mid-time of the first transit after time T₀ (an arbitrary time within the RV observation time window of $[\min ({\bm t}),\ \max ({\bm t})]$; T_c0 is also within this window). Then

$$\begin{equation} T_c = N\cdot P + T_{c0}, \end{equation} \tag{ A3 }$$

where P is the best-estimated period for this planet of interest and N is the smallest integer that is larger than (T − T_c0)/P. Next we compute T*_c0 for each bootstrap sample $\lbrace {\bm t},\ {\bm v}^*,\ {\bm \sigma }^* \rbrace$. Given that within the time window of RV observations ($[\min ({\bm t}),\ \max ({\bm t})]$), the phase of the planet should be known well enough, it is fair to assume that T_c0 is an unbiased estimator of the true transit mid-time. Therefore we assert that T*_c0 has to be well constrained and within the range of [T_c0 − P*/2, T_c0 + P*/2], where P* is the period estimated from this bootstrap sample. If not, then we subtract or add multiple P*'s until T*_c0 falls within the range. Then naturally

$$\begin{equation} T_c^* = N\cdot P^* + T_{c0}^*. \end{equation} \tag{ A4 }$$

The ensemble of T*_c's gives the sampling distribution of T_c and its standard error. Note that T*_c is not necessarily within the rage of [T_c − P/2, T_c + P/2].

Provided with the stellar mass M_⋆ and its uncertainty, we calculate, for each planet in the system, the standard errors for the semi-major axis a and the minimum mass of the planet M_{p, min} (denoted as Msin i in the main text as commonly seen in literature, but this is a somewhat imprecise notation). As the first step, the mass function is calculated for the best-fit ${\bm \beta }$ and each bootstrap sample ${\bm \beta }^*$,

$$\begin{equation} f(P,K,e) = \frac{P K^3 (1-e)^{3/2}}{2\pi G} = \frac{(M_p\cdot \sin i)^3}{(M_\star +M_p)^2}. \end{equation} \tag{ A5 }$$

The sampling distribution of f(P, K, e) then gives the standard error of the mass function. The minimum mass of the planet M_{p, min} is then calculated by assuming sin i = 1 and solving for M_p. Standard error of M_{p, min} is derived through simple propagation of error, as the covariance between M_⋆ and f(P, K, e) is probably negligible.

For the semi-major axis a,

$$\begin{equation} a^3 = \frac{P^2 G (M_\star +M_p)}{4\pi ^2} \approx \frac{P^2 G (M_\star + M_{\rm p,min})}{4\pi ^2}. \end{equation} \tag{ A6 }$$

The standard error of P² is calculated from its bootstrap sampling distribution, and via simple propagation of error we obtain the standard error of a (neglecting covariance between P², M_{p, min}, and M_⋆).