HAT-P-13b,c: a transiting hot Jupiter with a massive outer companion on an eccentric orbit1


We report on the discovery of a planetary system with a close-in transiting hot Jupiter on a near circular 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  days, with transit epoch (BJD) and duration  d. The outer planet, HAT-P-13c orbits the star with  days and nominal transit center (assuming zero impact parameter) of (BJD) or time of periastron passage (BJD). Transits of the outer planet have not been observed, and may not be present. The host star has a mass of , radius of , effective temperature  K, and is rather metal rich with . The inner planetary companion has a mass of , and radius of  yielding a mean density of . The outer companion has , and orbits on a highly eccentric orbit of . While we have not detected significant transit timing variations of HAT-P-13b, due to gravitational and light-travel time effects, future observations will constrain the orbital inclination of HAT-P-13c, along with its mutual inclination to HAT-P-13b. The HAT-P-13 (b,c) double-planet system may prove extremely valuable for theoretical studies of the formation and dynamics of planetary systems.

Subject headings:
planetary systems — stars: individual (HAT-P-13, GSC 3416-00543) techniques: spectroscopic, photometric

1. Introduction

Radial velocity (RV) surveys have shown that multiple-planet stellar systems are common. For example, Wright et al. (2007) concluded that the occurrence of additional planets among stars already having one known planet must be greater than 30%. Thus, one would expect that out of the published transiting extrasolar planet (TEP) systems24, there should be a number of systems with additional companions; these companions should make their presence known through radial velocity variations of the parent stars, even if they do not themselves transit. The fact that so far no TEPs in multiple planet systems have been reported is somewhat surprising based on the statistics from RV surveys (see also Fabrycky, 2009). Recently Smith et al. (2009) searched the light curves of 24 transiting planets for outer companions, finding no evidence for a double transiting system.

The Hungarian-made Automated Telescope Network (HATNet) survey (Bakos et al., 2002, 2004) has been a major contributor to the discovery of TEPs. Operational since 2003, it has covered approximately 10% of the Northern sky, searching for TEPs around bright stars ( mag). HATNet operates six wide field instruments: four at the Fred Lawrence Whipple Observatory (FLWO) in Arizona, and two on the roof of the Submillimeter Array hangar (SMA) of SAO at Hawaii. Since 2006, HATNet has announced and published 12 TEPs. A study similar to that of Smith et al. (2009) has been carried out on 9 known transiting planets from the HATNet project by Fabrycky (private communication) with a similar null result.

In this work we report on the 13th discovery of HATNet, the detection of the first known system with a transiting inner planet (HAT-P-13b) which also contains a second, outer planet (HAT-P-13c), as detected by the radial velocity variation of the host star. There have been examples of transiting systems where the RV variations do show a long term trend, such as HAT-P-7b (Winn et al., 2009) and HAT-P-11b (Bakos et al., 2009), but no orbital solution has yet been presented for any of these outer companions, simply because there has not been a long enough timespan to cover the period, or at least observe the long term trend changing sign. While no transits of HAT-P-13c have yet been detected, the probability that the companion actually transits the star is non-negligible if the orbits of the two planets are nearly co-planar (the pure geometric transit probability for HAT-P-13c is 1.3%, see § 4.3). The system is particularly interesting because the outer planet has both a high eccentricity and a very high mass. These properties, in turn, should induce transit timing variations (TTVs) of the inner planet of the order of 10 seconds (standard deviation, Agol et al., 2005). Such TTVs may be used to constrain the orbital parameters of the outer planet, including the inclination with respect to the line of sight and with respect to the orbital plane of the inner planet.

In § 2 of this paper we summarize the observations, including the photometric detection, and follow-up observations. § 3 describes the analysis of the data, such as the stellar parameter determination (§ 3.1), blend modeling (§ 3.2), and global modeling of the data (§ 3.3). We discuss our findings in § 4.

2. Observations

2.1. Photometric detection

The transits of HAT-P-13b were detected with the HAT-5 telescope. The region around GSC 3416-00543, a field internally labeled as 136, was observed on a nightly basis between 2005 November 25 and 2006 May 20, whenever weather conditions permitted. We gathered 4021 exposures of 5 minutes at a 5.5-minute cadence. Each image contained approximately 20000 stars down to . For the brightest stars in the field we achieved a per-image photometric precision of 3.1 mmag.

Figure 1.— The unbinned light curve of HAT-P-13 including all 4021 instrumental  band measurements obtained with the HAT-5 telescope of HATNet (see text for details), and folded with the period of  days (which is the result of the fit described in § 3).

The calibration of the HATNet frames was done utilizing standard procedures. The calibrated frames were then subjected to star detection and astrometry, as described in Pál & Bakos (2006). Aperture photometry was performed on each image at the stellar centroids derived from the 2MASS catalog (Skrutskie et al., 2006) and the individual astrometrical solutions. Description of numerous details related to the reduction are also given in Pál (2009). The resulting light curves were decorrelated against trends using the External Parameter Decorrelation technique in “constant” mode (EPD, see Bakos et al., 2009) and the Trend Filtering Algorithm (TFA, see Kovács et al., 2005). The light curves were searched for periodic box-like signals using the Box Least Squares method (BLS, see Kovács et al., 2002). We detected a significant signal in the light curve of GSC 3416-00543 (also known as 2MASS 08393180+4721073; , ; J2000), with a depth of  mmag, and a period of  days. The dip significance parameter (Kovács et al., 2002) of the signal was DSP = 17. The dip had a relative duration (first to last contact) of , corresponding to a total duration of  hours (see Fig. 1).

2.2. Reconnaissance Spectroscopy

One of the important tools used in our survey for establishing whether the transit-feature in the light curve of a candidate is due to astrophysical phenomena other than a planet transiting a star is the CfA Digital Speedometer (DS; Latham, 1992), mounted on the FLWO 1.5 m telescope. This yields high-resolution spectra with low signal-to-noise ratio sufficient to derive radial velocities with moderate precision (roughly 1 ), and to determine the effective temperature and surface gravity of the host star. With this facility we are able to weed out certain false alarms, such as eclipsing binaries and multiple stellar systems.

We obtained  spectra spanning  days with the DS. The RV measurements of HAT-P-13 showed an rms residual of , consistent with no detectable RV variation. The spectra were single-lined, showing no evidence for more than one star in the system. Atmospheric parameters for the star, including the initial estimates of effective temperature , surface gravity (log cgs), and projected rotational velocity , were derived as described by Torres et al. (2002). The effective temperature and surface gravity correspond to a G4 dwarf (Cox, 2000). The mean line-of-sight velocity of HAT-P-13 is .

2.3. High resolution, high S/N spectroscopy and the search for radial velocity signal components

Given the significant transit detection by HATNet, and the encouraging DS results, we proceeded with the follow-up of this candidate by obtaining high-resolution and high S/N spectra to characterize the radial velocity variations and to determine the stellar parameters with higher precision. Using the HIRES instrument (Vogt et al., 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 had low signal-to-noise ratio, thus we repeated the template observations during a later run, and acquired two high quality templates. We used the last two templates in the analysis. The observations were made between 2008 March 22 and 2009 June 5. The relative radial velocity (RV) measurements are listed in Tab. 1, and shown in Fig. 2.

Figure 2.— (Top) Radial-velocity measurements from Keck for HAT-P-13, along with the best 2-planet orbital fit, shown as a function of BJD (see § 3). The center-of-mass velocity has been subtracted. (Second panel) Phased residuals after subtracting the orbital fit (also see § 3). The rms variation of the residuals is about 3.4 , requiring a jitter of  to be added in quadrature to the individual errors to yield a reduced of 1.0. The error-bars in this panel have been inflated accordingly. Note the different vertical scale of the panels. (Third panel) Orbit of the outer planet after subtracting the orbital fit of the inner planet from the data. Zero phase is defined by the fictitious transit midpoint of the second planet (denoted as , where is the common subscript for “center”, and “2” refers to the second planet). The error-bars (inflated by the jitter) are smaller than the size of the points. (Fourth panel) Orbit of the inner planet after subtracting the orbital fit of the outer planet. Zero-phase is defined by the transit midpoint of HAT-P-13b (denoted as ). The error-bars are smaller than the size of the points. (Fifth panel) Bisector spans (BS) for the Keck spectra including the two template spectra (§ 3.2). The mean value has been subtracted. (Bottom) Relative S activity index values for the Keck spectra (see Hartman et al., 2009). Open circles in the phase-plots indicate data that appears twice due to plotting outside the [0,1] phase domain.
BJD RV25 26 BS S27
    (2,454,000) () () () ()
Table 1 Relative radial velocity, bisector and activity index measurements of HAT-P-13.

Observations and reductions have been carried out in an identical way to that described in earlier HATNet discovery papers, such as Bakos et al. (2009). References for the Keck iodine cell observations and the reduction of the radial velocities are given in Marcy & Butler (1992) and Butler et al. (1996).

Initial fits of these RVs to a single-planet Keplerian orbit were quite satisfactory but soon revealed a slight residual trend that became more significant with time. This is the reason for the observations extending over more than one year, as opposed to just a few months as necessary to confirm simpler purely sinusoidal variations seen 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 likely due to a more distant body in the same system, possibly a massive planet. Preliminary two-planet orbital solutions provided a much improved fit (although with a weakly determined outer period due to the short duration of the observations), and more importantly, held up after additional observations. With the data so far gathered and presented in this work, a false alarm probability (FAP) analysis for the Keplerian orbit of the second planet yielded an FAP of 0.00001.

A more thorough analysis using a two-component least squares Fourier fit (see Barning, 1963, for the single-component case), with one component fixed to the known frequency of the short-period inner planet, also confirmed the existence of the long period component, and indicated that the latter is due to a highly non-sinusoidal motion. A number of other low frequency peaks were eliminated as being aliases yielding worse quality fits. A three-harmonic representation of the fit yielded an RMS of , fairly close to the value expected from the formal errors. Full modeling of this complex motion of two Keplerian orbits superposed, simultaneously with the photometry, is described in our global fit of § 3.3.

2.4. Photometric follow-up observations

Figure 3.— Unbinned instrumental  band transit light curves, acquired with KeplerCam at the FLWO 1.2 m telescope on seven different dates. If the first full transit is assigned (2008 November 8/9 MST, the third event from the top), then the follow-up light curves have -68, -1, 0, 1, 24, 35, 62 from top to bottom. Superimposed are the best-fit transit model light curves as described in § 3. In the bottom of the figure we show the residuals from the fit. Error-bars represent the photon and background shot-noise, plus the readout noise.

We observed 7 transit events of HAT-P-13 with the KeplerCam CCD on the FLWO 1.2 m telescope between 2008 April 24 and 2009 May 8. All observations were carried out in the  band, and the typical exposure time was 15 seconds. The reduction of the images, including calibration, astrometry, and photometry, was performed as described in Bakos et al. (2009).

We performed EPD and TFA against trends simultaneously with the light curve modeling (for more details, see § 3). From the series of apertures, for each night, we chose the one yielding the smallest residual rms for producing the final light curve. These are shown in the upper plot of Fig. 3, superimposed with our best fit transit light curve models (see also § 3).

BJD Mag28 Mag(orig) Filter

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.

Table 2Photometry follow-up of HAT-P-13

3. Analysis

3.1. Properties of the parent star

We derived the initial stellar atmospheric parameters by using the first template spectrum obtained with the Keck/HIRES instrument. We used the SME package of Valenti & Piskunov (1996) along with the atomic-line database of Valenti & Fischer (2005), which yielded the following initial values and uncertainties (which we have conservatively increased, to include our estimates of the systematic errors): effective temperature  K, stellar surface gravity  (cgs), metallicity  dex, and projected rotational velocity .

Further analysis of the spectra has shown that the host star is chromospherically quiet with and (on an absolute scale).

To determine the stellar properties via a set of isochrones, we used three parameters: the stellar effective temperature, the metallicity, and the normalized semi-major axis  (or related mean stellar density ). We note that another possible parameter in place of  would be the stellar surface gravity, but in the case of planetary transits the  quantity typically imposes a stronger constraint on the stellar models (Sozzetti et al., 2007). (The validity of this assumption, namely that the adequate physical model describing our data is a planetary transit, as opposed to e.g. a blend, is shown later in § 3.2.) With quadratic limb darkening coefficients (listed in Tab. 3) from Claret (2004) appropriate for the values of , , and as derived from the SME analysis, we performed a global modeling of the data (§ 3.3), yielding a full Monte-Carlo distribution of . For  and we adopted normal distributions in the Monte Carlo analysis, with dispersions equal to the 1- uncertainties from the initial SME analysis.

For each combination within the large () set of , ,  values, we searched the stellar isochrones of the Yi et al. (2001) models for the best fit stellar model parameters (, , age, , etc.). We derived the mean values and uncertainties of the physical parameters based on their a posteriori distribution (see e.g. Pál et al., 2008).

We then repeated the SME analysis by fixing the stellar surface gravity to the refined value of based on the isochrone search, and only adjusting ,  and . The SME analysis was performed on the first, weaker template observation, and also on the second and third, higher S/N pair of template observations taken by Keck/HIRES. While the pair of high S/N templates were acquired very close in time, and the respective SME values were consistent to within a small fraction of the formal error-bars, they were also consistent to within 1- with the values based on the weaker template that was taken much earlier. This consistency reassured that our stellar atmospheric parameter determination is robust, and the error-bars are realistic. Because the second two templates were of better quality, we adopted the SME values found from these spectra with simple averaging, yielding  K, and . We adopted these as the final atmospheric parameters for the star. We then also repeated the isochrone search for stellar parameters, obtaining =, = and =. These are summarized in Tab. 3, along with other stellar properties. Model isochrones from Yi et al. (2001) for metallicity = are plotted in Fig. 4, with the final choice of effective temperature and  marked, and encircled by the 1- and 2- confidence ellipsoids. The second SME iteration at fixed stellar surface gravity (as determined from ) changed the metallicity and stellar temperature in such a way that the new (, ) values now fall in a more complex region of the isochrones, as compared to the initial SME values, allowing for multiple solutions (see Fig. 4, where the original SME values are marked with a triangle).29 The distribution of stellar age becomes bimodal with the dominant peak in the histogram (not shown) being at 5.0 Gyr, and a smaller peak (by a factor of 5) at 7.3 Gyr. This corresponds to a slightly bimodal mass distribution with the dominant peak at , and much smaller peak around . The asymmetric error-bars given in Tab. 3 for the mass and age account for the double-peaked distribution.

Figure 4.— Stellar isochrones from Yi et al. (2001) for metallicity = and ages between 3.6 and 8.0 Gyr in steps of 0.4 Gyr. The final choice of and  are marked and encircled by the 1- and 2- confidence ellipsoids. The open triangle denotes the , point found during the first SME iteration, before refining the stellar surface gravity.

The stellar evolution modeling also yields the absolute magnitudes and colors in various photometric passbands. We used the apparent magnitudes from the 2MASS catalogue (Skrutskie et al., 2006) to determine the distance of the system, after conversion to the ESO system of the models. The reported magnitudes for this star are , and , which we transformed to , and , following Carpenter (see 2001). These yield a color of , fully consistent with the expected, isochrone-based . We thus relied on the 2MASS apparent magnitude and the absolute magnitude derived from the above-mentioned modeling to determine the distance:  pc, assuming no reddening. The  band was chosen because it is the longest wavelength band-pass with the smallest expected discrepancies due to molecular lines in the spectrum of the star.

Parameter Value Source
(K) SME30
(dex) SME
() SME
() SME
() SME
() DS
() YY++SME31
() YY++SME
(cgs) YY++SME
() YY++SME
(mag) 10.622 TASS32
(mag) YY++SME
(mag) YY++SME
(mag,ESO) 2MASS33
(mag,ESO) YY++SME
Age (Gyr) YY++SME
Distance (pc) YY++SME
Table 3 Stellar parameters for HAT-P-13

3.2. Excluding blend scenarios

Following Torres et al. (2007), we explored the possibility that the measured radial velocities are not due to the (multiple) planet-induced orbital motion of the star, but are instead caused by distortions in the spectral line profiles. This could be due to contamination from a nearby unresolved eclipsing binary, in this case presumably with a second companion producing the RV signal corresponding to HAT-P-13c. A bisector analysis based on the Keck spectra was done as described in §3 of Hartman et al. (2009).

We detect no significant variation in the bisector spans (see Fig. 2, fifth panel). The correlation between the radial velocity and the bisector variations is also insignificant. Therefore, we conclude that the velocity variations of the host star are real, and can be interpreted as being due to a close-in planet, with the added complication from an outer object that we account for in the modeling that follows. Because of the negligible bisector variations that show no correlation with the radial velocities, we found no need to perform detailed blend modeling of hierarchical triple (quadruple) scenarios, such as that done for the case of HAT-P-12b (Hartman et al., 2009).

3.3. Global modeling of the data

This section presents a simultaneous fitting of the HATNet photometry, the follow-up light curves, and the RV 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) variations.

Our model for the follow-up light curves used analytic formulae based on Mandel & Agol (2002) for the eclipse of a star by a planet, where the stellar flux is described by quadratic limb-darkening. The limb darkening coefficients were taken from the tables by Claret (2004) for the  band, corresponding to the stellar properties determined from the SME analysis (§ 3.1). The transit shape was parametrized by the normalized planetary radius , the square of the impact parameter , and the reciprocal of the half duration of the transit . We chose these parameters because of their simple geometric meanings and the fact that they show negligible correlations (see e.g. Carter et al., 2008).

Our model for the HATNet data was the simplified “P1P3” version of the Mandel & Agol (2002) analytic functions (an expansion by Legendre polynomials), for the reasons described in Bakos et al. (2009). The depth of the HATNet transits was adjusted independently in the fit (the depth was “blending factor” times that of the follow-up data) to allow for possible contamination by nearby stars in the under-sampled images of HATNet.

As indicated earlier, initial modeling of the RV observations showed deviations from a Keplerian fit highly suggestive of a second body in the system with a much longer period than the transiting planet. Thus, in our global modeling, the RV curve was parametrized by the combination of an eccentric Keplerian orbit for the inner planet with semi-amplitude , and Lagrangian orbital elements , plus an eccentric Keplerian orbit for the outer object with , and , and a systemic RV zero-point . Throughout this paper the subscripts “1” and “2” will refer to HAT-P-13b and HAT-P-13c, respectively. If the subscript is omitted, we refer to HAT-P-13b.

In the past, for single transiting planet scenarios we have assumed strict periodicity in the individual transit times (Hartman et al., 2009, and references therein), even when a drift in the RV measurements indicated an outer companion (HAT-P-11b, Bakos et al., 2009). Since the expected transit timing variations (TTV) for these planets were negligible compared to the error-bars of the transit center measurements, the strict periodicity was a reasonable hypothesis. Those data were characterized by two transit centers ( and ), and all intermediate transits were interpolated using these two epochs and their corresponding transit numbers. The model for the RV data component also implicitly contained its ephemeris information through and , and thus was coupled with the photometry data in the time domain.

For HAT-P-13b, however, the disturbing force by the outer planet HAT-P-13c is expected to be not insignificant, because the RV semi-amplitude of the host star () indicates that HAT-P-13c is massive, and the 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 many variations of the global modeling, in our finally adopted physical model we assume strict periodicity only for the HATNet data, where the timing error on individual transits can be  seconds. In this final model we allow for departure from such a periodicity for the individual transit times for the seven follow-up photometry observations. In practice, we assigned the transit number to the first complete high quality follow-up light curve gathered on 2008 November 8/9 (MST). The HATNet data-set was characterized by and , covering all transit events observed by HATNet (here the subscript denotes “center” for the transits of HAT-P-13b). The transit follow-up observations were characterized by their respective times of transit center: , , , , , , . Initial estimates of the values yielded an initial epoch and period by linear fitting weighted by the respective error-bars of the transit centers. The model for the RV data component of the inner planet contained the ephemeris information through the and variables, i.e. it was coupled with the transit data. The global modeling was done in an iterative way. After an initial fit to the transit centers (and other parameters; see later), and were refined, and the fit was repeated.

The time dependence of the RV of the outer planet was described via its hypothetical transit time (as if its orbit were edge on), and its period . The time of periastron passage can be equivalently used in place of time of conjunction .

Altogether, the 21 parameters describing the physical model were , (HATNet), , , , , , , (FLWO 1.2 m), , , , , , , , , , , and . Two additional auxiliary parameters were the instrumental blend factor of HATNet, and the HATNet out-of-transit magnitude, .

Figure 5.— Part of the global modeling described in § 3 are corrections for systematic variations of the light curves via simultaneous fitting with the physical model of the transit. The figure shows the resulting EPD- and TFA-corrected light curves for all 7 follow-up events (small gray points), and the merged and binned light curve (bin-size 0.005 days).

We extended our physical model with an instrumental model that describes the systematic (non-physical) variations (such as trends) of the follow-up data. This was done in a similar fashion to the analysis presented in Bakos et al. (2009). The HATNet photometry had already been EPD- and TFA-corrected before the global modeling, so we only modeled systematic effects in the follow-up light curves. We chose the “ELTG” method, i.e. EPD was performed in “local” mode with EPD coefficients defined for each night, whereas TFA was performed in “global mode”, with the same coefficients describing the optimal weights for the selected template stars in the field. The five EPD parameters include the hour angle (characterizing a monotonic trend that changes linearly over a night), the square of the hour angle, and three stellar profile parameters (equivalent to FWHM, elongation, and position angle). The exact functional form of the above parameters contained 6 coefficients, including the auxiliary out-of-transit magnitude of the individual events. The EPD parameters were independent for all 7 nights, implying 42 additional coefficients in the global fit. For the global TFA analysis we chose 19 template stars that had good quality measurements for all nights and on all frames, implying 19 additional parameters in the fit. Thus, the total number of fitted parameters was 21 (physical model) + 2 (auxiliary) + 42 (local EPD) + 19 (TFA) = 84.

The joint fit was performed as described in Bakos et al. (2009), by minimizing  in the parameter space using a hybrid algorithm, combining the downhill simplex method (AMOEBA, see Press et al., 1992) with the classical linear least-squares algorithm. We used the partial derivatives of the model functions as derived by Pál (2008). Uncertainties in the parameters were derived using the Markov Chain Monte-Carlo method (MCMC, see Ford, 2006). Since the eccentricity of the inner system appeared as marginally significant (, , implying ), and also because the physical model dictates that in the presence of a massive outer companion, the inner planet could maintain a non-zero eccentricity, we did not fix and to zero. The best fit results for the relevant physical parameters are summarized in Tab. 4 and Tab. 5. The individual transit centers for the FLWO 1.2 m data are given in § 4.2. Other parameters fitted but not listed in the table were:  (BJD),  (BJD), , and ( band). The fit to the HATNet photometry data was shown earlier in Fig. 1, the orbital fit to the RV data is shown in Fig. 2, and the fit to the FLWO 1.2 m data is displayed in Fig. 3. The low rms of 3.4  around the orbital fit is due in part to the use of an iodine free Keck/HIRES template that was acquired with higher spectral resolution and higher S/N than usual. Note that the low rms implies the absence of any additional interior planets other than HAT-P-13b, consistent with expectations given that the massive and eccentric outer planet HAT-P-13c dynamically forbids such planets (see e.g. Wittenmyer et al., 2007). The stellar jitter required to reconcile the reduced  with for the RV data is . The low jitter is also consistent with HAT-P-13 being a chromospherically quiet star (based on  and ).

The planetary parameters and their uncertainties can be derived by the direct combination of the a posteriori distributions of the light curve, radial velocity and stellar parameters. We find that the mass of the inner planet is , the radius is , and its density is . The final planetary parameters are summarized at the bottom of Table 4. The simultaneous EPD- and TFA-corrected light curve of all photometry follow-up events is shown in Fig. 5.

The outer planet, HAT-P-13c, appears to be very massive, with , and orbits the star in a highly eccentric orbit with . The orientation of the orbit () is such that our line of sight is almost along the minor axis, coincidentally as in the HAT-P-2b system (Bakos et al., 2007). We note that the periastron passage of HAT-P-13c has not been well monitored, and thus the RV fit of the orbit suffers from a strong correlation between the quantities , and , leading to correlated and (Fig. 6).

Figure 6.— The orbital parameters, and as a consequence, the mass () and eccentricity of the outer planet are strongly correlated (see § 3). Gray dots represent the results of the MCMC runs (20,000 trials). The final solution is marked with an open circle. The bimodal distribution in is due to the similar distribution of the stellar mass (see § 3.1).
               Parameter Value
Light curve parameters, HAT-P-13b
    (days) 35
    (days) 36
RV parameters, as induced by HAT-P-13b
Other RV parameters
   Jitter ()
Secondary eclipse parameters for HAT-P-13b
Planetary parameters for HAT-P-13b
    () 39
    () 40
Table 4Orbital and planetary parameters for HAT-P-13b
               Parameter Value
RV parameters, as induced by HAT-P-13c
   41 (BJD)
Fictitious light curve parameters, HAT-P-13c42
   43 (days)
Fictitious secondary eclipse parameters for HAT-P-13c44
Planetary parameters for HAT-P-13c
    () 45
    () 46
    () 47
Table 5Orbital and planetary parameters for HAT-P-13c

4. Discussion

We present the discovery of HAT-P-13, the first detected multi-planet system with a transiting planet. The inner transiting planet, HAT-P-13b, is an inflated “hot Jupiter” in a nearly circular orbit. The outer planet, HAT-P-13c, is both extremely massive and highly eccentric, and orbits in a yr orbit . With an iron abundance of , the host star is also remarkable. As we describe below, this extraordinary system is a rich dynamical laboratory, the first to have an accurate clock (HAT-P-13b) and a known perturbing force (HAT-P-13c). The inner planet will help refine the orbital configuration (and thus true mass) of the outer planet through transit timing variations (TTVs). Conversely, the outer planet, through its known perturbation, may constrain structural parameters and the tidal dissipation rate of the inner planet (Batygin, Bodenheimer, & Laughlin, 2009), in addition to all the information that can be gleaned from the transits of HAT-P-13b.

4.1. The planet HAT-P-13b

The only other known planet with a host-star as metal rich as HAT-P-13 (=) is XO-2b (, Burke et al., 2007). XO-2b, however, has much smaller mass (0.56 ) and smaller radius (0.98 ) than HAT-P-13b (, ). Other planets similar to HAT-P-13b include HAT-P-9b (, , Shporer et al., 2008), and XO-1b (, , McCullough et al., 2006).

When compared to theoretical models, HAT-P-13b is a clearly inflated planet. The Baraffe et al. (2008) models with solar insolation (at 0.045 AU) are consistent with the observed mass and radius of HAT-P-13b either with overall metal content Z=0.02 and a very young age of 0.05–0.1 Gyr, or with metal content Z=0.1 and age 0.01–0.05 Gyr.48 Given the fact, however, that the host star is metal rich, and fairly old ( Gyr), it is unlikely that HAT-P-13b is newly formed (50 Myr) and very 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 planet at 0.02 AU solar distance and core-mass up to 25 , or a 1 Gyr core-less (pure H/He) planet. If HAT-P-13b has a significant rocky core, consistent with expectation given the high metallicity of HAT-P-13, it must somehow be inflated beyond models calculated for such a planet with age and insolation suggested by the data. Numerous explanations have been brought up in the past to explain the inflated radii of certain extrasolar planets (Miller, Fortney, & Jackson, 2009; Fabrycky, Johnson, & Goodman, 2007, and references therein). One among these has been the tidal heating due to non-zero eccentricity (Bodenheimer, Lin, & Mardling, 2001). No perturbing companion has been found for the inflated planets HD 209458b and HAT-P-1b (for an overview, see Mardling, 2007). The HAT-P-13 system may be the first, where the inflated radius can be explained by the non-zero eccentricity of HAT-P-13b, excited by the outer HAT-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 is non-zero at 2-), its pericenter is aligned (to within ) with that of the outer planet HAT-P-13c. Additional RV measurements and/or space-based timing of a secondary eclipse are necessary for a more accurate determination of the orbit. Nevertheless, if the apsidal lines of HAT-P-13b and HAT-P-13c are indeed aligned, and, if the configuration is co-planar, then the system is in a tidal fix-point configuration, as recently noted by Batygin, Bodenheimer, & Laughlin (2009). This configuration imposes constraints on the structure of the inner planet HAT-P-13b. In particular, Batygin, Bodenheimer, & Laughlin (2009) give limits on the tidal Love number , core-mass and tidal energy dissipation rate of HAT-P-13b.

4.2. Transit timing variations

It has long been known that multiple planets in the same planetary system perturb each other, and this may lead to detectable variations in the transit time of the transiting planet(s). Transit timing variations have been analytically described by Holman & Murray (2005) and Agol et al. (2005), and have been suggested as one of the most effective ways of detecting small mass perturbers of transiting 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), looking for companions of TEPs. However, for HAT-P-13b there is observational (spectroscopic) evidence for an outer component (HAT-P-13c), and thus transit timing variations of HAT-P-13b must be present. Transit timing variations can be used to constrain the mass and orbital elements of the perturbing planet, in our case those of HAT-P-13c.

The global modeling described in § 3.3 treats the transit centers of the FLWO 1.2 m telescope follow-up as independent variables, i.e., it automatically provides the basis for a TTV analysis. Throughout this work, by TTV we refer to the time difference between the observed transit center () and the calculated value based on a fixed epoch and period as given in Tab. 4, i.e. in the sense. The resulting TTVs are shown in Fig. 7. We believe the error-bars to be realistic as they are the result of a full MCMC analysis, where all parameters are varied (including the EPD and TFA parameters). As further support for this, the transits around (,,) show a standard deviation that is comparable to the error-bars (and we can safely assume zero TTV within a  day time range). It is also apparent from the plot that the smallest error-bars correspond to the full transit observations ( 0 and 24). TTVs of the order of 100 seconds are seen from the best fit period. Given the large error-bars on our transit centers (of the order of 100 seconds), in part due to possible remaining systematics in the partial transit light curve events, we consider these departures suggestive only.

Nevertheless, it is tempting to compare our results with analytic formulae presented by Agol et al. (2005) and Borkovits et al. (2003). The HAT-P-13(b,c) system corresponds to case “ii” of Agol et al. (2005), i.e. with an exterior planet on an eccentric orbit having a much larger semi-major axis than the inner planet. Formulae for the general (non co-planar) case are given by Borkovits et al. (2003, Eq. 46). The TTV effect will depend on the following known parameters for the HAT-P-13 system: , , , , , , , and (these are listed in Tab. 4 and Tab. 5). The TTV will also depend on the following unknown parameters: the true inclination of HAT-P-13c, , and the relative angle of the orbital normals projected onto the plane of the sky, (see Fig. 1 of Borkovits et al., 2003). In addition to the gravitational perturbation of HAT-P-13c on HAT-P-13b, the barycenter of the inner subsystem (composed of the host star HAT-P-13, and the inner planet HAT-P-13b) orbits about the three-body barycenter due to the massive HAT-P-13c companion. This leads to light-travel time effects in the transit times of HAT-P-13b (TTVl) that are of the same order of magnitude as the TTV effect due to the perturbation (TTVg).

We have evaluated the analytic formulae including both the TTVg and TTVl effects for cases with (corresponding to co-planar inner and outer orbits, and to ), and (an ad hoc value yielding an almost face-on orbit with ), and or . These four analytic models are illustrated in Fig. 7. The bottom panel of Fig. 7 shows the TTVl and TTVg effects separately for the and case. The () cases give TTV variations of the order of 15 seconds. The functional form of the and case appears to follow the observational values, albeit with much smaller amplitude. Curiously, for this case the TTVl effect cancels the TTVg effect to some extent (bottom panel of Fig. 7). Increasing the mass of the outer companion by decreasing at constant does increase the TTV amplitude up to 100 seconds at , but the functional form changes and no longer resembles the trend seen in the observational data.

In conclusion, while it is premature to fit the current data-set with analytic models because of the small number of data-points (7) and the large error-bars ( sec), these data are not inconsistent with the presence of TTVs, and will prove useful in later analyses. Future measurements of full transit events with high accuracy in principle can determine both and , i.e., HAT-P-13c may become an RV-based detection with known orbital orientation and a true mass (even if it does not transit).

Figure 7.— (Top) Transit times of the individual transit events of HAT-P-13b. The filled circles correspond to the global analysis described in § 4.2. The large filled square, the open circle and the open square correspond to the apastron, conjunction and periastron of the perturber HAT-P-13c. Overlaid are analytic models for the transit timing variation due to gravitational effects (TTVg) and light-travel time effects (TTVl) for four scenarios: i) the inclination of HAT-P-13c is , ii) (at fixed , corresponding to a large mass for HAT-P-13c), iii) and the mutual inclination of the orbital normals in the plane of the sky is (see Fig. 1 of Borkovits et al. (2003)), or iv) and . (Bottom) The selected case “i” from above with zoomed-in vertical scale, showing the TTVg (dashed line) and TTVl (dotted line) effects separately, and their net effect (solid line).
(BJD) (sec) (sec)
Table 6 Transit timing variations of HAT-P-13.

4.3. The planet HAT-P-13c

The probability of HAT-P-13c transiting the host star, as seen from the Earth, depends on , , , , and (Kane & von Braun, 2009). We evaluated the transit probability in a Monte Carlo fashion, as part of the global modeling, resulting in if ( if ). This derivation assumes an isotropic inclination distribution. However, dynamical constraints, such as precise measurements of TTV effects, or analysis of orbital stability, may further limit the allowed inclination range, and thus increase or decrease the chance of transits.

Unfortunately, our HATNet and FLWO 1.2 m datasets do not cover any time interval around the expected transit times of HAT-P-13c, thus we can not prove or refute the existence of such transits. Nevertheless, it is an interesting thought experiment to characterize the putative transits of HAT-P-13c, should they occur. HAT-P-13c would be the longest period transiting planet discovered to date, with a semi-major axis of over 1 AU, and a period of  days, about 4 times longer than the current record holder HD 80606 (Naef et al., 2001; Fossey, Waldmann, & Kipping, 2009; Moutou et al., 2009). With a true mass of , we have good reasons to believe that its radius would be around 1, based on heavy-mass transiting planets like HAT-P-2b (8.84 , Bakos et al., 2007), XO-3b (11.79 , Johns-Krull et al., 2008), Corot-3b (21.66 , Deleuil et al., 2008), all having radii around 1 . If the transits are full (i.e. not grazing), then the transit depth would be similar to that of the inner planet HAT-P-13b. The duration of the transit could be up to  hours. Follow-up observations of such transits would require either a world-wide effort, or uninterrupted space-based observations. The next transit center, according to the present analysis, will occur at 2010 April 12 9am UTC. Since we are looking at the orbit of HAT-P-13c along the semi-latus rectum (parallel to the minor-axis), the chance for an occultation of the planet by the star has a very similar chance of occurance as the primary eclipse. Observations of the secondary eclipse could greatly decrease the error on and .

As regards to the nature of HAT-P-13c, even at its minimal mass () it is the 10th most massive planet out of 327 planets listed on the Extrasolar Planet Encyclopedia as of 2009 July. Curiously, the recently announced Doppler-detection HD 126614b (Howard et al., 2009) appears to have similar orbital characteristics to HAT-P-13c. Both are Jovian planets in  yr, orbits around metal-rich () stars. As described earlier in § 4.2, we have good hopes that in the near future precise TTV measurements of the inner planet HAT-P-13b, will constrain the orbital inclination and thus the true mass of HAT-P-13c. Further, such TTV variations can also constrain the mutual inclination of HAT-P-13b and HAT-P-13c. Measuring the sky-projected angle between the stellar spin axis and the orbital normal of HAT-P-13b (the inner planet) via the Rossiter-McLaughlin effect will shed light on the migration history of HAT-P-13b, and by inference the scattering history between HAT-P-13b and HAT-P-13c.

HATNet operations have been funded by NASA grants NNG04GN74G, NNX08AF23G and SAO IR&D grants. Work of G.Á.B. and J. Johnson were supported by the Postdoctoral Fellowship of the NSF Astronomy and Astrophysics Program (AST-0702843 and AST-0702821, respectively). We acknowledge partial support also from the Kepler Mission under NASA Cooperative Agreement NCC2-1390 (D.W.L., PI). G.K. thanks the Hungarian Scientific Research Foundation (OTKA) for support through grant K-60750. G.T. acknowledges partial support from NASA Origins grant NNX09AF59G. This research has made use of Keck telescope time granted through NOAO (program A146Hr,A264Hr) and NASA (N128Hr,N145Hr). We are grateful to Josh Winn and Matthew Holman for their flexibility in swapping nights at the FLWO 1.2 m telescope. We thank the anonymous referee for the useful comments that improved this paper.


  1. affiliation: Based in part on observations obtained at the W. M. Keck Observatory, which is operated by the University of California and the California Institute of Technology. Keck time has been granted by NOAO (A146Hr,A264Hr) and NASA (N128Hr,N145Hr).
  2. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  3. affiliation: NSF Fellow
  4. affiliation: Department of Astronomy, University of California, Berkeley, CA
  5. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  6. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  7. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  8. affiliation: Konkoly Observatory, Budapest, Hungary
  9. affiliation: Department of Physics and Astronomy, San Francisco State University, San Francisco, CA
  10. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  11. affiliation: Institute for Astronomy, University of Hawaii, Honolulu, HI 96822; NSF Postdoctoral Fellow
  12. affiliation: Department of Astronomy, University of California, Berkeley, CA
  13. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  14. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  15. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  16. affiliation: Department of Astronomy, Eötvös Loránd University, Budapest, Hungary.
  17. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  18. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  19. affiliation: Konkoly Observatory, Budapest, Hungary
  20. affiliation: Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
  21. affiliation: Hungarian Astronomical Association, Budapest, Hungary
  22. affiliation: Hungarian Astronomical Association, Budapest, Hungary
  23. affiliation: Hungarian Astronomical Association, Budapest, Hungary
  24. www.exoplanet.eu
  25. The fitted zero-point that is on an arbitrary scale (denoted as in Tab. 4) has been subtracted from the velocities.
  26. The values for do not include the jitter.
  27. S values are on a relative scale.
  28. The out-of-transit level has been subtracted. These magnitudes have been subjected to the EPD and TFA procedures (in “ELTG” mode), carried out simultaneously with the transit fit.
  29. The reason for the intersecting isochrones is the “kink” on the evolutionary tracks for stars evolving off the main sequence.
  30. SME = “Spectroscopy Made Easy” package for analysis of high-resolution spectra Valenti & Piskunov (1996). These parameters depend primarily on SME, with a small dependence on the iterative analysis incorporating the isochrone search and global modeling of the data, as described in the text.
  31. YY++SME = Yi et al. (2001) isochrones,  luminosity indicator, and SME results.
  32. Based on the TASS catalogue (Droege et al., 2006).
  33. Based on the transformations from Carpenter (2001).
  34. footnotetext: SME+Claret = Based on SME analysis and tables from Claret (2004).
  35. : total transit duration, time between first to last contact; : ingress/egress time, time between first and second, or third and fourth contact.
  36. : total transit duration, time between first to last contact; : ingress/egress time, time between first and second, or third and fourth contact.
  37. Correlation coefficient between the planetary mass  and radius .
  38. The Safronov number is given by (see Hansen & Barman, 2007).
  39. Incoming flux per unit surface area.
  40. footnotemark:
  41. would be the center of transit of HAT-P-13c, if its (unknown) inclination was .
  42. Transits of HAT-P-13c have not been observed. The values are for guidance only, and assume zero impact parameter.
  43. : total transit duration, time between first to last contact, assuming zero impact parameter. : ingress/egress time, time between first and second, or third and fourth contact. Note that these values are fictitious, and transits of HAT-P-13c have not been observed.
  44. would be the center of transit of HAT-P-13c, if its (unknown) inclination was .
  45. Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
  46. Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
  47. Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
  48. The equivalent semi-major axis, at which HAT-P-13b would receive the same amount of insolation when orbiting our Sun, is  AU.


  1. Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567
  2. Cox, A. N. 2000, Allen’s Astrophysical Quantities
  3. Bakos, G. Á., Lázár, J., Papp, I., Sári, P. & Green, E. M. 2002, PASP, 114, 974
  4. Bakos, G. Á., Noyes, R. W., Kovács, G., Stanek, K. Z., Sasselov, D. D., & Domsa, I. 2004, PASP, 116, 266
  5. Bakos, G. Á., et al. 2007, ApJ, 670, 826
  6. Bakos, G. Á., et al. 2009, arXiv:0901.0282
  7. Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
  8. Barning, F. J. M. 1963, Bull. Astron. Inst. Netherlands, 17, 22
  9. Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, arXiv:0907.5019
  10. Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466
  11. Borkovits, T., Érdi, B., Forgács-Dajka, E., & Kovács, T. 2003, A&A, 398, 1091
  12. Burke, C. J., et al. 2007, ApJ, 671, 2115
  13. Butler, R. P. et al. 1996, PASP, 108, 500
  14. Carpenter, J. M. 2001, AJ, 121, 2851
  15. Carter, J. A., Yee, J. C., Eastman, J., Gaudi, B. S., & Winn, J. N. 2008, ApJ, 689, 499
  16. Claret, A. 2004, A&A, 428, 1001
  17. Deleuil, M., Deeg, H. J., Alonso, R., Bouchy, F., & Rouan, D. 2008, arXiv:0810.0919
  18. Deeming, T. J. 1975, Ap&SS, 36, 137
  19. Droege, T. F., Richmond,M. W., Sallman, M. P., & Creager, R. P. 2006, PASP, 118, 1666
  20. Fabrycky, D. C., Johnson, E. T., & Goodman, J. 2007, ApJ, 665, 754
  21. Fabrycky, D. C. 2009, IAU Symposium, 253, 173
  22. Ford, E. 2006, ApJ, 642, 505
  23. Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
  24. Fossey, S. J., Waldmann, I. P., & Kipping, D. M. 2009, MNRAS, 396, L16
  25. Hansen, B. M. S., & Barman, T. 2007, ApJ, 671, 861
  26. Hartman, J. D., et al. 2009, arXiv:0904.4704
  27. Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288
  28. Holman, M. J., et al. 2006, ApJ, 652, 1715
  29. Howard, W. A., et al. 2009, submitted to ApJ
  30. Kane, S. R., & von Braun, K. 2009, IAU Symposium, 253, 358
  31. Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  32. Kovács, G., Bakos, G. Á., & Noyes, R. W. 2005, MNRAS, 356, 557
  33. Johns-Krull, C. M., et al. 2008, ApJ, 677, 657
  34. Latham, D. W. 1992, in IAU Coll. 135, Complementary Approaches to Double and Multiple Star Research, ASP Conf. Ser. 32, eds. H. A. McAlister & W. I. Hartkopf (San Francisco: ASP), 110
  35. McCullough, P. R.,et al. 2006, ApJ, 648, 1228
  36. Pál, A., & Bakos, G. Á. 2006, PASP, 118, 1474
  37. Pál, A., Bakos, G. Á., Noyes, R. W. & Torres, G. 2008, proceedings of IAU Symp. 253 “Transiting Planets“, ed. by F. Pont (arXiv:0807.1530)
  38. Pál, A. 2008b, MNRAS, 390, 281
  39. Pál, A., PhD thesis, (arXiv:0906.3486)
  40. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P., 1992, Numerical Recipes in C: the art of scientific computing, Second Edition, Cambridge University Press
  41. Mandel, K., & Agol, E. 2002, ApJ, 580, L171
  42. Mardling, R. A. 2007, MNRAS, 382, 1768
  43. Marcy, G. W., & Butler, R. P. 1992, PASP, 104, 270
  44. McLaughlin, D. B. 1924, ApJ, 60, 22
  45. Miller, N., Fortney, J. J., & Jackson, B. 2009, arXiv:0907.1268
  46. Moutou, C., et al. 2009, A&A, 498, L5
  47. Naef, D., et al. 2001, A&A, 375, L27
  48. Shporer, A., et al. 2008, arXiv:0806.4008
  49. Skrutskie, M. F., et al. 2006, AJ, 131, 1163
  50. Smith, A. M. S., et al. 2009, arXiv:0906.3414
  51. Sozzetti, A. et al. 2007, ApJ, 664, 1190
  52. Torres, G., Boden, A. F., Latham, D. W., Pan, M. & Stefanik, R. P. 2002, AJ, 124, 1716
  53. Torres, G. et al. 2007, ApJ, 666, 121
  54. Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141
  55. Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595
  56. Vogt, S. S. et al. 1994, Proc. SPIE, 2198, 362
  57. Yi, S. K. et al. 2001, ApJS, 136, 417
  58. Winn, J. N., et al. 2007, AJ, 133, 1828
  59. Winn, J. N., Johnson, J. A., Albrecht, S., Howard, A. W., Marcy, G. W., Crossfield, I. J., & Holman, M. J. 2009, ApJ, 703, L99
  60. Wittenmyer, R. A., Endl, M., Cochran, W. D., & Levison, H. F. 2007, AJ, 134, 1276
  61. Wright, J. T., et al. 2007, ApJ, 657, 533
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.