ApJ, accepted
Preprint typeset using LATEX style emulateapj v. 05/04/06
HAT-P-13b,c: A TRANSITING HOT JUPITER WITH A MASSIVE OUTER COMPANION ON AN ECCENTRIC
ORBIT†
G. Á. Bakos1,2 , A. W. Howard3 , R. W. Noyes1 , J. Hartman1 , G. Torres1 , Géza Kovács4 , D. A. Fischer5 ,
D. W. Latham1 , J. A. Johnson6 , G. W. Marcy3 , D. D. Sasselov1 , R. P. Stefanik1 , B. Sipőcz1,7 , Gábor Kovács1 ,
G. A. Esquerdo1 , A. Pál4,1 , J. Lázár8 , I. Papp8 , P. Sári8
arXiv:0907.3525v3 [astro-ph.EP] 12 Oct 2009
ApJ, accepted
ABSTRACT
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 P = 2.916260±0.000010 days, with
transit epoch Tc = 2454779.92979 ± 0.00038 (BJD) and duration 0.1345 ± 0.0017 d. The outer planet,
HAT-P-13c orbits the star with P2 = 428.5±3.0 days and nominal transit center (assuming zero impact
parameter) of T2c = 2454870.4 ± 1.8 (BJD) or time of periastron passage T2,peri = 2454890.05 ± 0.48
(BJD). Transits of the outer planet have not been observed, and may not be present. The host star
has a mass of 1.22+0.05
−0.10 M⊙ , radius of 1.56 ± 0.08 R⊙ , effective temperature 5653 ± 90 K, and is rather
metal rich with [Fe/H] = +0.41 ± 0.08. The inner planetary companion has a mass of 0.853+0.029
−0.046 MJ ,
−3
and radius of 1.281±0.079 RJ yielding a mean density of 0.498+0.103
g
cm
.
The
outer
companion
has
−0.069
m2 sin i2 = 15.2 ± 1.0 MJ , and orbits on a highly eccentric orbit of e2 = 0.691 ± 0.018. 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 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 ∼ 50 published transiting extrasolar planet
(TEP) systems10 , 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 compan1 Harvard-Smithsonian Center for Astrophysics, Cambridge,
MA, gbakos@cfa.harvard.edu
2 NSF Fellow
3 Department of Astronomy, University of California, Berkeley,
CA
4 Konkoly Observatory, Budapest, Hungary
5 Department of Physics and Astronomy, San Francisco State
University, San Francisco, CA
6 Institute for Astronomy, University of Hawaii, Honolulu, HI
96822; NSF Postdoctoral Fellow
7 Department of Astronomy, Eötvös Loránd University, Budapest, Hungary.
8 Hungarian Astronomical Association, Budapest, Hungary
† 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).
10 www.exoplanet.eu
ions, 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
(8 . I . 12.5 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 par-
2
Bakos et al.
ticularly 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
I ∼ 14.0. For the brightest stars in the field we achieved
a per-image photometric precision of 3.1 mmag.
2.3. High resolution, high S/N spectroscopy and the
9.77
I-band magnitude
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 km s−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 7 spectra spanning 37 days with the DS.
The RV measurements of HAT-P-13 showed an rms residual of 0.68 km s−1 , 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 Teff⋆ = 5500 ± 100 K, surface gravity log g⋆ = 4.0 ± 0.25 (log cgs), and projected rotational
velocity v sin i = 0.5 ± 1.0 km s−1 , 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
+14.69 ± 0.68 km s−1 .
9.775
9.78
9.785
9.79
9.795
9.8
9.805
9.81
9.815
-0.4
-0.2
0
0.2
0.4
Orbital phase
Fig. 1.— The unbinned light curve of HAT-P-13 including all
4021 instrumental I band measurements obtained with the HAT-5
telescope of HATNet (see text for details), and folded with the period of P = 2.9162595 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 341600543 (also known as 2MASS 08393180+4721073; α =
08h 39m31.82s, δ = +47◦21′ 07.4′′ ; J2000), with a depth
of ∼ 6.2 mmag, and a period of P = 2.9163 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 q ≈ 0.0461 ± 0.0006, corresponding to a
total duration of P q ≈ 3.228 ± 0.040 hours (see Fig. 1).
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.
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
HAT-P-13b,c
3
TABLE 1
Relative radial velocity, bisector and activity index measurements of
HAT-P-13.
BJD
(2,454,000+)
RVa
(m s−1 )
σRV b
(m s−1 )
BS
(m s−1 )
σBS
(m s−1 )
547.89534 . . . . . . .
548.80173 . . . . . . . .
602.73395 . . . . . . . .
602.84690 . . . . . . . .
603.73414 . . . . . . . .
603.84323 . . . . . . . .
633.77240 . . . . . . . .
634.75907 . . . . . . . .
635.75475 . . . . . . . .
636.74968 . . . . . . . .
727.13851 . . . . . . . .
728.13190 . . . . . . . .
778.07302 . . . . . . . .
779.08375 . . . . . . . .
780.09369 . . . . . . . .
791.11130 . . . . . . . .
809.99576 . . . . . . . .
810.92159 . . . . . . . .
839.06085 . . . . . . . .
865.02660 . . . . . . . .
867.90311 . . . . . . . .
928.83635 . . . . . . . .
929.84447 . . . . . . . .
955.86964 . . . . . . . .
956.86327 . . . . . . . .
963.85163 . . . . . . . .
983.74976 . . . . . . . .
984.76460 . . . . . . . .
985.73856 . . . . . . . .
985.74584 . . . . . . . .
985.75333 . . . . . . . .
986.76358 . . . . . . . .
988.74066 . . . . . . . .
21.42
···
17.80
22.56
182.48
202.33
211.86
40.20
183.01
204.23
215.53
38.36
42.19
219.58
87.11
193.44
−14.92
157.29
−124.84
−350.51
−387.07
−188.44
−189.69
−90.43
95.55
−20.24
139.67
−35.66
120.74
···
···
125.20
149.15
4.30
···
1.43
1.61
1.35
2.02
1.93
1.95
2.19
1.71
1.81
1.64
1.45
1.72
1.59
1.66
2.33
3.43
1.49
1.41
2.55
1.35
2.88
1.58
1.81
1.62
1.47
1.40
1.39
···
···
1.73
1.64
1.37
−1.90
12.35
−6.03
6.95
−0.77
−4.24
−4.20
−4.47
−2.98
13.22
10.16
−0.61
6.92
11.80
−0.17
0.27
−2.44
−8.53
0.09
−12.77
4.00
−19.31
−13.92
−5.32
−5.27
4.02
6.16
5.20
6.77
6.37
3.95
−7.87
5.85
7.27
5.35
6.73
5.76
6.52
7.08
7.62
7.63
6.98
5.05
5.60
6.46
5.77
5.04
6.67
6.54
6.77
7.60
6.71
9.37
6.19
8.75
8.19
7.22
7.21
6.21
5.92
6.05
5.66
5.88
6.31
7.22
Sc
σS
0.0041
0.0044
0.0042
0.0040
0.0042
0.0040
0.0039
0.0040
0.0041
0.0041
0.0041
0.0043
0.0041
0.0041
0.0042
0.0042
0.0041
0.0046
0.0041
0.0041
0.0042
0.0042
0.0040
0.0040
0.0039
0.0041
0.0041
0.0040
0.0041
0.0040
0.0042
0.0042
0.0041
0.00005
0.00004
0.00002
0.00002
0.00002
0.00003
0.00003
0.00003
0.00003
0.00002
0.00003
0.00003
0.00003
0.00003
0.00003
0.00002
0.00003
0.00004
0.00002
0.00002
0.00003
0.00002
0.00004
0.00003
0.00003
0.00002
0.00002
0.00002
0.00001
0.00002
0.00002
0.00002
0.00002
a The fitted zero-point that is on an arbitrary scale (denoted as γ
rel in Tab. 4)
has been subtracted from the velocities.
b The values for σ
RV do not include the jitter.
c S values are on a relative scale.
weakly determined outer period due to the short 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 singlecomponent 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 5 m s−1 , 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
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 i 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).
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 Teff⋆ = 5760 ±
90 K, stellar surface gravity log g⋆ = 4.28 ± 0.10 (cgs),
metallicity [Fe/H] = +0.50 ± 0.08 dex, and projected rotational velocity v sin i = 1.9 ± 0.5 km s−1 .
Further analysis of the spectra has shown that the host
′
star is chromospherically quiet with log RHK
= −5.10
and S = 0.14 (on an absolute scale).
Bakos et al.
400
200
RV (m/s)
0
-200
-400
-600
-800
Residual
-1000
15
10
5
0
-5
-10
550
550
600
600
650
650
700
700
750
800
750
800
BJD - 2454000
850
850
900
900
950
950
1000
i-band i-band i-band i-band i-band i-band i-band
4
0
0.05
0.1
0.15
1000
0.2
200
RV (m/s)
0
-200
0.25
-400
-600
-800
0.3
Outer phase with respect to Tc2
-0.15
100
RV (m/s)
-100
Inner phase with respect to Tc1
Bisec. (m/s)
0
0.05
0.1
0.15
Fig. 3.— Unbinned instrumental i 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 Ntr = 0 (2008
November 8/9 MST, the third event from the top), then the followup light curves have Ntr = -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 shotnoise, plus the readout noise.
0
-50
S activity
-0.05
Time from transit center (days)
50
40
20
0
-20
-40
-60
0.006
0.005
0.004
0.003
0.002
0.001
-0.1
TABLE 2
Photometry follow-up of HAT-P-13
Inner phase with respect to Tc1
BJD
(2,400,000+)
0.0
0.2
0.4
0.6
0.8
Inner phase with respect to Tc1
1.0
Fig. 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 m s−1 , requiring a jitter of 3.0 m s−1 to be added in quadrature
to the individual errors to yield a reduced χ2 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 Tc2 , where c 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 Tc1 ). 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.
54581.66038
54581.66070
54581.66104
54581.66138
54581.66171
54581.66205
54581.66237
54581.66271
Maga
σMag
0.00550
0.00763
0.00877
0.00725
0.00670
0.01004
0.00673
0.00787
0.00080
0.00080
0.00080
0.00080
0.00080
0.00080
0.00080
0.00080
Mag(orig)
Filter
0.03620
0.03652
0.03686
0.03720
0.03753
0.03787
0.03819
0.03853
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.
a 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.
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 a/R⋆ (or related mean stellar density ρ⋆ ). We
note that another possible parameter in place of a/R⋆
would be the stellar surface gravity, but in the case of
planetary transits the a/R⋆ quantity typically imposes a
stronger constraint on the stellar models (Sozzetti et al.
2007). (The validity of this assumption, namely that the
HAT-P-13b,c
11 The reason for the intersecting isochrones is the “kink” on the
evolutionary tracks for M⋆ & 1.0 M⊙ stars evolving off the main
sequence.
3.0
4.0
5.0
6.0
a/R*
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 Teff⋆ , [Fe/H], and log g⋆ as derived from the SME analysis, we performed a global modeling of the data (§ 3.3),
yielding a full Monte-Carlo distribution of a/R⋆ . For
Teff⋆ and [Fe/H]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 (∼ 20000) set
of a/R⋆ , Teff⋆ , [Fe/H] values, we searched the stellar
isochrones of the Yi et al. (2001) models for the best fit
stellar model parameters (M⋆ , R⋆ , age, log g⋆ , 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 log g⋆ = 4.13 ± 0.04
based on the isochrone search, and only adjusting Teff⋆ ,
[Fe/H] and v sin i. 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 Teff⋆ = 5653 ± 90 K,
[Fe/H] = +0.41 ± 0.08 and v sin i = 2.9 ± 1.0 km s−1 .
We adopted these as the final atmospheric parameters
for the star. We then also repeated the isochrone search
for stellar parameters, obtaining M⋆ =1.219+0.050
−0.099 M⊙ ,
R⋆ =1.559 ± 0.082 R⊙ and L⋆ =2.22 ± 0.31 L⊙ . These are
summarized in Tab. 3, along with other stellar properties. Model isochrones from Yi et al. (2001) for metallicity [Fe/H]=+0.41 are plotted in Fig. 4, with the final
choice of effective temperature Teff⋆ and a/R⋆ marked,
and encircled by the 1-σ and 2-σ confidence ellipsoids.
The second SME iteration at fixed stellar surface gravity
(as determined from a/R⋆ ) changed the metallicity and
stellar temperature in such a way that the new (Teff⋆ ,
a/R⋆ ) 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).11 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 ∼ 1.23 M⊙ , and much smaller peak
around ∼ 1.13 M⊙. The asymmetric error-bars given
in Tab. 3 for the mass and age account for the doublepeaked distribution.
The stellar evolution modeling also yields the abso-
5
7.0
8.0
9.0
10.0
6000
5500
5000
Effective temperature [K]
Fig. 4.— Stellar isochrones from Yi et al. (2001) for metallicity
[Fe/H]=+0.41 and ages between 3.6 and 8.0 Gyr in steps of 0.4 Gyr.
The final choice of Teff⋆ and a/R⋆ are marked and encircled by
the 1-σ and 2-σ confidence ellipsoids. The open triangle denotes
the Teff⋆ ,a/R⋆ point found during the first SME iteration, before
refining the stellar surface gravity.
lute 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 J2MASS = 9.328 ± 0.018, H2MASS = 9.058 ± 0.017
and K2MASS = 8.975 ± 0.017, which we transformed to
JESO = 9.396±0.022, HESO = 9.070±0.021 and KESO =
9.018±0.018, following Carpenter (see 2001). These yield
a color of (J − K) = 0.378 ± 0.028, fully consistent with
the expected, isochrone-based (J −K)iso = 0.374±0.015.
We thus relied on the 2MASS K apparent magnitude and
the MK = 2.36 ± 0.12 absolute magnitude derived from
the above-mentioned modeling to determine the distance:
214 ± 12 pc, assuming no reddening. The K 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.
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
6
Bakos et al.
TABLE 3
Stellar parameters for HAT-P-13
Parameter
Value
Source
Teff⋆ (K). . . . . . .
[Fe/H] (dex) . . .
v sin i (km s−1 )
vmac (km s−1 ) .
vmic (km s−1 ) . .
γRV (km s−1 ) . .
γ1(i) . . . . . . . . . . .
γ2(i) . . . . . . . . . . .
M⋆ (M⊙ ) . . . . . .
R⋆ (R⊙ ) . . . . . . .
log g⋆ (cgs) . . . .
L⋆ (L⊙ ) . . . . . . .
V (mag) . . . . . . .
B − V (mag) . .
MV (mag) . . . . .
K (mag,ESO)
MK (mag,ESO)
Age (Gyr) . . . . .
Distance (pc) . .
5653 ± 90
+0.41 ± 0.08
2.9 ± 1.0
3.83
0.85
+14.69 ± 0.68
0.3060
0.3229
1.22+0.05
−0.10
1.56 ± 0.08
4.13 ± 0.04
2.22 ± 0.31
10.622
0.73 ± 0.03
3.97 ± 0.17
8.999 ± 0.017
2.36 ± 0.12
5.0+2.5
−0.7
214 ± 12
SMEa
SME
SME
SME
SME
DS
SME+Claret
SME+Claret
YY+a/R⋆ +SMEc
YY+a/R⋆ +SME
YY+a/R⋆ +SME
YY+a/R⋆ +SME
TASSd
YY+a/R⋆ +SME
YY+a/R⋆ +SME
2MASSe
YY+a/R⋆ +SME
YY+a/R⋆ +SME
YY+a/R⋆ +SME
a 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.
b SME+Claret = Based on SME analysis and tables
from Claret (2004).
c YY+a/R +SME = Yi et al. (2001) isochrones, a/R
⋆
⋆
luminosity indicator, and SME results.
d Based on the TASS catalogue (Droege et al. 2006).
e Based on the transformations from Carpenter (2001).
This section presents a simultaneous fitting of 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 i band, corresponding to the stellar properties determined from the SME analysis (§ 3.1). The transit
shape was parametrized by the normalized planetary radius p ≡ Rp /R⋆ , the square of the impact parameter
b2 , and the reciprocal of the half duration of the transit ζ/R⋆ . 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 Binst “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 K, and Lagrangian orbital
elements (k, h) = e × (cos ω, sin ω), plus an eccentric Keplerian orbit for the outer object with K2 , k2 and h2 ,
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 (TA and TB ), and all intermediate transits were interpolated using these two epochs
and their corresponding Ntr transit numbers. The model
for the RV data component also implicitly contained its
ephemeris information through TA and TB , 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
(∼ 0.5 km s−1) 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-P13b 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 ∼ 1000 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 Ntr = 0 to the first complete high quality followup light curve gathered on 2008 November 8/9 (MST).
The HATNet data-set was characterized by Tc,−370 and
Tc,−309 , covering all transit events observed by HATNet
(here the c subscript denotes “center” for the transits of
HAT-P-13b). The transit follow-up observations were
characterized by their respective times of transit center: Tc,−68 , Tc,−1 , Tc,0 , Tc,1 , Tc,24 , Tc,35 , Tc,62 . Initial
estimates of the Tc,i values yielded an initial epoch Tc
and period P 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 Tc and P 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), Tc
and P were refined, and the fit was repeated.
The time dependence of the RV of the outer planet
was described via its hypothetical transit time T2c (as
if its orbit were edge on), and its period P2 . The time
of periastron passage T2peri can be equivalently used in
place of time of conjunction T2c .
Altogether, the 21 parameters describing the physical model were Tc,−370 , Tc,−309 (HATNet), Tc,−68 , Tc,−1 ,
HAT-P-13b,c
Tc,0 , Tc,1 , Tc,24 , Tc,35 , Tc,62 (FLWO 1.2 m), Rp /R⋆ , b2 ,
ζ/R⋆ , K, γ, k = e cos ω, h = e sin ω, K2 , k2 , h2 , P2 and
T2c . Two additional auxiliary parameters were the instrumental blend factor Binst of HATNet, and the HATNet out-of-transit magnitude, M0,HATNet .
-0.01
-0.005
Relative i-mag
0
0.005
0.01
0.015
0.02
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
Time from transit center (days)
Fig. 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 χ2 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 (k = −0.016 ± 0.008,
h = 0.008 ± 0.012, implying e = 0.021 ± 0.009), and also
7
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 k and h 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: Tc,−370 = 2453700.9182 ± 0.0068 (BJD), Tc,−309 =
2453878.8216 ± 0.0093 (BJD), Binst = 0.82 ± 0.06, and
M0,HATNet = 9.7966 ± 0.0001 (I 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 m s−1 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
χ2 with 1.0 for the RV data is 3.0 m s−1 . The low jitter is
also consistent with HAT-P-13 being a chromospherically
′
quiet star (based on log RHK
and S).
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 Mp = 0.853+0.029
−0.046 MJ , the radius is Rp = 1.281±
−3
0.079 RJ, and its density is ρp = 0.498+0.103
. The
−0.069 g cm
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 m2 sin i2 = 15.2±1.0 MJ, and orbits the star in
a highly eccentric orbit with e2 = 0.691 ± 0.018. The orientation of the orbit (ω2 = 176.7 ± 0.5◦ ) 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 K2 , e2
and γ, leading to correlated m2 sin i and e2 (Fig. 6).
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 P > 1 yr orbit . With an iron
abundance of [Fe/H] = +0.41 ± 0.08, 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 perturba-
8
Bakos et al.
0.74
TABLE 4
Orbital and planetary parameters for HAT-P-13b
0.72
Parameter
0.7
Value
e2
Best fit
0.68
0.66
0.64
0.62
11
12
13
14
15
16
17
18
19
m2 sin(i)
Fig. 6.— The orbital parameters, and as a consequence, the
mass (m2 sin i) and eccentricity e2 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 m2 sin i is due to the similar
distribution of the stellar mass (see § 3.1).
tion, 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 ([Fe/H]=+0.41 ± 0.08) is XO2b ([F e/H] = +0.45, Burke et al. 2007). XO-2b, however, has much smaller mass (0.56 MJ ) and smaller radius
(0.98 RJ ) than HAT-P-13b (0.85 MJ , 1.28 RJ ). Other
planets similar to HAT-P-13b include HAT-P-9b (Mp =
0.78 MJ, Rp = 1.4 RJ , Shporer et al. 2008), and XO-1b
(Mp = 0.92 MJ, Rp = 1.21 RJ , 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.12 Given the fact, however, that the
host star is metal rich, and fairly old (5.0+2.5
−0.7 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, HATP-13b is broadly consistent only with a 300 Myr planet
at 0.02 AU solar distance and core-mass up to 25 M⊕ ,
or a 1 Gyr core-less (pure H/He) planet. If HAT-P-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
12 The equivalent semi-major axis, at which HAT-P-13b would
receive the same amount of insolation when orbiting our Sun, is
aequiv = 0.0285 ± 0.0016 AU.
Light curve parameters, HAT-P-13b
P (days) . . . . . . . . . . . . . . . . . . . .
2.916260 ± 0.000010
Tc (BJD) . . . . . . . . . . . . . . . . . . .
2454779.92979 ± 0.00038
T14 (days) a . . . . . . . . . . . . . . . . .
0.1345 ± 0.0017
T12 = T34 (days) a . . . . . . . . . .
0.0180 ± 0.0018
a/R⋆ . . . . . . . . . . . . . . . . . . . . . . . .
5.84 ± 0.26
ζ/R⋆ . . . . . . . . . . . . . . . . . . . . . . . .
17.07 ± 0.16
Rp /R⋆ . . . . . . . . . . . . . . . . . . . . . .
0.0844 ± 0.0013
b2 . . . . . . . . . . . . . . . . . . . . . . . . . . .
0.447+0.044
−0.056
b ≡ a cos i/R⋆ . . . . . . . . . . . . . . .
0.668+0.032
−0.045
i (deg) . . . . . . . . . . . . . . . . . . . . . .
83.4 ± 0.6
Tperi (days) . . . . . . . . . . . . . . . . .
2454780.64 ± 0.42
RV parameters, as induced by HAT-P-13b
K (m s−1 ) . . . . . . . . . . . . . . . . . .
k ............................
h ...........................
e ............................
ω ...........................
106.1 ± 1.4
−0.016 ± 0.008
0.008 ± 0.012
0.021 ± 0.009
181 ± 45◦
Other RV parameters
γrel (m s−1 ) . . . . . . . . . . . . . . . . .
Jitter (m s−1 ) . . . . . . . . . . . . . . .
−51.3 ± 3.8
3.0
Secondary eclipse parameters for HAT-P-13b
Ts (BJD) . . . . . . . . . . . . . . . . . . .
2454781.359 ± 0.014
Ts,14 . . . . . . . . . . . . . . . . . . . . . . . .
0.1345 ± 0.0069
Ts,12 . . . . . . . . . . . . . . . . . . . . . . . .
0.0180 ± 0.0018
Planetary parameters for HAT-P-13b
Mp (MJ ) . . . . . . . . . . . . . . . . . . . .
Rp (RJ ) . . . . . . . . . . . . . . . . . . . . .
C(Mp , Rp ) b . . . . . . . . . . . . . . . .
ρp (g cm−3 ) . . . . . . . . . . . . . . . . .
a (AU) . . . . . . . . . . . . . . . . . . . . . .
log gp (cgs) . . . . . . . . . . . . . . . . . .
Teq (K) . . . . . . . . . . . . . . . . . . . . .
Θ c .........................
Fper (erg s−1 cm−2 ) d . . . . . . .
Fap (erg s−1 cm−2 ) e . . . . . . . .
hF i (erg s−1 cm−2 ) . . . . . . . . . .
0.853+0.029
−0.046
1.281 ± 0.079
0.56
0.498+0.103
−0.069
0.0427+0.0006
−0.0012
3.11 ± 0.05
1653 ± 45
0.046 ± 0.003
(1.76 ± 0.20) · 1013
(1.62 ± 0.18) · 1013
(1.69 ± 0.18) · 1013
a T : total transit duration, time between first to last contact;
14
T12 = T34 : ingress/egress time, time between first and second,
or third and fourth contact.
b Correlation coefficient between the planetary mass M and
p
radius Rp .
c The Safronov number is given by Θ = 1 (V
2
esc /Vorb ) =
2
(a/Rp )(Mp /M⋆ ) (see Hansen & Barman 2007).
d Incoming flux per unit surface area.
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 nonzero only at the 2-σ level (because k is non-zero at 2-σ),
its pericenter is aligned (to within 4◦ ± 40◦ ) 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 par-
HAT-P-13b,c
TABLE 5
Orbital and planetary parameters for HAT-P-13c
Parameter
Value
RV parameters, as induced by HAT-P-13c
P2 (days) . . . . . . . . . . . . . . . . . . .
428.5 ± 3.0
T2c a (BJD) . . . . . . . . . . . . . . . . .
2454870.4 ± 1.8
K2 (m s−1 ) . . . . . . . . . . . . . . . . .
502 ± 37
k2 . . . . . . . . . . . . . . . . . . . . . . . . . .
−0.690 ± 0.018
h2 . . . . . . . . . . . . . . . . . . . . . . . . . .
0.039 ± 0.005
e2 . . . . . . . . . . . . . . . . . . . . . . . . . . .
0.691 ± 0.018
ω2 . . . . . . . . . . . . . . . . . . . . . . . . . .
176.7 ± 0.5◦
T2,peri (days) . . . . . . . . . . . . . . .
2454890.05 ± 0.48
Fictitious light curve parameters, HAT-P-13cb
T2,14 c (days) . . . . . . . . . . . . . . . .
0.621 ± 0.029
T2,12 = T34 (days) . . . . . . . . . .
0.0355 ± 0.0012
Fictitious secondary eclipse parameters for HAT-P-13ca
T2s (BJD) . . . . . . . . . . . . . . . . . .
2454912.7 ± 1.8
T2s,14 (days) . . . . . . . . . . . . . . . .
0.574 ± 0.028
T2s,12 (days) . . . . . . . . . . . . . . . .
0.0355 ± 0.0012
Planetary parameters for HAT-P-13c
m2 sin i2 (MJ ) . . . . . . . . . . . . . .
15.2 ± 1.0
a2 (AU) . . . . . . . . . . . . . . . . . . . . .
1.188+0.018
−0.033
T2,eq (K) . . . . . . . . . . . . . . . . . . . .
340 ± 9
−1
−2
d
F2,per (erg s cm )
.....
(2.26 ± 0.35) · 1011
F2,ap (erg s−1 cm−2 ) d . . . . . .
(7.61 ± 0.86) · 109
hF2 i (erg s−1 cm−2 ) d . . . . . . .
(3.00 ± 0.33) · 1010
a T
2c would be the center of transit of HAT-P-13c, if its
(unknown) inclination was 90◦ .
b Transits of HAT-P-13c have not been observed. The
values are for guidance only, and assume zero impact parameter.
c T : total transit duration, time between first to last
14
contact, assuming zero impact parameter. T12 = T34 :
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.
d Incoming flux per unit surface area in periastron, apastron, and averaged over the orbit.
ticular, Batygin, Bodenheimer, & Laughlin (2009) give
limits on the tidal Love number k2b , core-mass and tidal
energy dissipation rate Qb 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 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 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
9
we refer to the time difference between the observed transit center (Tc,i ) and the calculated value based on a fixed
epoch and period as given in Tab. 4, i.e. in the O − C
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 Ntr = 0
(Tc,−1 ,Tc,0 ,Tc,1 ) show a standard deviation that is comparable to the error-bars (and we can safely assume zero
TTV within a ±2.9163 day time range). It is also apparent from the plot that the smallest error-bars correspond
to the full transit observations (Ntr = 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:
M⋆ , Mp , P , ip , P2 , e2 , ω2 , m2 sin i2 and T2,peri (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, i2 , and the relative angle of
the orbital normals projected onto the plane of the sky,
D = Ω1 − Ω2 (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-P13b (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 i2 =
83.4 ± 0.6◦ (corresponding to co-planar inner and outer
orbits, and to M2 = 15.2 MJ), and i2 = 8.1◦ (an ad
hoc value yielding an almost face-on orbit with M2 =
105 MJ), and D = 0◦ or D = 90◦ . These four analytic
models are illustrated in Fig. 7. The bottom panel of
Fig. 7 shows the TTVl and TTVg effects separately for
the i2 = 83.4±0.6◦ and D = 0◦ case. The i2 = 83.4±0.6◦
(M2 = 15.2 MJ) cases give TTV variations of the order of
15 seconds. The functional form of the i2 = 83.4 ± 0.6◦
and D = 0◦ 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 i2 at constant
m2 sin i2 does increase the TTV amplitude up to 100 seconds at i2 = 8.1◦ , 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 (∼
10
Bakos et al.
TABLE 6
Transit timing variations of
HAT-P-13.
Ntr
Tc
(BJD)
σTc
(sec)
O-C
(sec)
−68
−1
0
1
24
35
62
2454581.62406
2454777.01287
2454779.92953
2454782.84357
2454849.92062
2454882.00041
2454960.73968
105.7
86.7
54.3
133.6
65.1
129.5
154.1
−10.9
−58.5
−24.2
−215.6
50.7
132.2
156.1
100 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 i2 and D, i.e., HATP-13c may become an RV-based detection with known
orbital orientation and a true mass (even if it does not
transit).
i2=i1,Ω1-Ω2=0
i2=i1,Ω2-Ω1=90
i2=8.1,Ω2-Ω1=0
i2=8.1,Ω2-Ω1=90
200
TTV [sec]
100
0
-100
-200
-80
-60
-40
-20
0
20
40
60
80
40
60
80
Ntr transit number
40
TTVg + TTVl
TTVg
TTVl
TTV [sec]
20
0
-20
-40
-80
-60
-40
-20
0
20
Ntr transit number
Fig. 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 i2 = 83.1◦ = i1 , ii) i2 = 8.1◦ (at fixed m2 sin i2 ,
corresponding to a large mass for HAT-P-13c), iii) i2 = 83.1◦ and
the mutual inclination of the orbital normals in the plane of the
sky is D = Ω1 − Ω2 = 90◦ (see Fig. 1 of Borkovits et al. (2003)), or
iv) i2 = 8.1◦ and D = 90◦ . (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).
4.3. The planet HAT-P-13c
The probability of HAT-P-13c transiting the host star,
as seen from the Earth, depends on R⋆ , R2p , e2 , ω2 ,
and a2 (Kane & von Braun 2009). We evaluated the
transit probability in a Monte Carlo fashion, as part of
the global modeling, resulting in P2,tr = 0.0130 ± 0.0008
if R2,p = 1 RJ (P2,tr = 0.0122 ± 0.0007 if R2,p → 0).
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 ∼ 428 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 15.2 MJ , we have good reasons to believe
that its radius would be around 1RJ , based on heavymass transiting planets like HAT-P-2b (8.84 MJ , Bakos
et al. 2007), XO-3b (11.79 MJ , Johns-Krull et al. 2008),
Corot-3b (21.66 MJ , Deleuil et al. 2008), all having radii
around 1 RJ . 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 14.9 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 e2 and ω2 .
As regards to the nature of HAT-P-13c, even at its
minimal mass (15.2 ± 1.0 MJ ) 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 P > 1 yr,
e = 0.7 orbits around metal-rich ([Fe/H] ≈ 0.5) 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 Hun-
HAT-P-13b,c
garian 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
11
(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.
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ács-Dajka, 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
View publication stats
Kovács, G., Bakos, G. Á., & Noyes, R. W. 2005, MNRAS, 356, 557
Johns-Krull, 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