HATP13b,c: a transiting hot Jupiter with a massive
outer companion on an eccentric orbit^{1}
Abstract
We report on the discovery of a planetary system with a closein transiting hot Jupiter on a near circular orbit and a massive outer planet on a highly eccentric orbit. The inner planet, HATP13b, transits the bright V=10.622 G4 dwarf star GSC 341600543 every days, with transit epoch (BJD) and duration d. The outer planet, HATP13c 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 HATP13b, due to gravitational and lighttravel time effects, future observations will constrain the orbital inclination of HATP13c, along with its mutual inclination to HATP13b. The HATP13 (b,c) doubleplanet system may prove extremely valuable for theoretical studies of the formation and dynamics of planetary systems.
Subject headings:
planetary systems — stars: individual (HATP13, GSC 341600543) techniques: spectroscopic, photometric1. Introduction
Radial velocity (RV) surveys have shown that multipleplanet 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)
systems
The Hungarianmade 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 (HATP13b) which also contains a second, outer planet (HATP13c), 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 HATP7b (Winn et al., 2009) and HATP11b (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 HATP13c have yet been detected, the probability that the companion actually transits the star is nonnegligible if the orbits of the two planets are nearly coplanar (the pure geometric transit probability for HATP13c 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 followup 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 HATP13b were detected with the HAT5 telescope. The region around GSC 341600543, 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.5minute cadence. Each image contained approximately 20000 stars down to . For the brightest stars in the field we achieved a perimage photometric precision of 3.1 mmag.
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 boxlike 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 341600543 (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 transitfeature 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 highresolution spectra with low signaltonoise 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 HATP13 showed an rms residual of , consistent with no detectable RV variation. The spectra were singlelined, 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 lineofsight velocity of HATP13 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 followup of this candidate by obtaining highresolution 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 iodinefree templates. The first template had low signaltonoise 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.
BJD  RV 

BS  S 


(2,454,000)  ()  ()  ()  ()  
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 singleplanet 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 twoplanet 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 twocomponent least squares Fourier fit (see Barning, 1963, for the singlecomponent case), with one component fixed to the known frequency of the shortperiod inner planet, also confirmed the existence of the long period component, and indicated that the latter is due to a highly nonsinusoidal motion. A number of other low frequency peaks were eliminated as being aliases yielding worse quality fits. A threeharmonic 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 followup observations
We observed 7 transit events of HATP13 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  Mag 
Mag(orig)  Filter  

(2,400,000)  
i  
i  
i  
i  
i  
i  
i  
i 
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.
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 atomicline 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 semimajor 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 MonteCarlo 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 errorbars, 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
errorbars 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).
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, isochronebased . We thus relied on the 2MASS apparent magnitude and the absolute magnitude derived from the abovementioned modeling to determine the distance: pc, assuming no reddening. The band was chosen because it is the longest wavelength bandpass with the smallest expected discrepancies due to molecular lines in the spectrum of the star.
Parameter  Value  Source 

(K)  SME 

(dex)  SME  
()  SME  
()  SME  
()  SME  
()  DS  
SME+Claret  
SME+Claret  
()  YY++SME 

()  YY++SME  
(cgs)  YY++SME  
()  YY++SME  
(mag)  10.622  TASS 
(mag)  YY++SME  
(mag)  YY++SME  
(mag,ESO)  2MASS 

(mag,ESO)  YY++SME  
Age (Gyr)  YY++SME  
Distance (pc)  YY++SME 
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) planetinduced 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 HATP13c. 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 closein 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 HATP12b (Hartman et al., 2009).
3.3. Global modeling of the data
This section presents a simultaneous fitting of the HATNet photometry, the followup 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 followup 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 limbdarkening. 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 followup data) to allow for possible contamination by nearby stars in the undersampled 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 semiamplitude , and Lagrangian orbital elements , plus an eccentric Keplerian orbit for the outer object with , and , and a systemic RV zeropoint . Throughout this paper the subscripts “1” and “2” will refer to HATP13b and HATP13c, respectively. If the subscript is omitted, we refer to HATP13b.
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 (HATP11b, Bakos et al., 2009). Since the expected transit timing variations (TTV) for these planets were negligible compared to the errorbars 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 HATP13b, however, the disturbing force by the outer planet HATP13c is expected to be not insignificant, because the RV semiamplitude of the host star () indicates that HATP13c is massive, and the relatively short period and eccentric orbit (see later) indicate that it moves in relatively close to HATP13b. Thus, the assumption of strict periodicity for HATP13b 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 followup photometry observations. In practice, we assigned the transit number to the first complete high quality followup light curve gathered on 2008 November 8/9 (MST). The HATNet dataset was characterized by and , covering all transit events observed by HATNet (here the subscript denotes “center” for the transits of HATP13b). The transit followup 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 errorbars 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 outoftransit magnitude, .
We extended our physical model with an instrumental model that describes the systematic (nonphysical) variations (such as trends) of the followup 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 TFAcorrected before the global modeling, so we only modeled systematic effects in the followup 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 outoftransit 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 leastsquares 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 MonteCarlo 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 nonzero 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 HATP13b, consistent with expectations given that the massive and eccentric outer planet HATP13c 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 HATP13 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 TFAcorrected light curve of all photometry followup events is shown in Fig. 5.
The outer planet, HATP13c, 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 HATP2b system (Bakos et al., 2007). We note that the periastron passage of HATP13c 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).
Parameter  Value 

Light curve parameters, HATP13b  
(days)  
()  
(days)


(days)


(deg)  
(days)  
RV parameters, as induced by HATP13b  
()  
Other RV parameters  
()  
Jitter ()  
Secondary eclipse parameters for HATP13b  
(BJD)  
Planetary parameters for HATP13b  
()  
()  


()  
(AU)  
(cgs)  
(K)  


() 

() 

() 
Parameter  Value 

RV parameters, as induced by HATP13c  
(days)  


()  
(days)  
Fictitious light curve parameters, HATP13c 



(days)  
Fictitious secondary eclipse parameters for HATP13c 

(BJD)  
(days)  
(days)  
Planetary parameters for HATP13c  
()  
(AU)  
(K)  
() 

() 

() 
4. Discussion
We present the discovery of HATP13, the first detected multiplanet system with a transiting planet. The inner transiting planet, HATP13b, is an inflated “hot Jupiter” in a nearly circular orbit. The outer planet, HATP13c, 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 (HATP13b) and a known perturbing force (HATP13c). 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 HATP13b.
4.1. The planet HATP13b
The only other known planet with a hoststar as metal rich as HATP13 (=) is XO2b (, Burke et al., 2007). XO2b, however, has much smaller mass (0.56 ) and smaller radius (0.98 ) than HATP13b ( , ). Other planets similar to HATP13b include HATP9b (, , Shporer et al., 2008), and XO1b (, , McCullough et al., 2006).
When compared to theoretical models, HATP13b 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
HATP13b 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.
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 followup 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 HATP13b there is observational (spectroscopic) evidence for an outer component (HATP13c), and thus transit timing variations of HATP13b 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 HATP13c.
The global modeling described in § 3.3 treats the transit centers of the FLWO 1.2 m telescope followup 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 errorbars 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 errorbars (and we can safely assume zero TTV within a day time range). It is also apparent from the plot that the smallest errorbars 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 errorbars 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 HATP13(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 semimajor axis than the inner planet. Formulae for the general (non coplanar) case are given by Borkovits et al. (2003, Eq. 46). The TTV effect will depend on the following known parameters for the HATP13 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 HATP13c, , 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 HATP13c on HATP13b, the barycenter of the inner subsystem (composed of the host star HATP13, and the inner planet HATP13b) orbits about the threebody barycenter due to the massive HATP13c companion. This leads to lighttravel time effects in the transit times of HATP13b (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 coplanar inner and outer orbits, and to ), and (an ad hoc value yielding an almost faceon 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 dataset with analytic models because of the small number of datapoints (7) and the large errorbars ( 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., HATP13c may become an RVbased detection with known orbital orientation and a true mass (even if it does not transit).
OC  

(BJD)  (sec)  (sec)  
4.3. The planet HATP13c
The probability of HATP13c 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 HATP13c, 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 HATP13c, should they occur. HATP13c would be the longest period transiting planet discovered to date, with a semimajor 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 heavymass transiting planets like HATP2b (8.84 , Bakos et al., 2007), XO3b (11.79 , JohnsKrull et al., 2008), Corot3b (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 HATP13b. The duration of the transit could be up to hours. Followup observations of such transits would require either a worldwide effort, or uninterrupted spacebased 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 HATP13c along the semilatus rectum (parallel to the minoraxis), 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 HATP13c, 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 Dopplerdetection HD 126614b (Howard et al., 2009) appears to have similar orbital characteristics to HATP13c. Both are Jovian planets in yr, orbits around metalrich () stars. As described earlier in § 4.2, we have good hopes that in the near future precise TTV measurements of the inner planet HATP13b, will constrain the orbital inclination and thus the true mass of HATP13c. Further, such TTV variations can also constrain the mutual inclination of HATP13b and HATP13c. Measuring the skyprojected angle between the stellar spin axis and the orbital normal of HATP13b (the inner planet) via the RossiterMcLaughlin effect will shed light on the migration history of HATP13b, and by inference the scattering history between HATP13b and HATP13c.
Footnotes
 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).
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: NSF Fellow
 affiliation: Department of Astronomy, University of California, Berkeley, CA
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: Konkoly Observatory, Budapest, Hungary
 affiliation: Department of Physics and Astronomy, San Francisco State University, San Francisco, CA
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: Institute for Astronomy, University of Hawaii, Honolulu, HI 96822; NSF Postdoctoral Fellow
 affiliation: Department of Astronomy, University of California, Berkeley, CA
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: Department of Astronomy, Eötvös Loránd University, Budapest, Hungary.
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: Konkoly Observatory, Budapest, Hungary
 affiliation: HarvardSmithsonian Center for Astrophysics, Cambridge, MA, gbakos@cfa.harvard.edu
 affiliation: Hungarian Astronomical Association, Budapest, Hungary
 affiliation: Hungarian Astronomical Association, Budapest, Hungary
 affiliation: Hungarian Astronomical Association, Budapest, Hungary
 www.exoplanet.eu
 The fitted zeropoint that is on an arbitrary scale (denoted as in Tab. 4) has been subtracted from the velocities.
 The values for do not include the jitter.
 S values are on a relative scale.
 The outoftransit 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.
 The reason for the intersecting isochrones is the “kink” on the evolutionary tracks for stars evolving off the main sequence.
 SME = “Spectroscopy Made Easy” package for analysis of highresolution 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.
 YY++SME = Yi et al. (2001) isochrones, luminosity indicator, and SME results.
 Based on the TASS catalogue (Droege et al., 2006).
 Based on the transformations from Carpenter (2001).
 footnotetext: SME+Claret = Based on SME analysis and tables from Claret (2004).
 : total transit duration, time between first to last contact; : ingress/egress time, time between first and second, or third and fourth contact.
 : total transit duration, time between first to last contact; : ingress/egress time, time between first and second, or third and fourth contact.
 Correlation coefficient between the planetary mass and radius .
 The Safronov number is given by (see Hansen & Barman, 2007).
 Incoming flux per unit surface area.
 footnotemark:
 would be the center of transit of HATP13c, if its (unknown) inclination was .
 Transits of HATP13c have not been observed. The values are for guidance only, and assume zero impact parameter.
 : 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 HATP13c have not been observed.
 would be the center of transit of HATP13c, if its (unknown) inclination was .
 Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
 Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
 Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
 The equivalent semimajor axis, at which HATP13b would receive the same amount of insolation when orbiting our Sun, is AU.
References
 Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567
 Cox, A. N. 2000, Allen’s Astrophysical Quantities
 Bakos, G. Á., Lázár, J., Papp, I., Sári, P. & Green, E. M. 2002, PASP, 114, 974
 Bakos, G. Á., Noyes, R. W., Kovács, G., Stanek, K. Z., Sasselov, D. D., & Domsa, I. 2004, PASP, 116, 266
 Bakos, G. Á., et al. 2007, ApJ, 670, 826
 Bakos, G. Á., et al. 2009, arXiv:0901.0282
 Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
 Barning, F. J. M. 1963, Bull. Astron. Inst. Netherlands, 17, 22
 Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, arXiv:0907.5019
 Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466
 Borkovits, T., Érdi, B., ForgácsDajka, E., & Kovács, T. 2003, A&A, 398, 1091
 Burke, C. J., et al. 2007, ApJ, 671, 2115
 Butler, R. P. et al. 1996, PASP, 108, 500
 Carpenter, J. M. 2001, AJ, 121, 2851
 Carter, J. A., Yee, J. C., Eastman, J., Gaudi, B. S., & Winn, J. N. 2008, ApJ, 689, 499
 Claret, A. 2004, A&A, 428, 1001
 Deleuil, M., Deeg, H. J., Alonso, R., Bouchy, F., & Rouan, D. 2008, arXiv:0810.0919
 Deeming, T. J. 1975, Ap&SS, 36, 137
 Droege, T. F., Richmond,M. W., Sallman, M. P., & Creager, R. P. 2006, PASP, 118, 1666
 Fabrycky, D. C., Johnson, E. T., & Goodman, J. 2007, ApJ, 665, 754
 Fabrycky, D. C. 2009, IAU Symposium, 253, 173
 Ford, E. 2006, ApJ, 642, 505
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
 Fossey, S. J., Waldmann, I. P., & Kipping, D. M. 2009, MNRAS, 396, L16
 Hansen, B. M. S., & Barman, T. 2007, ApJ, 671, 861
 Hartman, J. D., et al. 2009, arXiv:0904.4704
 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288
 Holman, M. J., et al. 2006, ApJ, 652, 1715
 Howard, W. A., et al. 2009, submitted to ApJ
 Kane, S. R., & von Braun, K. 2009, IAU Symposium, 253, 358
 Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
 Kovács, G., Bakos, G. Á., & Noyes, R. W. 2005, MNRAS, 356, 557
 JohnsKrull, C. M., et al. 2008, ApJ, 677, 657
 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
 McCullough, P. R.,et al. 2006, ApJ, 648, 1228
 Pál, A., & Bakos, G. Á. 2006, PASP, 118, 1474
 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)
 Pál, A. 2008b, MNRAS, 390, 281
 Pál, A., PhD thesis, (arXiv:0906.3486)
 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
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171
 Mardling, R. A. 2007, MNRAS, 382, 1768
 Marcy, G. W., & Butler, R. P. 1992, PASP, 104, 270
 McLaughlin, D. B. 1924, ApJ, 60, 22
 Miller, N., Fortney, J. J., & Jackson, B. 2009, arXiv:0907.1268
 Moutou, C., et al. 2009, A&A, 498, L5
 Naef, D., et al. 2001, A&A, 375, L27
 Shporer, A., et al. 2008, arXiv:0806.4008
 Skrutskie, M. F., et al. 2006, AJ, 131, 1163
 Smith, A. M. S., et al. 2009, arXiv:0906.3414
 Sozzetti, A. et al. 2007, ApJ, 664, 1190
 Torres, G., Boden, A. F., Latham, D. W., Pan, M. & Stefanik, R. P. 2002, AJ, 124, 1716
 Torres, G. et al. 2007, ApJ, 666, 121
 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141
 Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595
 Vogt, S. S. et al. 1994, Proc. SPIE, 2198, 362
 Yi, S. K. et al. 2001, ApJS, 136, 417
 Winn, J. N., et al. 2007, AJ, 133, 1828
 Winn, J. N., Johnson, J. A., Albrecht, S., Howard, A. W., Marcy, G. W., Crossfield, I. J., & Holman, M. J. 2009, ApJ, 703, L99
 Wittenmyer, R. A., Endl, M., Cochran, W. D., & Levison, H. F. 2007, AJ, 134, 1276
 Wright, J. T., et al. 2007, ApJ, 657, 533