The Discovery of HD 37605 and a Dispositive Null Detection of Transits of HD 3760511affiliation: Based on observations obtained with the Hobby-Eberly Telescope, which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig Maximilians Universität München, and Georg August Universität Göttingen, and observations obtained at the Keck Observatory, which is operated by the University of California. The Keck Observatory was made possible by the generous financial support of the W. M. Keck Foundation.
We report the radial-velocity discovery of a second planetary mass companion to the K0 V star HD 37605, which was already known to host an eccentric, days Jovian planet, HD 37605. This second planet, HD 37605, has a period of years with a low eccentricity and an of . Our discovery was made with the nearly 8 years of radial velocity follow-up at the Hobby-Eberly Telescope and Keck Observatory, including observations made as part of the Transit Ephemeris Refinement and Monitoring Survey (TERMS) effort to provide precise ephemerides to long-period planets for transit follow-up. With a total of 137 radial velocity observations covering almost eight years, we provide a good orbital solution of the HD 37605 system, and a precise transit ephemeris for HD 37605. Our dynamic analysis reveals very minimal planet-planet interaction and an insignificant transit time variation. Using the predicted ephemeris, we performed a transit search for HD 37605 with the photometric data taken by the T12 0.8-m Automatic Photoelectric Telescope (APT) and the Microvariability and Oscillations of Stars (MOST) satellite. Though the APT photometry did not capture the transit window, it characterized the stellar activity of HD 37605, which is consistent of it being an old, inactive star, with a tentative rotation period of days. The MOST photometry enabled us to report a dispositive null detection of a non-grazing transit for this planet. Within the predicted transit window, we exclude an edge-on predicted depth of 1.9% at the level, and exclude any transit with an impact parameter at greater than . We present the BOOTTRAN package for calculating Keplerian orbital parameter uncertainties via bootstrapping. We made a comparison and found consistency between our orbital fit parameters calculated by the RVLIN package and error bars by BOOTTRAN with those produced by a Bayesian analysis using MCMC.
Subject headings:planetary systems — stars: individual (HD 37605)
Jupiter analogs orbiting other stars represent the first signposts of true Solar System analogs, and the eccentricity distribution of these planets with AU will reveal how rare or frequent true Jupiter analogs are. To date, only 9 “Jupiter analogs” have been well-characterized in the peer reviewed literature111HD 13931 (Howard et al., 2010), HD 72659 (Moutou et al., 2011), 55 Cnc (Marcy et al., 2002), HD 134987 (Jones et al., 2010), HD 154345 (Wright et al., 2008, but with possibility of being an activity cycle-induced signal), Ara (Pepe et al., 2007), HD 183263 (Wright et al., 2009), HD 187123 (Wright et al., 2009), and GJ 832 (Bailey et al., 2009). (defined here as years, , and ; Wright et al., 2011, exoplanets.org). As the duration of existing planet searches approach 10–20 years, more and 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 many thousands more candidates have been discovered by the Kepler telescope. Of all of these planets, only three orbit stars with 22255 Cnc (McArthur et al., 2004; Demory et al., 2011), HD 189733 (Bouchy et al., 2005), and HD 209458 (Henry et al., 2000; Charbonneau et al., 2000). and all have 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 2 transiting planets with days have been discovered around stars with , and both have (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 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.
1.2. Initial Discovery and Followup
The inner planet in the system, HD 37605, was the first planet discovered with the Hobby-Eberly Telescope (HET) at McDonald Observatory (Cochran et al., 2004). It is a super Jupiter () on an eccentric orbit with an orbital period in the “period valley” ( 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 radial velocity (RV) residuals (i.e., a none-zero ; 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 37605, and by early 2011 the orbital parameters of the component were clear, and the Texas team was preparing the system for publication.
1.3. TERMS Data
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., 2011a; 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 radial velocity observations of HD 37605 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 radial velocity residuals to the original Cochran et al. (2004) solution. In October 2010 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 January 2011 meeting of the American Astronomical Society (Kane et al., 2011c).
1.4. Synthesis and Outline
In early 2011, the Texas and TERMS teams combined efforts and began joint radial velocity 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 component to be calculated for a thorough transit search. We herein report the transit exclusion of HD 37605 and a stable dynamical solution to the system.
In § 2, we describe our spectroscopic observations and analysis, which provided the radial velocities and the stellar properties of HD 37605. § 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 detection333A dispositive null detection is one that disposes of the question of whether an effect is present, as opposed to one that merely fails to detect a purported or hypothetical effect that may yet lie beneath the detection threshold. The paragon of dispositive null detections is the Michelson-Morley demonstration that the luminiferous ether does not exist (Michelson & Morley, 1887). of non-grazing transits of HD 37605 in § 4. After § 5, Summary and Conclusion, we present updates on of two previously published systems (HD 114762 and HD 168443) in § 6. In the Appendix we describe the algorithm used in the package BOOTTRAN (for calculating orbital parameter error bars; see § 3.2).
2. Spectroscopic Observations and Analysis
2.1. HET and Keck Observations
Observations on HD 37605 at HET started December of 2003. In total, 101 RV observations took place over the course of almost eight 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 next allocated night or after TAC 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 60,000. For more details, see Cochran et al. (2004).
Observations at Keck were taken starting August 2006. A set of 33 observations spanning over five years were made through the HIRES spectrometer (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 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 function444Some authors refer to this as the “point spread function” or the “instrumental profile” of the spectrograph. (IRF) for radial velocity extraction (Marcy & Butler, 1992; Butler et al., 1996). In addition, we also have observations taken without 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 radial velocity extraction and stellar property analysis. The typical working wavelength range for this technique is roughly 5000 Å– 6000 Å.
2.2. Data Reduction and Doppler Analysis
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. (2011b).
We have constructed a complete pipeline for analyzing HET data – from raw data reduction to radial velocity extraction. The raw reduction is done using the REDUCE package by Piskunov & Valenti (2002). This package is designed to optimally extract echelle spectra from 2-D 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 radial velocity extraction for each observation. We have adopted the Doppler code from CPS (e.g. Howard et al., 2009, 2011; Johnson et al., 2011b). 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 TERMS program. We have performed re-reduction on these 44 observations together with all the rest 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 (), 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 m s against best Keplerian fit, much larger than the 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 7.
2.3. Stellar Analysis
|Spectral typeaaESA (1997); van Leeuwen (2008).||K0 V|
|Distance (pc)aaESA (1997); van Leeuwen (2008).||44.0 2.1|
|()||0.901 0.045cc5% relative errors, not the SME intrinsic errors. See footnote 5 for details.|
|()||1.000 0.050cc5% relative errors, not the SME intrinsic errors. See footnote 5 for details.|
|AgebbIsaacson & Fischer (2010), see § 4.1.||Gyr|
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 1 lists the results of our analysis555Note that the errors on the stellar radius and mass listed in Table 1 are not intrinsic to the SME code, but are 5% and 5%. This is because the intrinsic errors reported by SME do not include the errors stemming from the adopted stellar models, and a more realistic precision for and would be around %. Intrinsic errors reported by SME are for and for ., including the effective temperature , surface gravity , iron abundance , projected rotational velocity , bolometric correction BC, bolometric magnitude , stellar luminosity , stellar radius , stellar mass and age. HD 37605 is found to be a metal rich star () with and .
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 , , , , and so on, which are used to derive the bolometric correction BC and 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 and values (along with age and stellar radius). Next, Valenti et al. (2009) introduced an outside loop which re-runs SME with fixed at , followed by another isochrone fit deriving a new log using the updated SME results. The loop continues until values converge. This additional iterative procedure to enforce self-consistency on log is shown to improve the accuracy of other derived stellar parameters (Valenti et al., 2009). The stellar radius and reported here in Table 1 are derived from the final isochrone fit, which are consistent with the purely spectroscopic results. The gravity () is also consistent with the purely spectroscopic gravity () based on strong Mg b damping wings, so for HD 37605 the iteration process is optional.
Cochran et al. (2004) reported the values of , , and Fe/H for HD 37605, and their estimates agree with ours within 1 uncertainty. Santos et al. (2005) also estimated , , Fe/H, and , 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 () is likely km s (corresponding to rotation period 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 , 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 APT data (See §4.1).
3. Orbital Solution
3.1. Transit Ephemeris
The traditional parameters for reporting the ephemerides of spectroscopic binaries are , and , 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 and can be highly covariant. In such cases the circular case is often not excluded by the data, and so the estimation of includes the case , where is undefined. If the best or most likely value of in this case is small but not zero, then it is associated with some nominal value of , but will be very large (approaching ). Since represents the epoch at which the true anomaly equals , will have a similarly large uncertainty (approaching ), 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 lack of knowledge of the covariance in parameters, in particular and (whose covariance is sensitive to the approximate epoch chosen for ). 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 radial velocity data. We used bootstrapping (see Appendix) with the time of conjunction, (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.
3.2. The 37605 System
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 §2.2), we also included six666The velocity from observation on BJD was rejected as it was from a twilight observation, which had both low precision ( m s) and low accuracy (having a residual against the best Keplerian fit of over 100 m s). 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 bootstrap replicates) using the BOOTTRAN package, which is described in detail in the Appendix777The BOOTTRAN package is made publicly available online at http://exoplanets.org/code/ and the Astrophysics Source Code Library..
The best-fit Keplerian parameters are listed in Table 2. The joint Keplerian fit for HD 37605 and HD 37605 has 13 free parameters: the orbital period , time of periastron passage , velocity semi-amplitude , eccentricity , 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 ( and , 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 radial velocity trend in the fit (i.e., fixed to zero). The radial velocity signals and the best Keplerian fits for the system, HD 37605 only, and HD 37605 only are plotted in the three panels of Fig. 1, respectively.
Adopting a stellar mass of (as in Table 1), we estimated the minimum mass () for HD 37605 to be and for HD 37605. While HD 37605 is on a close-in orbit at AU that is highly eccentric (), HD 37605 is found to be on a nearly circular orbit () out at AU, which qualifies it as one of the “Jupiter analogs”.
In order to see whether the period and mass of the outer planet, HD 37605, are well constrained, we mapped out the values for the best Keplerian fit in the - space (subscript ‘’ denoting parameters for the outer planet, HD 37605). Each value on the - grid was obtained by searching for the best-fit model while fixing the period for the outer planet and requiring constraints on and to maintain fixed. As shown in Fig. 2, our data are sufficient to have both and well-constrained. This is also consistent with the tight sampling distributions for and found in our bootstrapping results.
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 ( m s and m s, respectively). This might be due to stellar jitter or underestimated systematic errors in the velocities. We note that the is reduced to if we introduce a stellar jitter of 3.6 m s (added in quadrature to all the RV errors).
|Parameter||HD 37605||HD 37605|
|(days)||55.01307 0.00064||2720 57|
|(BJD)aaTime of Periastron passage.||2453378.241 0.020||2454838 581|
|(BJD)bbTime of conjunction (mid-transit, if the system transits).||2455901.361 0.069|
|(m s)||202.99 0.72||48.90 0.86|
|0.6767 0.0019||0.013 0.015|
|(deg)||220.86 0.28||221 78|
|()||2.802 0.011||3.366 0.072|
|(AU)||0.2831 0.0016||3.814 0.058|
|(m s)||50.7 4.6|
|(m s)ccOffset with respect to the velocities from the 2.1 m telescope.||55.1 4.7|
|(m s)ccOffset with respect to the velocities from the 2.1 m telescope.||36.7 4.7|
|rms (m s)||7.61|
|Jitter (m s)ddIf a jitter of m s is added in quadrature to all RV errors, becomes .||3.6|
|Parameter||HD 37605||HD 37605|
|RVLINBOOTTRAN||MCMCaaMedian values of the marginalized posterior distributions and the 68.27% (‘1’) confidence intervals.||RVLINBOOTTRAN||MCMCaaMedian values of the marginalized posterior distributions and the 68.27% (‘1’) confidence intervals.|
|(days)||55.01307 0.00064||2720 57|
|(BJD)||2453378.243 0.020||2454838 581|
|(m s)||202.99 0.72||48.90 0.86|
|0.6767 0.0019||0.013 0.015|
|(deg)||220.86 0.28||221 78|
|M (deg)bbMean anomaly of the first observation (BJD ).||62.31 0.15||117 78|
|()||2.802 0.011||3.366 0.072|
|(AU)||0.2831 0.0016||3.814 0.058|
|Jitter (m s)ccLike RVLIN, BOOTTRAN assumes no jitter or fixes jitter to a certain value, while MCMC treats it as a free parameter. See § 3.3.||3.6|
3.3. Comparison with MCMC Results
We compared our best Keplerian fit from RVLIN and uncertainties derived from BOOTTRAN (abbreviated as RVLINBOOTTRAN hereafter) with that from a Bayesian framework following Ford (2005) and Ford (2006) (referred to as the MCMC analysis hereafter). Table 3 lists the major orbital parameters from both methods for a direct comparison. Fig. 3 illustrates this comparison, but with the MCMC results presented in terms of 2-D confidence contours for , , , , and of both planets, as well as for of HD 37605.
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 () and jitter (), we adopted a prior of the form , with 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); Ford & Gregory (2007). The likelihood for radial velocity terms assumes that each radial velocity observation () is independent and normally distributed about the true radial velocity with a variance of , where is the published measurement uncertainty. 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 5 Markov chains, each with states. We discarded the first half of the chains and calculate Gelman-Rubin test statistics for each model parameter and several ancillary variables. We found no indications of non-convergence amongst the individual chains. We randomly drew 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 values reported in Table 3. 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 (), 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., 2011a; Payne & Ford, 2011) as the non-Keplerian perturbations should be tiny (detail on the magnitude of the perturbations is provided in §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 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 §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 3 and Fig. 3, the parameter estimates from RVLINBOOTTRAN 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. m s) unless one goes to an extremely low jitter value (e.g. 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 and semi-major axes , the MCMC analysis does not incorporate the errors on the stellar mass estimate. Note here, as previously mentioned in § 3.1, that the “best-fit” parameters reported by the MCMC analysis here listed in Table 3 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 , while the MCMC parameters listed in Table 3 give a higher value of 2.91.
3.4. Dynamical Analysis
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: (i) are long-term stable; (ii) 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); (iii) 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 AU will have companion orbits which are Hill stable for AU (Gladman, 1993), so while Hill stability does not preclude outward scatter of the outer planet, the fact that AU suggests that the system will be far from any such instability.
We integrated the systems for years ( the orbital period of the outer planet and the secular period of the system), and plot in Fig. 4 the evolution of the orbital elements , & . On the left-hand side of the plot we provide short-term detail, illustrating that over the 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 out of years of system evolution, demonstrating that (i) the secular variation in some of the elements (particularly the eccentricity of the outer planet; see in red) over a time span of years can be significant: in this case we see , but (ii) the system appears completely stable, as one would expect for planets with a period ratio . Finally, at the bottom of the figure we display the range of parameter space covered by the parameters ( in blue for inner planet and 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.
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 year simulation of a system in which , pushing the planetary masses to . Even in such a pathological system the eccentricity oscillations are only increased by a factor of 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 s, or day, which is much smaller than the error bar on ( day). Therefore we did not take into account the effect of TTV when performing our transit analysis in the next section.
3.5. Activity Cycles and Jupiter Analogs
The coincidence of the Solar activity cycle period of 11 years and Jupiter’s orbital period near 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 radial velocity 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 radial velocity 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 ), but that this star also shows activity variations in phase with the radial velocity variations. They noted that many Sun-like stars, including the precise radial velocity 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 radial velocity variations (Santos et al., 2010; Lovis et al., 2011; 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 radial velocity in the average measurements of a sample of HARPS stars, and 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 semiamplitude 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 Fig. 5 plots the T12 APT observations from all five observing seasons (data provided in Table 4; see details on APT photometry in § 4.1). The dashed line marks the mean relative magnitude () of the first season. The seasonal mean brightness of the star increases gradually from year to year by a total of 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 & 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 m s.
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.
4. The Dispositive Null Detection of Transits of HD 37605
We have performed a transit search for the inner planet of the system, HD 37605. This planet has a transit probability of 1.595% and a predicted transit duration of 0.352 day, as derived from the stellar parameters listed in Table 1 and the orbital parameters given in Table 2. From the minimum planet mass (; see Table 2) and the models of Bodenheimer et al. (2003), we estimate its radius to be . Combined with the stellar radius of HD 37605 listed in Table 1, , we estimate the transit depth to be 1.877% (for an edge-on transit, ). We used both ground-based (APT; §4.1) and space-based (MOST; §4.2) facilities in our search.
4.1. APT Observations and Analysis
The T12 0.8-m Automatic Photoelectric Telescope (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 , the mean of the differential magnitudes acquired simultaneously in the Strömgren and 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, BV = 0.996, K0 III), HD 38145 (V = 7.89, BV = 0.326, F0 V), and HD 38779 (V = 7.08, BV = 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).
|Heliocentric Julian Date|
Note. – This table is presented in its entirety in the electronic edition of the Astrophysical Journal. A portion is shown here for guidance regarding its form and content.
Fig. 5 illustrates the APT photometric data and our transit search. As mentioned in § 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 37605, 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 37605, days, and the predicted time of mid-transit, , defined as Phase 0. The normalized and phased data are plotted in the middle panel of Fig. 5. The solid line is the predicted transit light curve, with the predicted transit duration (0.352 day or 0.0064 phase unit) and transit depth (1.877% or mag) as estimated above. The scatter of the phased data from their mean is 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 five years. A least-squares sine fit of the phased data gives a very small semi-amplitude of 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 Fig. 5 plots the phased data around the predicted time of mid-transit, , at an expanded scale on the abscissa. The horizontal bar below the transit window represents the 1 uncertainty on (0.138 day or 0.0025 phase unit for ’s near BJD ; see § 3.2). The light curve appears to be highly clustered, or binned, due to the near integral orbital period ( 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 m s RV signal (for HD 37605) is unlikely to be the result of stellar activity.
Analysis of the complete, normalized data set, however, suggests a week periodicity of days with a peak-to-peak amplitude of just mag (see Fig. 6). We tentatively identify this as the stellar rotation period. This period is consistent with the projected rotational velocity of km sderived from our stellar analysis described in §2.3. It is also consistent with the analysis of Isaacson & Fischer (2010), who derived a Mount Wilson chromospheric Ca ii H & K index of , corresponding to . Together, these results imply a rotation period days and an age of Gyr (see Table 1). Similarly, Ibukiyama & Arimoto (2002) find an age of Gyr using isochrones along with the Hipparcos parallax and space motion, supporting HD 37605’s low activity and long rotation period.
4.2. MOST Observations and Analysis
|Heliocentric Julian Date||Relative Magnitude|
Note. – This table is presented in its entirety in the electronic edition of the Astrophysical Journal. A portion is shown here for guidance regarding its form and content.
As noted earlier, the near-integer period of HD 37605 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 at BJD as listed in Table 2) with the MOST (Microvariability and Oscillations of Stars) 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 PSF 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 per 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 5.
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 (0.069 day) 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 a 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% 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 37605’s transit with extremely high confidence (149).
Note that the above significance is for an edge-on transit with an impact parameter of . 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 . To estimate a more conservative lower limit for , 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 non-linear 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 at (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 . If in reality the planet is so small that even a transit would fall below our detection threshold, it would mean that the planet has a radius of (a density of 74.50 gcm), which seems unlikely. It is also very unlikely that our MOST photometry has missed the transit window completely due to an ill-predicted . In the sampling distribution of from BOOTTRAN (with 1000 replicates; see § 3.2 and Appendix), there is no that would put the transit window completely off the MOST coverage. In the marginalized posterior distribution of calculated via MCMC (see § 3.3 and Fig. 3), there is only 1 such out of (0.003%).
5. Summary and Conclusion
In this paper, we report the discovery of HD 37605 and the dispositive null detection of non-grazing transits of HD 37605, the first planet discovered by HET. HD 37605 is the outer planet of the system with a period of years on a nearly circular orbit () at AU. It is a “Jupiter analog” with , which adds one more sample to the currently still small inventory of such planets (only 10 including HD 37605; see §1). The discovery and characterization of “Jupiter analogs” will help understanding the formation of gas giants as well as the frequency of true solar system analogs. This discovery is a testimony to the power of continued observation of planet-bearing stars.
Using our RV data with nearly 8-year long baseline, we refined the orbital parameters and transit ephemerides of HD 37605. The uncertainty on the predicted mid-transit time was constrained down to 0.069 day (at and near in BJD), which is small compared to the transit duration (0.352 day). In fact, just the inclusion of the two most recent points in our RV data have reduced the uncertainty on by over 10%. We have performed transit search with APT and the MOST satellite. Because of the near-integer period of HD 37605 and the longitude of Fairborn Observatory, the APT photometry was unable to cover the transit window. However, its excellent photometric precision over five observing seasons enabled us to rule out the possibility of the RV signal being induced by stellar activity. The MOST photometric data, on the other hand, were able to rule out an edge-on transit with a predicted depth of 1.877% at a 10 level, with a lower limit on the impact parameter of . This transit exclusion is a further demonstration of the TERMS strategy, where follow-up RV observations help to reduce the uncertainty on transit timing and enable transit searches.
Our best-fit orbital parameters and errors from RVLINBOOTTRAN were found to be consistent with those derived from a Bayesian analysis using MCMC. Based on the best-fit MCMC systems, we performed dynamic and TTV analysis on the HD 37605 system. Dynamic analysis shows no sign of orbital resonance and very minimal planet-planet interaction. We derived a TTV of 100 s, which is much smaller than .
We have also performed a stellar analysis on HD 37605, which shows that it is a metal rich star () with a stellar mass of with a radius of . The small variation seen in our photometric data (amplitude mag over the course of four years) suggests that HD 37605 is consistent as being an old, inactive star that is probably slowly rotating. We tentatively propose that the rotation period of the star is days, based on a weak periodic signal seen in our APT photometry.
6. Note on Previously Published Orbital Fits
In early 2012, we repaired a minor bug in the BOOTTRAN package, mostly involving the calculation and error bar estimation of . As a result, the values and their errors for two previously published systems (three planets) need to be updated. They are: HD 114762 (Kane et al., 2011a), HD 168443, and HD 168443 (Pilyavsky et al., 2011b). Table 6 lists the updated and error bars.
One additional system, HD 63454 (Kane et al., 2011d), was also analyzed using BOOTTRAN. However, the mass of HD 63454 is small enough compared to its host mass and thus was not affected by this change.
|Planet||std. error ()|
|HD 114762aaFor best orbital fit with RV trend ().|
|HD 114762bbFor best orbital fit without RV trend ().|
- Bailey et al. (2009) Bailey, J., Butler, R. P., Tinney, C. G., Jones, H. R. A., O’Toole, S., Carter, B. D., & Marcy, G. W. 2009, ApJ, 690, 743
- Baliunas et al. (1995) Baliunas, S. L., et al. 1995, ApJ, 438, 269
- Baluev (2009) Baluev, R. V. 2009, MNRAS, 393, 969
- Barbieri et al. (2007) Barbieri, M., et al. 2007, A&A, 476, L13
- Bender et al. (2012) Bender, C. F., et al. 2012, ApJ, 751, L31
- Bodenheimer et al. (2003) Bodenheimer, P., Laughlin, G., & Lin, D. N. C. 2003, ApJ, 592, 555
- Boisse et al. (2012) Boisse, I., et al. 2012, A&A, 545, A55
- Bouchy et al. (2005) Bouchy, F., et al. 2005, A&A, 444, L15
- Butler et al. (1996) Butler, R. P., Marcy, G. W., Williams, E., McCarthy, C., Dosanjh, P., & Vogt, S. S. 1996, PASP, 108, 500
- Campbell et al. (1991) Campbell, B., Yang, S., Irwin, A. W., & Walker, G. A. H. 1991, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 390, Bioastronomy: The Search for Extraterrestial Life — The Exploration Broadens, ed. J. Heidmann & M. J. Klein, 19–20
- Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
- Charbonneau et al. (2000) Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45
- Cochran et al. (2004) Cochran, W. D., et al. 2004, ApJ, 611, L133
- Davison & Hinkley (1997) Davison, A., & Hinkley, D. 1997, Bootstrap Methods and Their Application, Cambridge Series on Statistical and Probabilistic Mathematics (Cambridge University Press)
- Demarque et al. (2004) Demarque, P., Woo, J.-H., Kim, Y.-C., & Yi, S. K. 2004, ApJS, 155, 667
- Deming et al. (1987) Deming, D., Espenak, F., Jennings, D. E., Brault, J. W., & Wagner, J. 1987, ApJ, 316, 771
- Demory et al. (2011) Demory, B.-O., et al. 2011, A&A, 533, A114
- Dragomir et al. (2011) Dragomir, D., et al. 2011, AJ, 142, 115
- Dravins (1985) Dravins, D. 1985, in Stellar Radial Velocities, ed. A. G. D. Philip & D. W. Latham, 311–320
- Dumusque et al. (2011) Dumusque, X., et al. 2011, A&A, 535, A55
- Dumusque et al. (2012) —. 2012, Nature, advanced online publication, doi:10.1038/nature11572
- ESA (1997) ESA. 1997, VizieR Online Data Catalog, 1239, 0
- Fischer et al. (2007) Fischer, D. A., et al. 2007, ApJ, 669, 1336
- Ford (2005) Ford, E. B. 2005, AJ, 129, 1706
- Ford (2006) —. 2006, ApJ, 642, 505
- Ford & Gregory (2007) Ford, E. B., & Gregory, P. C. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 371, Statistical Challenges in Modern Astronomy IV, ed. G. J. Babu & E. D. Feigelson, 189
- Fossey et al. (2009) Fossey, S. J., Waldmann, I. P., & Kipping, D. M. 2009, MNRAS, 396, L16
- Freedman (1981) Freedman, D. A. 1981, The Annals of Statistics, 9, pp. 1218
- Gladman (1993) Gladman, B. 1993, Icarus, 106, 247
- Henry (1999) Henry, G. W. 1999, PASP, 111, 845
- Henry et al. (2000) Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000, ApJ, 529, L41
- Hoaglin et al. (1983) Hoaglin, D. C., Mosteller, F., & Tukey, J. W. 1983, Understanding Robust and Exploratory Data Anlysis, Wiley Classics Library Editions (Wiley)
- Horne (1986) Horne, K. 1986, PASP, 98, 609
- Howard et al. (2009) Howard, A. W., et al. 2009, ApJ, 696, 75
- Howard et al. (2010) —. 2010, ApJ, 721, 1467
- Howard et al. (2011) —. 2011, ApJ, 726, 73
- Ibukiyama & Arimoto (2002) Ibukiyama, A., & Arimoto, N. 2002, A&A, 394, 927
- Isaacson & Fischer (2010) Isaacson, H., & Fischer, D. 2010, ApJ, 725, 875
- Johnson et al. (2011a) Johnson, J. A., et al. 2011a, AJ, 141, 16
- Johnson et al. (2011b) —. 2011b, ApJS, 197, 26
- Jones et al. (2010) Jones, H. R. A., et al. 2010, MNRAS, 403, 1703
- Kane et al. (2012) Kane, S. R., Ephemeris Refinement, T., & Survey(TERMS), M. 2012, in American Astronomical Society Meeting Abstracts, Vol. 219, American Astronomical Society Meeting Abstracts, #228.07
- Kane et al. (2011a) Kane, S. R., Henry, G. W., Dragomir, D., Fischer, D. A., Howard, A. W., Wang, X., & Wright, J. T. 2011a, ApJ, 735, L41
- Kane et al. (2009) Kane, S. R., Mahadevan, S., von Braun, K., Laughlin, G., & Ciardi, D. R. 2009, PASP, 121, 1386
- Kane et al. (2011b) Kane, S. R., et al. 2011b, Detection and Dynamics of Transiting Exoplanets, St. Michel l’Observatoire, France, Edited by F. Bouchy; R. Díaz; C. Moutou; EPJ Web of Conferences, Volume 11, id.06005, 11, 6005
- Kane et al. (2011c) Kane, S. R., et al. 2011c, in American Astronomical Society Meeting Abstracts, #128.03
- Kane et al. (2011d) —. 2011d, ApJ, 737, 58
- Laughlin et al. (2009) Laughlin, G., Deming, D., Langton, J., Kasen, D., Vogt, S., Butler, P., Rivera, E., & Meschiari, S. 2009, Nature, 457, 562
- Lovis et al. (2011) Lovis, C., et al. 2011, A&A submitted
- Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171
- Marcy & Butler (1992) Marcy, G. W., & Butler, R. P. 1992, PASP, 104, 270
- Marcy et al. (2002) Marcy, G. W., Butler, R. P., Fischer, D. A., Laughlin, G., Vogt, S. S., Henry, G. W., & Pourbaix, D. 2002, ApJ, 581, 1375
- Matthews et al. (2004) Matthews, J. M., Kuschnig, R., Guenther, D. B., Walker, G. A. H., Moffat, A. F. J., Rucinski, S. M., Sasselov, D., & Weiss, W. W. 2004, Nature, 430, 51
- McArthur et al. (2004) McArthur, B. E., et al. 2004, ApJ, 614, L81
- Michelson & Morley (1887) Michelson, A. A., & Morley, E. 1887, American Journal of Science, 34, 333
- Moutou et al. (2011) Moutou, C., et al. 2011, A&A, 527, A63
- Payne & Ford (2011) Payne, M. J., & Ford, E. B. 2011, ApJ, 729, 98
- Pepe et al. (2007) Pepe, F., et al. 2007, A&A, 462, 769
- Pilyavsky et al. (2011a) Pilyavsky, G., et al. 2011a, ApJ, 743, 162
- Pilyavsky et al. (2011b) —. 2011b, ApJ, 743, 162
- Piskunov & Valenti (2002) Piskunov, N. E., & Valenti, J. A. 2002, A&A, 385, 1095
- Rowe et al. (2006) Rowe, J. F., et al. 2006, ApJ, 646, 1241
- Santos et al. (2010) Santos, N. C., Gomes da Silva, J., Lovis, C., & Melo, C. 2010, A&A, 511, A54
- Santos et al. (2005) Santos, N. C., Israelian, G., Mayor, M., Bento, J. P., Almeida, P. C., Sousa, S. G., & Ecuvillon, A. 2005, A&A, 437, 1127
- Schneider et al. (2011) Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79
- Shetrone et al. (2007) Shetrone, M., et al. 2007, PASP, 119, 556
- Torres et al. (2010) Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67
- Tull (1998) Tull, R. G. 1998, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 3355, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. S. D’Odorico, 387–398
- Valenti & Fischer (2005) Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141
- Valenti & Piskunov (1996) Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595
- Valenti et al. (2009) Valenti, J. A., et al. 2009, 702, 989
- van Leeuwen (2008) van Leeuwen, F. 2008, VizieR Online Data Catalog, 1311, 0
- Vogt et al. (1994) Vogt, S. S., et al. 1994, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 2198, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. D. L. Crawford & E. R. Craine, 362
- Walker et al. (2003) Walker, G., et al. 2003, PASP, 115, 1023
- Wittenmyer et al. (2007) Wittenmyer, R. A., Endl, M., Cochran, W. D., & Levison, H. F. 2007, AJ, 134, 1276
- Wittenmyer et al. (2012) Wittenmyer, R. A., et al. 2012, ApJ, 753, 169
- Wright & Howard (2009) Wright, J. T., & Howard, A. W. 2009, ApJS, 182, 205
- Wright et al. (2008) Wright, J. T., Marcy, G. W., Butler, R. P., Vogt, S. S., Henry, G. W., Isaacson, H., & Howard, A. W. 2008, ApJ, 683, L63
- Wright et al. (2009) Wright, J. T., Upadhyay, S., Marcy, G. W., Fischer, D. A., Ford, E. B., & Johnson, J. A. 2009, ApJ, 693, 1084
- Wright et al. (2011) Wright, J. T., et al. 2011, PASP, 123, 412
Uncertainties via Bootstrapping
The uncertainties listed for the orbital parameter estimates888Through out the paper and sometimes in this Appendix, we refer to the “estimates of the parameters” (as distinguished from the “true parameters”, which are not known and can only be estimated) simply as the “parameters”. and transit mid-time are calculated via bootstrapping (Freedman, 1981; Davison & Hinkley, 1997) using the package BOOTTRAN, which we have made publicly available (see § 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 . Nevertheless, we find good agreement between the uncertainty estimates derived from bootstrap and MCMC calculations.
The radial velocity data are denoted as , where each , , represents radial velocity observed at time (BJD) with velocity uncertainty . 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):
where is the best fitted orbital parameters999As described in §3.2, this includes the , , , , and for each planet, as well as , (if applicable), and velocity offsets between instruments/telescopes (if applicable) for the system.. From , we derive , the best-fit model (here 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, , with equal probability for each . This new ensemble of residuals, denoted as , is then added to the best-fit model to create one bootstrap sample, 101010We simply use the raw residual instead of any form of modified residual, because the RV data for any single instrument or telescope are usually close enough to homoscedasticity.. Associated with , the uncertainties are also re-assigned to – that is, if is drawn as and added to to generate , then the uncertainty for is set to be .
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 radial velocity precision. Therefore, the resampling process is done by breaking down the data into different groups, , 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 .
To construct the approximate sampling distribution of the orbital parameter estimates , we compute
for each bootstrap sample, . The sampling distribution for each orbital parameter estimate can be constructed from the multiple sets of calculated from multiple bootstrap samples ( from ). The standard errors for are simply the standard deviations of the sampling distributions111111The standard deviation of a sampling distribution is estimated in a robust way using the IDL function robust_sigma, which is written by H. Fruedenreich based on the principles of robust estimation outlined in Hoaglin et al. (1983)..
The sampling distribution of the estimated transit mid-time, , is calculated likewise. Here 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 . However, the situation is complicated by the periodic nature of . Our approach is to first calculate, based on the original RV data, , the estimated mid-time of the first transit after time (an arbitrary time within the RV observation time window of ; is also within this window). Then
where is the best-estimated period for this planet of interest, and is the smallest integer that is larger than . Next we compute for each bootstrap sample . Given that within the time window of radial velocity observations (), the phase of the planet should be known well enough, it is fair to assume that is an unbiased estimator of the true transit mid-time. Therefore we assert that has to be well constrained and within the range of , where is the period estimated from this bootstrap sample. If not, then we subtract or add multiple ’s until falls within the range. Then naturally
The ensemble of ’s gives the sampling distribution of and its standard error. Note that is not necessarily within the rage of .
Provided with the stellar mass and its uncertainty, we calculate, for each planet in the system, the standard errors for the semi-major axis and the minimum mass of the planet (denoted as 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 and each bootstrap sample ,
The sampling distribution of then gives the standard error of the mass function. The minimum mass of the planet is then calculated by assuming and solving for . Standard error of is derived through simple propagation of error, as the covariance between and is probably negligible.
For the semi-major axis ,
The standard error of is calculated from its bootstrap sampling distribution, and via simple propagation of error we obtain the standard error of (neglecting covariance between , , and ).
|BJD||(m s)||(m s)||Telescope|