HAT-P-13b,c: a transiting hot Jupiter with a massive outer companion on an eccentric orbit
G. A. Bakos, A. W. Howard, R. W. Noyes, J. Hartman, G. Torres, Geza Kovacs, D. A. Fischer, D. W. Latham, J. A. Johnson, G. W. Marcy, D. D. Sasselov, R. P. Stefanik, B. Sipocz, Gabor Kovacs, G. A. Esquerdo, A. Pal, J. Lazar, I. Papp
aa r X i v : . [ a s t r o - ph . E P ] O c t ApJ, accepted
Preprint typeset using L A TEX style emulateapj v. 05/04/06
HAT-P-13b,c: A TRANSITING HOT JUPITER WITH A MASSIVE OUTER COMPANION ON AN ECCENTRICORBIT † G. ´A. Bakos , A. W. Howard , R. W. Noyes , J. Hartman , G. Torres , G´eza Kov´acs , D. A. Fischer ,D. W. Latham , J. A. Johnson , G. W. Marcy , D. D. Sasselov , R. P. Stefanik , B. Sip˝ocz , G´abor Kov´acs ,G. A. Esquerdo , A. P´al , J. L´az´ar , I. Papp , P. S´ari ApJ, accepted
ABSTRACTWe report on the discovery of a planetary system with a close-in transiting hot Jupiter on a nearcircular orbit and a massive outer planet on a highly eccentric orbit. The inner planet, HAT-P-13b,transits the bright V=10.622 G4 dwarf star GSC 3416-00543 every P = 2 . ± . T c = 2454779 . ± . . ± . P = 428 . ± . T c = 2454870 . ± . T ,peri = 2454890 . ± . . +0 . − . M ⊙ , radius of 1 . ± . R ⊙ , effective temperature 5653 ±
90 K, and is rathermetal rich with [Fe / H] = +0 . ± .
08. The inner planetary companion has a mass of 0 . +0 . − . M J ,and radius of 1 . ± . R J yielding a mean density of 0 . +0 . − . g cm − . The outer companion has m sin i = 15 . ± . M J , and orbits on a highly eccentric orbit of e = 0 . ± . Subject headings: planetary systems — stars: individual (HAT-P-13, GSC 3416-00543) techniques:spectroscopic, photometric INTRODUCTION
Radial velocity (RV) surveys have shown that multiple-planet stellar systems are common. For example, Wrightet al. (2007) concluded that the occurrence of additionalplanets among stars already having one known planetmust be greater than 30%. Thus, one would expect thatout of the ∼
50 published transiting extrasolar planet(TEP) systems , there should be a number of systemswith additional companions; these companions shouldmake their presence known through radial velocity vari-ations of the parent stars, even if they do not them-selves transit. The fact that so far no TEPs in multipleplanet systems have been reported is somewhat surpris-ing based on the statistics from RV surveys (see also Fab-rycky 2009). Recently Smith et al. (2009) searched thelight curves of 24 transiting planets for outer compan- Harvard-Smithsonian Center for Astrophysics, Cambridge,MA, [email protected] NSF Fellow Department of Astronomy, University of California, Berkeley,CA Konkoly Observatory, Budapest, Hungary Department of Physics and Astronomy, San Francisco StateUniversity, San Francisco, CA Institute for Astronomy, University of Hawaii, Honolulu, HI96822; NSF Postdoctoral Fellow Department of Astronomy, E¨otv¨os Lor´and University, Bu-dapest, Hungary. Hungarian Astronomical Association, Budapest, Hungary † Based in part on observations obtained at the W. M. Keck Ob-servatory, which is operated by the University of California and theCalifornia Institute of Technology. Keck time has been granted byNOAO (A146Hr,A264Hr) and NASA (N128Hr,N145Hr). ions, finding no evidence for a double transiting system.The Hungarian-made Automated Telescope Network(HATNet) survey (Bakos et al. 2002, 2004) has beena major contributor to the discovery of TEPs. Oper-ational since 2003, it has covered approximately 10% ofthe Northern sky, searching for TEPs around bright stars(8 . I . . § § § § § § § OBSERVATIONS
Photometric detection
The transits of HAT-P-13b were detected with theHAT-5 telescope. The region around GSC 3416-00543, afield internally labeled as 136, was observed on a nightlybasis between 2005 November 25 and 2006 May 20,whenever weather conditions permitted. We gathered4021 exposures of 5 minutes at a 5.5-minute cadence.Each image contained approximately 20000 stars down to I ∼ .
0. For the brightest stars in the field we achieveda per-image photometric precision of 3.1 mmag. I - band m agn i t ude Orbital phase
Fig. 1.—
The unbinned light curve of HAT-P-13 including all4021 instrumental I band measurements obtained with the HAT-5telescope of HATNet (see text for details), and folded with the pe-riod of P = 2 . § The calibration of the HATNet frames was done uti-lizing standard procedures. The calibrated frames werethen subjected to star detection and astrometry, as de-scribed in P´al & Bakos (2006). Aperture photometrywas performed on each image at the stellar centroids de-rived from the 2MASS catalog (Skrutskie et al. 2006)and the individual astrometrical solutions. Descriptionof numerous details related to the reduction are alsogiven in P´al (2009). The resulting light curves weredecorrelated against trends using the External Parame-ter Decorrelation technique in “constant” mode (EPD,see Bakos et al. 2009) and the Trend Filtering Algorithm(TFA, see Kov´acs et al. 2005). The light curves weresearched for periodic box-like signals using the Box LeastSquares method (BLS, see Kov´acs et al. 2002). We de-tected a significant signal in the light curve of GSC 3416-00543 (also known as 2MASS 08393180+4721073; α =08 h m . δ = +47 ◦ ′ . ′′ ; J2000), with a depthof ∼ . P = 2 . q ≈ . ± . P q ≈ . ± .
040 hours (see Fig. 1).
Reconnaissance Spectroscopy
One of the important tools used in our survey for estab-lishing whether the transit-feature in the light curve of acandidate is due to astrophysical phenomena other thana planet transiting a star is the CfA Digital Speedome-ter (DS; Latham 1992), mounted on the FLWO 1.5 mtelescope. This yields high-resolution spectra with lowsignal-to-noise ratio sufficient to derive radial velocitieswith moderate precision (roughly 1 km s − ), and to de-termine the effective temperature and surface gravity ofthe host star. With this facility we are able to weedout certain false alarms, such as eclipsing binaries andmultiple stellar systems.We obtained 7 spectra spanning 37 days with the DS.The RV measurements of HAT-P-13 showed an rms resid-ual of 0 .
68 km s − , consistent with no detectable RV vari-ation. The spectra were single-lined, showing no evidencefor more than one star in the system. Atmospheric pa-rameters for the star, including the initial estimates ofeffective temperature T eff ⋆ = 5500 ± K , surface grav-ity log g ⋆ = 4 . ± .
25 (log cgs), and projected rotationalvelocity v sin i = 0 . ± . − , were derived as de-scribed by Torres et al. (2002). The effective tempera-ture and surface gravity correspond to a G4 dwarf (Cox2000). The mean line-of-sight velocity of HAT-P-13 is+14 . ± .
68 km s − . High resolution, high S/N spectroscopy and thesearch for radial velocity signal components
Given the significant transit detection by HATNet,and the encouraging DS results, we proceeded with thefollow-up of this candidate by obtaining high-resolutionand high S/N spectra to characterize the radial velocityvariations and to determine the stellar parameters withhigher precision. Using the HIRES instrument (Vogt etal. 1994) on the Keck I telescope located on Mauna Kea,Hawaii, we obtained 30 exposures with an iodine cell,plus three iodine-free templates. The first template hadlow signal-to-noise ratio, thus we repeated the templateobservations during a later run, and acquired two highquality templates. We used the last two templates inthe analysis. The observations were made between 2008March 22 and 2009 June 5. The relative radial velocity(RV) measurements are listed in Tab. 1, and shown inFig. 2.Observations and reductions have been carried out inan identical way to that described in earlier HATNetdiscovery papers, such as Bakos et al. (2009). Referencesfor the Keck iodine cell observations and the reduction ofthe radial velocities are given in Marcy & Butler (1992)and Butler et al. (1996).Initial fits of these RVs to a single-planet Keplerianorbit were quite satisfactory but soon revealed a slightresidual trend that became more significant with time.This is the reason for the observations extending overmore than one year, as opposed to just a few months asnecessary to confirm simpler purely sinusoidal variationsseen in other transiting planets discovered by HATNet.Eventually we noticed that the residual trend reversed(Fig. 2, top), a clear sign of a coherent motion most likelydue to a more distant body in the same system, possi-bly a massive planet. Preliminary two-planet orbital so-lutions provided a much improved fit (although with aAT-P-13b,c 3
TABLE 1Relative radial velocity, bisector and activity index measurements ofHAT-P-13.
BJD RV a σ RVb BS σ BS S c σ S (2,454,000+) (m s − ) (m s − ) (m s − ) (m s − )547 . .
42 4 .
30 1 .
37 5 .
85 0 . . . · · · · · · − .
90 7 .
27 0 . . . .
80 1 .
43 12 .
35 5 .
35 0 . . . .
56 1 . − .
03 6 .
73 0 . . . .
48 1 .
35 6 .
95 5 .
76 0 . . . .
33 2 . − .
77 6 .
52 0 . . . .
86 1 . − .
24 7 .
08 0 . . . .
20 1 . − .
20 7 .
62 0 . . . .
01 2 . − .
47 7 .
63 0 . . . .
23 1 . − .
98 6 .
98 0 . . . .
53 1 .
81 13 .
22 5 .
05 0 . . . .
36 1 .
64 10 .
16 5 .
60 0 . . . .
19 1 . − .
61 6 .
46 0 . . . .
58 1 .
72 6 .
92 5 .
77 0 . . . .
11 1 .
59 11 .
80 5 .
04 0 . . . .
44 1 . − .
17 6 .
67 0 . . . − .
92 2 .
33 0 .
27 6 .
54 0 . . . .
29 3 . − .
44 6 .
77 0 . . . − .
84 1 . − .
53 7 .
60 0 . . . − .
51 1 .
41 0 .
09 6 .
71 0 . . . − .
07 2 . − .
77 9 .
37 0 . . . − .
44 1 .
35 4 .
00 6 .
19 0 . . . − .
69 2 . − .
31 8 .
75 0 . . . − .
43 1 . − .
92 8 .
19 0 . . . .
55 1 . − .
32 7 .
22 0 . . . − .
24 1 . − .
27 7 .
21 0 . . . .
67 1 .
47 4 .
02 6 .
21 0 . . . − .
66 1 .
40 6 .
16 5 .
92 0 . . . .
74 1 .
39 5 .
20 6 .
05 0 . . . · · · · · · .
77 5 .
66 0 . . . · · · · · · .
37 5 .
88 0 . . . .
20 1 .
73 3 .
95 6 .
31 0 . . . .
15 1 . − .
87 7 .
22 0 . . a The fitted zero-point that is on an arbitrary scale (denoted as γ rel in Tab. 4)has been subtracted from the velocities. b The values for σ RV do not include the jitter. c S values are on a relative scale. weakly determined outer period due to the short dura-tion of the observations), and more importantly, held upafter additional observations. With the data so far gath-ered and presented in this work, a false alarm probabil-ity (FAP) analysis for the Keplerian orbit of the secondplanet yielded an FAP of 0.00001.A more thorough analysis using a two-component leastsquares Fourier fit (see Barning 1963, for the single-component case), with one component fixed to the knownfrequency of the short-period inner planet, also confirmedthe existence of the long period component, and indi-cated that the latter is due to a highly non-sinusoidalmotion. A number of other low frequency peaks wereeliminated as being aliases yielding worse quality fits. Athree-harmonic representation of the fit yielded an RMSof 5 m s − , fairly close to the value expected from theformal errors. Full modeling of this complex motionof two Keplerian orbits superposed, simultaneously withthe photometry, is described in our global fit of § Photometric follow-up observations
We observed 7 transit events of HAT-P-13 with theKeplerCam CCD on the FLWO 1.2 m telescope between2008 April 24 and 2009 May 8. All observations werecarried out in the i band, and the typical exposure timewas 15 seconds. The reduction of the images, including calibration, astrometry, and photometry, was performedas described in Bakos et al. (2009).We performed EPD and TFA against trends simulta-neously with the light curve modeling (for more details,see § § ANALYSIS
Properties of the parent star
We derived the initial stellar atmospheric parametersby using the first template spectrum obtained with theKeck/HIRES instrument. We used the SME package ofValenti & Piskunov (1996) along with the atomic-linedatabase of Valenti & Fischer (2005), which yielded thefollowing initial values and uncertainties (which we haveconservatively increased, to include our estimates of thesystematic errors): effective temperature T eff ⋆ = 5760 ±
90 K, stellar surface gravity log g ⋆ = 4 . ± .
10 (cgs),metallicity [Fe / H] = +0 . ± .
08 dex, and projected ro-tational velocity v sin i = 1 . ± . − .Further analysis of the spectra has shown that the hoststar is chromospherically quiet with log R ′ HK = − . S = 0 .
14 (on an absolute scale). Bakos et al. -1000-800-600-400-200 0 200 400 550 600 650 700 750 800 850 900 950 1000 R V ( m / s ) -10-5 0 5 10 15 550 600 650 700 750 800 850 900 950 1000 R e s i dua l BJD - 2454000-800-600-400-200 0 200 R V ( m / s ) Outer phase with respect to T c2 -100-50 0 50 100 R V ( m / s ) Inner phase with respect to T c1 -60-40-20 0 20 40 B i s e c . ( m / s ) Inner phase with respect to T c1 S a c t i v i t y Inner phase with respect to T c1 Fig. 2.— (Top) Radial-velocity measurements from Keck forHAT-P-13, along with the best 2-planet orbital fit, shown as afunction of BJD (see § § − , requiring a jitter of 3 . − to be added in quadratureto the individual errors to yield a reduced χ of 1.0. The error-barsin this panel have been inflated accordingly. Note the different ver-tical scale of the panels. (Third panel) Orbit of the outer planetafter subtracting the orbital fit of the inner planet from the data.Zero phase is defined by the fictitious transit midpoint of the sec-ond planet (denoted as T c , where c is the common subscript for“center”, and “2” refers to the second planet). The error-bars (in-flated by the jitter) are smaller than the size of the points. (Fourthpanel) Orbit of the inner planet after subtracting the orbital fit ofthe outer planet. Zero-phase is defined by the transit midpoint ofHAT-P-13b (denoted as T c ). The error-bars are smaller than thesize of the points. (Fifth panel) Bisector spans (BS) for the Keckspectra including the two template spectra ( § Time from transit center (days) i - band i - band i - band i - band i - band i - band i - band Fig. 3.—
Unbinned instrumental i band transit light curves,acquired with KeplerCam at the FLWO 1.2 m telescope on sevendifferent dates. If the first full transit is assigned N tr = 0 (2008November 8/9 MST, the third event from the top), then the follow-up light curves have N tr = -68, -1, 0, 1, 24, 35, 62 from top tobottom. Superimposed are the best-fit transit model light curves asdescribed in §
3. In the bottom of the figure we show the residualsfrom the fit. Error-bars represent the photon and background shot-noise, plus the readout noise.
TABLE 2Photometry follow-up of HAT-P-13
BJD Mag a σ Mag
Mag(orig) Filter(2,400,000+)54581 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Note . — This table is presented in its entirety in the elec-tronic edition of the Astrophysical Journal. A portion is shownhere for guidance regarding its form and content. a The out-of-transit level has been subtracted. These magni-tudes have been subjected to the EPD and TFA procedures (in“ELTG” mode), carried out simultaneously with the transit fit.
To determine the stellar properties via a set ofisochrones, we used three parameters: the stellar effectivetemperature, the metallicity, and the normalized semi-major axis a/R ⋆ (or related mean stellar density ρ ⋆ ). Wenote that another possible parameter in place of a/R ⋆ would be the stellar surface gravity, but in the case ofplanetary transits the a/R ⋆ quantity typically imposes astronger constraint on the stellar models (Sozzetti et al.2007). (The validity of this assumption, namely that theAT-P-13b,c 5adequate physical model describing our data is a plane-tary transit, as opposed to e.g. a blend, is shown later in § T eff ⋆ , [Fe / H], and log g ⋆ as derived from the SME anal-ysis, we performed a global modeling of the data ( § a/R ⋆ . For T eff ⋆ and [Fe / H]we adopted normal distributions in theMonte Carlo analysis, with dispersions equal to the 1- σ uncertainties from the initial SME analysis.For each combination within the large ( ∼ a/R ⋆ , T eff ⋆ , [Fe / H] values, we searched the stellarisochrones of the Yi et al. (2001) models for the best fitstellar model parameters ( M ⋆ , R ⋆ , age, log g ⋆ , etc.). Wederived the mean values and uncertainties of the physicalparameters based on their a posteriori distribution (seee.g. P´al et al. 2008).We then repeated the SME analysis by fixing the stellarsurface gravity to the refined value of log g ⋆ = 4 . ± . T eff ⋆ ,[Fe / H] and v sin i . The SME analysis was performed onthe first, weaker template observation, and also on thesecond and third, higher S/N pair of template observa-tions taken by Keck/HIRES. While the pair of high S/Ntemplates were acquired very close in time, and the re-spective SME values were consistent to within a smallfraction of the formal error-bars, they were also consis-tent to within 1- σ with the values based on the weakertemplate that was taken much earlier. This consistencyreassured that our stellar atmospheric parameter deter-mination is robust, and the error-bars are realistic. Be-cause the second two templates were of better quality,we adopted the SME values found from these spectrawith simple averaging, yielding T eff ⋆ = 5653 ±
90 K,[Fe / H] = +0 . ± .
08 and v sin i = 2 . ± . − .We adopted these as the final atmospheric parametersfor the star. We then also repeated the isochrone searchfor stellar parameters, obtaining M ⋆ =1 . +0 . − . M ⊙ , R ⋆ =1 . ± . R ⊙ and L ⋆ =2 . ± . L ⊙ . These aresummarized in Tab. 3, along with other stellar proper-ties. Model isochrones from Yi et al. (2001) for metal-licity [Fe / H]=+0 .
41 are plotted in Fig. 4, with the finalchoice of effective temperature T eff ⋆ and a/R ⋆ marked,and encircled by the 1- σ and 2- σ confidence ellipsoids.The second SME iteration at fixed stellar surface gravity(as determined from a/R ⋆ ) changed the metallicity andstellar temperature in such a way that the new ( T eff ⋆ , a/R ⋆ ) values now fall in a more complex region of theisochrones, as compared to the initial SME values, allow-ing for multiple solutions (see Fig. 4, where the originalSME values are marked with a triangle). The distribu-tion of stellar age becomes bimodal with the dominantpeak in the histogram (not shown) being at 5.0 Gyr, anda smaller peak (by a factor of 5) at 7.3 Gyr. This corre-sponds to a slightly bimodal mass distribution with thedominant peak at ∼ . M ⊙ , and much smaller peakaround ∼ . M ⊙ . The asymmetric error-bars givenin Tab. 3 for the mass and age account for the double-peaked distribution.The stellar evolution modeling also yields the abso- The reason for the intersecting isochrones is the “kink” on theevolutionary tracks for M ⋆ & . M ⊙ stars evolving off the mainsequence. a / R * Effective temperature [K]
Fig. 4.—
Stellar isochrones from Yi et al. (2001) for metallicity[Fe / H]=+0 .
41 and ages between 3.6 and 8.0 Gyr in steps of 0.4 Gyr.The final choice of T eff ⋆ and a/R ⋆ are marked and encircled bythe 1- σ and 2- σ confidence ellipsoids. The open triangle denotesthe T eff ⋆ , a/R ⋆ point found during the first SME iteration, beforerefining the stellar surface gravity. lute magnitudes and colors in various photometric pass-bands. We used the apparent magnitudes from the2MASS catalogue (Skrutskie et al. 2006) to determinethe distance of the system, after conversion to the ESOsystem of the models. The reported magnitudes for thisstar are J = 9 . ± . H = 9 . ± . K = 8 . ± . J ESO = 9 . ± . H ESO = 9 . ± .
021 and K ESO =9 . ± . J − K ) = 0 . ± . J − K ) iso = 0 . ± . K apparent magnitude andthe M K = 2 . ± .
12 absolute magnitude derived fromthe above-mentioned modeling to determine the distance:214 ±
12 pc, assuming no reddening. The K band waschosen because it is the longest wavelength band-passwith the smallest expected discrepancies due to molecu-lar lines in the spectrum of the star. Excluding blend scenarios
Following Torres et al. (2007), we explored the pos-sibility that the measured radial velocities are not dueto the (multiple) planet-induced orbital motion of thestar, but are instead caused by distortions in the spec-tral line profiles. This could be due to contaminationfrom a nearby unresolved eclipsing binary, in this casepresumably with a second companion producing the RVsignal corresponding to HAT-P-13c. A bisector analysisbased on the Keck spectra was done as described in § Global modeling of the data
Bakos et al.
TABLE 3Stellar parameters for HAT-P-13
Parameter Value Source T eff ⋆ (K). . . . . . . 5653 ±
90 SME a [Fe / H] (dex) . . . +0 . ± .
08 SME v sin i (km s − ) 2 . ± . v mac (km s − ) . 3 .
83 SME v mic (km s − ). . 0 .
85 SME γ RV (km s − ) . . +14 . ± .
68 DS γ i ) . . . . . . . . . . . 0 . γ i ) . . . . . . . . . . . 0 . M ⋆ ( M ⊙ ) . . . . . . 1 . +0 . − . YY+ a/R ⋆ +SME c R ⋆ ( R ⊙ ) . . . . . . . 1 . ± .
08 YY+ a/R ⋆ +SMElog g ⋆ (cgs) . . . . 4 . ± .
04 YY+ a/R ⋆ +SME L ⋆ ( L ⊙ ) . . . . . . . 2 . ± .
31 YY+ a/R ⋆ +SME V (mag) . . . . . . . 10.622 TASS d B − V (mag) . . 0 . ± .
03 YY+ a/R ⋆ +SME M V (mag) . . . . . 3 . ± .
17 YY+ a/R ⋆ +SME K (mag,ESO) 8 . ± .
017 2MASS e M K (mag,ESO) 2 . ± .
12 YY+ a/R ⋆ +SMEAge (Gyr) . . . . . 5 . +2 . − . YY+ a/R ⋆ +SMEDistance (pc) . . 214 ±
12 YY+ a/R ⋆ +SME a SME = “Spectroscopy Made Easy” package foranalysis of high-resolution spectra Valenti & Piskunov(1996). These parameters depend primarily on SME,with a small dependence on the iterative analysis incor-porating the isochrone search and global modeling ofthe data, as described in the text. b SME+Claret = Based on SME analysis and tablesfrom Claret (2004). c YY+ a/R ⋆ +SME = Yi et al. (2001) isochrones, a/R ⋆ luminosity indicator, and SME results. d Based on the TASS catalogue (Droege et al. 2006). e Based on the transformations from Carpenter (2001).
This section presents a simultaneous fitting of theHATNet photometry, the follow-up light curves, and theRV observations, which we refer to as “global” modeling.It incorporates not only a physical model of the system,but also a description of systematic (instrumental) vari-ations.Our model for the follow-up light curves used analyticformulae based on Mandel & Agol (2002) for the eclipseof a star by a planet, where the stellar flux is describedby quadratic limb-darkening. The limb darkening coef-ficients were taken from the tables by Claret (2004) forthe i band, corresponding to the stellar properties de-termined from the SME analysis ( § p ≡ R p /R ⋆ , the square of the impact parameter b , and the reciprocal of the half duration of the tran-sit ζ/R ⋆ . We chose these parameters because of theirsimple geometric meanings and the fact that they shownegligible correlations (see e.g. Carter et al. 2008).Our model for the HATNet data was the simplified“P1P3” version of the Mandel & Agol (2002) analyticfunctions (an expansion by Legendre polynomials), forthe reasons described in Bakos et al. (2009). The depthof the HATNet transits was adjusted independently inthe fit (the depth was B inst “blending factor” times thatof the follow-up data) to allow for possible contaminationby nearby stars in the under-sampled images of HATNet.As indicated earlier, initial modeling of the RV ob-servations showed deviations from a Keplerian fit highlysuggestive of a second body in the system with a muchlonger period than the transiting planet. Thus, in our global modeling, the RV curve was parametrized by thecombination of an eccentric Keplerian orbit for the innerplanet with semi-amplitude K , and Lagrangian orbitalelements ( k, h ) = e × (cos ω, sin ω ), plus an eccentric Ke-plerian orbit for the outer object with K , k and h ,and a systemic RV zero-point γ . Throughout this paperthe subscripts “1” and “2” will refer to HAT-P-13b andHAT-P-13c, respectively. If the subscript is omitted, werefer to HAT-P-13b.In the past, for single transiting planet scenarios wehave assumed strict periodicity in the individual tran-sit times (Hartman et al. 2009, and references therein),even when a drift in the RV measurements indicated anouter companion (HAT-P-11b, Bakos et al. 2009). Sincethe expected transit timing variations (TTV) for theseplanets were negligible compared to the error-bars of thetransit center measurements, the strict periodicity wasa reasonable hypothesis. Those data were characterizedby two transit centers ( T A and T B ), and all interme-diate transits were interpolated using these two epochsand their corresponding N tr transit numbers. The modelfor the RV data component also implicitly contained itsephemeris information through T A and T B , and thus wascoupled with the photometry data in the time domain.For HAT-P-13b, however, the disturbing force by theouter planet HAT-P-13c is expected to be not insignif-icant, because the RV semi-amplitude of the host star( ∼ . − ) indicates that HAT-P-13c is massive, andthe relatively short period and eccentric orbit (see later)indicate that it moves in relatively close to HAT-P-13b.Thus, the assumption of strict periodicity for HAT-P-13b is not precisely correct. While we performed manyvariations of the global modeling, in our finally adoptedphysical model we assume strict periodicity only for theHATNet data, where the timing error on individual tran-sits can be ∼ N tr = 0 to the first complete high quality follow-up light curve gathered on 2008 November 8/9 (MST).The HATNet data-set was characterized by T c, − and T c, − , covering all transit events observed by HATNet(here the c subscript denotes “center” for the transits ofHAT-P-13b). The transit follow-up observations werecharacterized by their respective times of transit cen-ter: T c, − , T c, − , T c, , T c, , T c, , T c, , T c, . Initialestimates of the T c,i values yielded an initial epoch T c and period P by linear fitting weighted by the respec-tive error-bars of the transit centers. The model for theRV data component of the inner planet contained theephemeris information through the T c and P variables,i.e. it was coupled with the transit data. The global mod-eling was done in an iterative way. After an initial fit tothe transit centers (and other parameters; see later), T c and P were refined, and the fit was repeated.The time dependence of the RV of the outer planetwas described via its hypothetical transit time T c (asif its orbit were edge on), and its period P . The timeof periastron passage T peri can be equivalently used inplace of time of conjunction T c .Altogether, the 21 parameters describing the physi-cal model were T c, − , T c, − (HATNet), T c, − , T c, − ,AT-P-13b,c 7 T c, , T c, , T c, , T c, , T c, (FLWO 1.2 m), R p /R ⋆ , b , ζ/R ⋆ , K , γ , k = e cos ω , h = e sin ω , K , k , h , P and T c . Two additional auxiliary parameters were the in-strumental blend factor B inst of HATNet, and the HAT-Net out-of-transit magnitude, M , HATNet . -0.01-0.005 0 0.005 0.01 0.015 0.02 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 R e l a t i v e i - m ag Time from transit center (days)
Fig. 5.—
Part of the global modeling described in § We extended our physical model with an instrumen-tal model that describes the systematic (non-physical)variations (such as trends) of the follow-up data. Thiswas done in a similar fashion to the analysis presentedin Bakos et al. (2009). The HATNet photometry hadalready been EPD- and TFA-corrected before the globalmodeling, so we only modeled systematic effects in thefollow-up light curves. We chose the “ELTG” method,i.e. EPD was performed in “local” mode with EPD co-efficients defined for each night, whereas TFA was per-formed in “global mode”, with the same coefficients de-scribing the optimal weights for the selected templatestars in the field. The five EPD parameters includethe hour angle (characterizing a monotonic trend thatchanges linearly over a night), the square of the hourangle, and three stellar profile parameters (equivalent toFWHM, elongation, and position angle). The exact func-tional form of the above parameters contained 6 coeffi-cients, including the auxiliary out-of-transit magnitudeof the individual events. The EPD parameters were in-dependent for all 7 nights, implying 42 additional coef-ficients in the global fit. For the global TFA analysiswe chose 19 template stars that had good quality mea-surements for all nights and on all frames, implying 19additional parameters in the fit. Thus, the total numberof fitted parameters was 21 (physical model) + 2 (auxil-iary) + 42 (local EPD) + 19 (TFA) = 84.The joint fit was performed as described in Bakos etal. (2009), by minimizing χ in the parameter space us-ing a hybrid algorithm, combining the downhill simplexmethod (AMOEBA, see Press et al. 1992) with the clas-sical linear least-squares algorithm. We used the par-tial derivatives of the model functions as derived by P´al(2008). Uncertainties in the parameters were derived us-ing the Markov Chain Monte-Carlo method (MCMC, seeFord 2006). Since the eccentricity of the inner systemappeared as marginally significant ( k = − . ± . h = 0 . ± . e = 0 . ± . k and h tozero. The best fit results for the relevant physical param-eters are summarized in Tab. 4 and Tab. 5. The individ-ual transit centers for the FLWO 1.2 m data are given in § T c , − = 2453700 . ± . T c , − =2453878 . ± . B inst = 0 . ± .
06, and M , HATNet = 9 . ± . I band). The fit to theHATNet photometry data was shown earlier in Fig. 1,the orbital fit to the RV data is shown in Fig. 2, andthe fit to the FLWO 1.2 m data is displayed in Fig. 3.The low rms of 3.4 m s − around the orbital fit is due inpart to the use of an iodine free Keck/HIRES templatethat was acquired with higher spectral resolution andhigher S/N than usual. Note that the low rms impliesthe absence of any additional interior planets other thanHAT-P-13b, consistent with expectations given that themassive and eccentric outer planet HAT-P-13c dynam-ically forbids such planets (see e.g. Wittenmyer et al.2007). The stellar jitter required to reconcile the reduced χ with 1 . . − . The low jitter isalso consistent with HAT-P-13 being a chromosphericallyquiet star (based on log R ′ HK and S ).The planetary parameters and their uncertainties canbe derived by the direct combination of the a posteri-ori distributions of the light curve, radial velocity andstellar parameters. We find that the mass of the innerplanet is M p = 0 . +0 . − . M J , the radius is R p = 1 . ± . R J , and its density is ρ p = 0 . +0 . − . g cm − . Thefinal planetary parameters are summarized at the bottomof Table 4. The simultaneous EPD- and TFA-correctedlight curve of all photometry follow-up events is shownin Fig. 5.The outer planet, HAT-P-13c, appears to be very mas-sive, with m sin i = 15 . ± . M J , and orbits the star ina highly eccentric orbit with e = 0 . ± . ω = 176 . ± . ◦ ) is such that ourline of sight is almost along the minor axis, coincidentallyas in the HAT-P-2b system (Bakos et al. 2007). We notethat the periastron passage of HAT-P-13c has not beenwell monitored, and thus the RV fit of the orbit suffersfrom a strong correlation between the quantities K , e and γ , leading to correlated m sin i and e (Fig. 6). DISCUSSION
We present the discovery of HAT-P-13, the first de-tected multi-planet system with a transiting planet. Theinner transiting planet, HAT-P-13b, is an inflated “hotJupiter” in a nearly circular orbit. The outer planet,HAT-P-13c, is both extremely massive and highly eccen-tric, and orbits in a
P > / H] = +0 . ± .
08, the host star is alsoremarkable. As we describe below, this extraordinarysystem is a rich dynamical laboratory, the first to havean accurate clock (HAT-P-13b) and a known perturbingforce (HAT-P-13c). The inner planet will help refine theorbital configuration (and thus true mass) of the outerplanet through transit timing variations (TTVs). Con-versely, the outer planet, through its known perturba- Bakos et al. e m sin(i) Best fit
Fig. 6.—
The orbital parameters, and as a consequence, themass ( m sin i ) and eccentricity e of the outer planet are stronglycorrelated (see § m sin i is due to the similardistribution of the stellar mass (see § tion, may constrain structural parameters and the tidaldissipation rate of the inner planet (Batygin, Boden-heimer, & Laughlin 2009), in addition to all the informa-tion that can be gleaned from the transits of HAT-P-13b. The planet HAT-P-13b
The only other known planet with a host-star asmetal rich as HAT-P-13 ([Fe / H]=+0 . ± .
08) is XO-2b ([
F e/H ] = +0 .
45, Burke et al. 2007). XO-2b, how-ever, has much smaller mass (0.56 M J ) and smaller radius(0.98 R J ) than HAT-P-13b (0 . M J , 1 . R J ). Otherplanets similar to HAT-P-13b include HAT-P-9b ( M p =0 . M J , R p = 1 . R J , Shporer et al. 2008), and XO-1b( M p = 0 . M J , R p = 1 . R J , McCullough et al. 2006).When compared to theoretical models, HAT-P-13b isa clearly inflated planet. The Baraffe et al. (2008) mod-els with solar insolation (at 0.045 AU) are consistentwith the observed mass and radius of HAT-P-13b ei-ther with overall metal content Z=0.02 and a very youngage of 0.05–0.1 Gyr, or with metal content Z=0.1 andage 0.01–0.05 Gyr. Given the fact, however, that thehost star is metal rich, and fairly old (5 . +2 . − . Gyr), it isunlikely that HAT-P-13b is newly formed (50 Myr) andvery metal poor. Comparison with Fortney et al. (2007)leads to similar conclusions. Using these models, HAT-P-13b is broadly consistent only with a 300 Myr planetat 0.02 AU solar distance and core-mass up to 25 M ⊕ ,or a 1 Gyr core-less (pure H/He) planet. If HAT-P-13bhas a significant rocky core, consistent with expectationgiven the high metallicity of HAT-P-13, it must somehowbe inflated beyond models calculated for such a planetwith age and insolation suggested by the data. Nu-merous explanations have been brought up in the pastto explain the inflated radii of certain extrasolar plan-ets (Miller, Fortney, & Jackson 2009; Fabrycky, John-son, & Goodman 2007, and references therein). Oneamong these has been the tidal heating due to non-zeroeccentricity (Bodenheimer, Lin, & Mardling 2001). Noperturbing companion has been found for the inflatedplanets HD 209458b and HAT-P-1b (for an overview,see Mardling 2007). The HAT-P-13 system may be the The equivalent semi-major axis, at which HAT-P-13b wouldreceive the same amount of insolation when orbiting our Sun, is a equiv = 0 . ± . TABLE 4Orbital and planetary parameters for HAT-P-13b
Parameter ValueLight curve parameters, HAT-P-13b P (days) . . . . . . . . . . . . . . . . . . . . 2 . ± . T c (BJD) . . . . . . . . . . . . . . . . . . . 2454779 . ± . T (days) a . . . . . . . . . . . . . . . . . 0 . ± . T = T (days) a . . . . . . . . . . 0 . ± . a/R ⋆ . . . . . . . . . . . . . . . . . . . . . . . . 5 . ± . ζ/R ⋆ . . . . . . . . . . . . . . . . . . . . . . . . 17 . ± . R p /R ⋆ . . . . . . . . . . . . . . . . . . . . . . 0 . ± . b . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . +0 . − . b ≡ a cos i/R ⋆ . . . . . . . . . . . . . . . 0 . +0 . − . i (deg) . . . . . . . . . . . . . . . . . . . . . . 83 . ± . T peri (days) . . . . . . . . . . . . . . . . . 2454780 . ± . K (m s − ) . . . . . . . . . . . . . . . . . . 106 . ± . k . . . . . . . . . . . . . . . . . . . . . . . . . . . . − . ± . h . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . e . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . ω . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 ± ◦ Other RV parameters γ rel (m s − ) . . . . . . . . . . . . . . . . . − . ± . − ) . . . . . . . . . . . . . . . 3 . T s (BJD) . . . . . . . . . . . . . . . . . . . 2454781 . ± . T s, . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . T s, . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . M p ( M J ) . . . . . . . . . . . . . . . . . . . . 0 . +0 . − . R p ( R J ) . . . . . . . . . . . . . . . . . . . . . 1 . ± . C ( M p , R p ) b . . . . . . . . . . . . . . . . 0 . ρ p (g cm − ) . . . . . . . . . . . . . . . . . 0 . +0 . − . a (AU) . . . . . . . . . . . . . . . . . . . . . . 0 . +0 . − . log g p (cgs) . . . . . . . . . . . . . . . . . . 3 . ± . T eq (K) . . . . . . . . . . . . . . . . . . . . . 1653 ± c . . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . F per (erg s − cm − ) d . . . . . . . (1 . ± . · F ap (erg s − cm − ) e . . . . . . . . (1 . ± . · h F i (erg s − cm − ) . . . . . . . . . . (1 . ± . · T : total transit duration, time between first to last contact; T = T : ingress/egress time, time between first and second,or third and fourth contact. b Correlation coefficient between the planetary mass M p andradius R p . c The Safronov number is given by Θ = ( V esc /V orb ) =( a/R p )( M p /M ⋆ ) (see Hansen & Barman 2007). d Incoming flux per unit surface area. first, where the inflated radius can be explained by thenon-zero eccentricity of HAT-P-13b, excited by the outerHAT-P-13c planet orbiting on a highly eccentric orbit.We note that while the eccentricity of HAT-P-13b is non-zero only at the 2- σ level (because k is non-zero at 2- σ ),its pericenter is aligned (to within 4 ◦ ± ◦ ) with thatof the outer planet HAT-P-13c. Additional RV measure-ments and/or space-based timing of a secondary eclipseare necessary for a more accurate determination of the or-bit. Nevertheless, if the apsidal lines of HAT-P-13b andHAT-P-13c are indeed aligned, and, if the configurationis co-planar, then the system is in a tidal fix-point config-uration, as recently noted by Batygin, Bodenheimer, &Laughlin (2009). This configuration imposes constraintson the structure of the inner planet HAT-P-13b. In par-AT-P-13b,c 9 TABLE 5Orbital and planetary parameters for HAT-P-13c
Parameter ValueRV parameters, as induced by HAT-P-13c P (days) . . . . . . . . . . . . . . . . . . . 428 . ± . T c a (BJD) . . . . . . . . . . . . . . . . . 2454870 . ± . K (m s − ) . . . . . . . . . . . . . . . . . 502 ± k . . . . . . . . . . . . . . . . . . . . . . . . . . − . ± . h . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . e . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . ± . ω . . . . . . . . . . . . . . . . . . . . . . . . . . 176 . ± . ◦ T ,peri (days) . . . . . . . . . . . . . . . 2454890 . ± . b T , (days) . . . . . . . . . . . . . . . . 0 . ± . T , = T (days) . . . . . . . . . . 0 . ± . a T s (BJD) . . . . . . . . . . . . . . . . . . 2454912 . ± . T s, (days) . . . . . . . . . . . . . . . . 0 . ± . T s, (days) . . . . . . . . . . . . . . . . 0 . ± . m sin i ( M J ) . . . . . . . . . . . . . . 15 . ± . a (AU) . . . . . . . . . . . . . . . . . . . . . 1 . +0 . − . T , eq (K) . . . . . . . . . . . . . . . . . . . . 340 ± F ,per (erg s − cm − ) d . . . . . (2 . ± . · F ,ap (erg s − cm − ) d . . . . . . (7 . ± . · h F i (erg s − cm − ) d . . . . . . . (3 . ± . · T c would be the center of transit of HAT-P-13c, if its(unknown) inclination was 90 ◦ . b Transits of HAT-P-13c have not been observed. Thevalues are for guidance only, and assume zero impact pa-rameter. c T : total transit duration, time between first to lastcontact, assuming zero impact parameter. T = T :ingress/egress time, time between first and second, or thirdand fourth contact. Note that these values are fictitious,and transits of HAT-P-13c have not been observed. d Incoming flux per unit surface area in periastron, apas-tron, and averaged over the orbit. ticular, Batygin, Bodenheimer, & Laughlin (2009) givelimits on the tidal Love number k b , core-mass and tidalenergy dissipation rate Q b of HAT-P-13b. Transit timing variations
It has long been known that multiple planets in thesame planetary system perturb each other, and this maylead to detectable variations in the transit time of thetransiting planet(s). Transit timing variations have beenanalytically described by Holman & Murray (2005) andAgol et al. (2005), and have been suggested as one of themost effective ways of detecting small mass perturbers oftransiting planets. This has motivated extensive follow-up of known TEPs e.g. by the Transit Light Curve (TLC)project (Holman et al. 2006; Winn et al. 2007), lookingfor companions of TEPs. However, for HAT-P-13b thereis observational (spectroscopic) evidence for an outercomponent (HAT-P-13c), and thus transit timing vari-ations of HAT-P-13b must be present. Transit timingvariations can be used to constrain the mass and orbitalelements of the perturbing planet, in our case those ofHAT-P-13c.The global modeling described in § T c,i ) and the calculated value based on a fixedepoch and period as given in Tab. 4, i.e. in the O − C sense. The resulting TTVs are shown in Fig. 7. Webelieve the error-bars to be realistic as they are the re-sult of a full MCMC analysis, where all parameters arevaried (including the EPD and TFA parameters). Asfurther support for this, the transits around N tr = 0( T c, − , T c, , T c, ) show a standard deviation that is com-parable to the error-bars (and we can safely assume zeroTTV within a ± . N tr = 0 and 24). TTVsof the order of 100 seconds are seen from the best fit pe-riod. Given the large error-bars on our transit centers (ofthe order of 100 seconds), in part due to possible remain-ing systematics in the partial transit light curve events,we consider these departures suggestive only.Nevertheless, it is tempting to compare our results withanalytic formulae presented by Agol et al. (2005) andBorkovits et al. (2003). The HAT-P-13(b,c) system cor-responds to case “ii” of Agol et al. (2005), i.e. with anexterior planet on an eccentric orbit having a much largersemi-major axis than the inner planet. Formulae for thegeneral (non co-planar) case are given by Borkovits etal. (2003, Eq. 46). The TTV effect will depend on thefollowing known parameters for the HAT-P-13 system: M ⋆ , M p , P , i p , P , e , ω , m sin i and T ,peri (theseare listed in Tab. 4 and Tab. 5). The TTV will alsodepend on the following unknown parameters: the trueinclination of HAT-P-13c, i , and the relative angle ofthe orbital normals projected onto the plane of the sky, D = Ω − Ω (see Fig. 1 of Borkovits et al. 2003). In ad-dition to the gravitational perturbation of HAT-P-13c onHAT-P-13b, the barycenter of the inner subsystem (com-posed of the host star HAT-P-13, and the inner planetHAT-P-13b) orbits about the three-body barycenter dueto the massive HAT-P-13c companion. This leads tolight-travel time effects in the transit times of HAT-P-13b (TTVl) that are of the same order of magnitude asthe TTV effect due to the perturbation (TTVg).We have evaluated the analytic formulae includingboth the TTVg and TTVl effects for cases with i =83 . ± . ◦ (corresponding to co-planar inner and outerorbits, and to M = 15 . M J ), and i = 8 . ◦ (an adhoc value yielding an almost face-on orbit with M =105 M J ), and D = 0 ◦ or D = 90 ◦ . These four analyticmodels are illustrated in Fig. 7. The bottom panel ofFig. 7 shows the TTVl and TTVg effects separately forthe i = 83 . ± . ◦ and D = 0 ◦ case. The i = 83 . ± . ◦ ( M = 15 . M J ) cases give TTV variations of the order of15 seconds. The functional form of the i = 83 . ± . ◦ and D = 0 ◦ case appears to follow the observationalvalues, albeit with much smaller amplitude. Curiously,for this case the TTVl effect cancels the TTVg effectto some extent (bottom panel of Fig. 7). Increasing themass of the outer companion by decreasing i at constant m sin i does increase the TTV amplitude up to 100 sec-onds at i = 8 . ◦ , but the functional form changes andno longer resembles the trend seen in the observationaldata.In conclusion, while it is premature to fit the cur-rent data-set with analytic models because of the smallnumber of data-points (7) and the large error-bars ( ∼ TABLE 6Transit timing variations ofHAT-P-13. N tr T c σ T c O-C(BJD) (sec) (sec) −
68 2454581 . . − . − . . − .
50 2454779 . . − .
21 2454782 . . − .
624 2454849 . . .
735 2454882 . . .
262 2454960 . . .
100 sec), these data are not inconsistent with the pres-ence of TTVs, and will prove useful in later analyses. Fu-ture measurements of full transit events with high accu-racy in principle can determine both i and D , i.e., HAT-P-13c may become an RV-based detection with knownorbital orientation and a true mass (even if it does nottransit). -200-100 0 100 200-80 -60 -40 -20 0 20 40 60 80 TT V [ s e c ] N tr transit number i2=i1, Ω - Ω =0i2=i1, Ω - Ω =90i2=8.1, Ω - Ω =0i2=8.1, Ω - Ω =90-40-20 0 20 40-80 -60 -40 -20 0 20 40 60 80 TT V [ s e c ] N tr transit number TTVg + TTVlTTVgTTVl
Fig. 7.— (Top) Transit times of the individual transit eventsof HAT-P-13b. The filled circles correspond to the global analysisdescribed in § i = 83 . ◦ = i , ii) i = 8 . ◦ (at fixed m sin i ,corresponding to a large mass for HAT-P-13c), iii) i = 83 . ◦ andthe mutual inclination of the orbital normals in the plane of thesky is D = Ω − Ω = 90 ◦ (see Fig. 1 of Borkovits et al. (2003)), oriv) i = 8 . ◦ and D = 90 ◦ . (Bottom) The selected case “i” fromabove with zoomed-in vertical scale, showing the TTVg (dashedline) and TTVl (dotted line) effects separately, and their net effect(solid line). The planet HAT-P-13c
The probability of HAT-P-13c transiting the host star,as seen from the Earth, depends on R ⋆ , R p , e , ω ,and a (Kane & von Braun 2009). We evaluated thetransit probability in a Monte Carlo fashion, as part ofthe global modeling, resulting in P ,tr = 0 . ± . R ,p = 1 R J ( P ,tr = 0 . ± . R ,p → ∼
428 days, about 4 times longer than the currentrecord holder HD 80606 (Naef et al. 2001; Fossey, Wald-mann, & Kipping 2009; Moutou et al. 2009). With atrue mass of 15 . M J , we have good reasons to believethat its radius would be around 1 R J , based on heavy-mass transiting planets like HAT-P-2b (8.84 M J , Bakoset al. 2007), XO-3b (11.79 M J , Johns-Krull et al. 2008),Corot-3b (21.66 M J , Deleuil et al. 2008), all having radiiaround 1 R J . If the transits are full (i.e. not grazing),then the transit depth would be similar to that of theinner planet HAT-P-13b. The duration of the transitcould be up to 14 . e and ω .As regards to the nature of HAT-P-13c, even at itsminimal mass (15 . ± . M J ) it is the 10th most massiveplanet out of 327 planets listed on the Extrasolar PlanetEncyclopedia as of 2009 July. Curiously, the recentlyannounced Doppler-detection HD 126614b (Howard etal. 2009) appears to have similar orbital characteristicsto HAT-P-13c. Both are Jovian planets in P > e = 0 . / H] ≈ .
5) stars.As described earlier in §4.2, we have good hopes thatin the near future precise TTV measurements of the in-ner planet HAT-P-13b, will constrain the orbital incli-nation and thus the true mass of HAT-P-13c. Further,such TTV variations can also constrain the mutual in-clination of HAT-P-13b and HAT-P-13c. Measuring thesky-projected angle between the stellar spin axis and theorbital normal of HAT-P-13b (the inner planet) via theRossiter-McLaughlin effect will shed light on the migra-tion history of HAT-P-13b, and by inference the scatter-ing history between HAT-P-13b and HAT-P-13c.HATNet operations have been funded by NASA grantsNNG04GN74G, NNX08AF23G and SAO IR&D grants.Work of G. ´A.B. and J. Johnson were supported by thePostdoctoral Fellowship of the NSF Astronomy and As-trophysics Program (AST-0702843 and AST-0702821, re-spectively). We acknowledge partial support also fromthe Kepler Mission under NASA Cooperative Agree-ment NCC2-1390 (D.W.L., PI). G.K. thanks the Hun-AT-P-13b,c 11garian Scientific Research Foundation (OTKA) for sup-port through grant K-60750. G.T. acknowledges partialsupport from NASA Origins grant NNX09AF59G. Thisresearch has made use of Keck telescope time grantedthrough NOAO (program A146Hr,A264Hr) and NASA (N128Hr,N145Hr). We are grateful to Josh Winn andMatthew Holman for their flexibility in swapping nightsat the FLWO 1.2 m telescope. We thank the anonymousreferee for the useful comments that improved this paper.