Photometric follow-up of the transiting planetary system TrES-3: transit timing variation and long-term stability of the system^★

Abstract

INTRODUCTION

Since the first discovery of a transiting planet around HD 209458 (Charbonneau et al. 2000), 235 transiting extrasolar systems have already been confirmed up to 2012 December 4.¹ Whilst transiting exoplanets offer unique scientific possibilities, their study involves several complications. In general, it is impossible to measure the mass and radius of a planet based on a data set obtained with one observational technique. Transit light curves allow us to determine just the relative size of a star and planet, the orbital inclination and the stellar limb-darkening coefficients. By combining this with radial-velocity measurements, the observations offer the opportunity to measure the precise stellar and planetary parameters. In order to obtain such parameters, some constraints are needed, and are usually provided by forcing the properties of the host stars to match theoretical expectations (Southworth 2010). Significant uncertainties remain in the stellar mass and radius determinations of many systems. In some cases, this is due to poorly determined photospheric properties (i.e. effective temperature and metallicity), and in other cases due to a lack of an accurate luminosity estimate (Sozzetti et al. 2009). In addition, the different methods used for these determinations as well as different approaches towards systematic errors are leading to rather inhomogeneous set of planet properties. Because of such inhomogenities, recently a few papers were published where authors re-analysed a large subset of known transiting planets, and applied a uniform methodology to all systems (e.g. Torres, Winn & Holman 2008; Southworth 2010, 2012).

In this paper, we focus on the transiting system TrES-3. The system consists of a nearby G-type dwarf and a massive hot Jupiter with an orbital period of 1.3 d. It was discovered by O’Donovan et al. (2007) and also detected by the SuperWASP (Wide Angle Search for Planets) survey (Collier Cameron et al. 2007). Later, Sozzetti et al. (2009) presented new spectroscopic and photometric observations of the host star. A detailed abundance analysis based on high-resolution spectra yields [Fe/H] = −0.19 ± 0.08, T_eff = 5650 ± 75 K and log g = 4.0 ± 0.1. The spectroscopic orbital solution was improved with new radial velocity measurements obtained by Sozzetti et al. (2009). Moreover, these authors redetermined the stellar parameters (i.e. M_* = 0.928|$^{+0.028}_{-0.048}$| M_⊙ and R_* = 0.829|$^{+0.015}_{-0.022}$| R_⊙) and finally, the new values of the planetary mass and radius were determined (see Table 3). They also studied the transit timing variations (TTVs) of TrES-3 and noted significant outliers from a constant period. In the same year, Gibson et al. (2009) presented the follow-up transit photometry. It consisted of nine transits of TrES-3, taken as part of transit timing programme using the RISE instrument on the Liverpool Telescope. These transits, together with eight transit times published before (Sozzetti et al. 2009), were used to place upper mass limit as a function of the period ratio of a potential perturbing planet and transiting planet. It was shown that timing residuals are sufficiently sensitive to probe sub-Earth mass planet in both interior and exterior 2:1 resonances, assuming that the additional planet is in an initially circular orbit. Christiansen, Ballard & Charbonneau (2011) have observed TrES-3 as a part of the NASA EPOXI Mission of Opportunity. They detected a long-term variability in the TrES-3 light curve, which may be due to star spots. They also confirmed that the planetary atmosphere does not have a temperature inversion. Later, Turner et al. (2013) observed nine primary transits of the hot Jupiter TrES-3b in several optical and near-UV photometric bands from 2009 June to 2012 April in an attempt to detect its magnetic field. Authors determined an upper limit of TrES-3b's magnetic field strength between 0.013 and 1.3 G using a timing difference of 138 s derived from the Nyquist–Shannon sampling theorem. They also presented a refinement of the physical parameters of TrES-3b, an updated ephemeris and its first published near-UV light curve. The near-UV planetary radius of R_p = 1.386|$^{+0.248}_{-0.144}$|R_J was also determined. This value is consistent with the planet's optical radius. Recently, Kundurthy et al. (2013) observed eleven transits of TrES-3b over a two-year period in order to constrain system parameters and look for transit timing and depth variations. They also estimated the system parameters for TrES-3b and found consistency with previous estimates. Their analysis of the transit timing data show no evidence for transit timing variations and timing measurements are able to rule out Super-Earth and Gas Giant companions in low order mean motion resonance (MMR) with TrES-3b.

The main aims of this study can be summarized in the following items: (i) determination of the system parameters for TrES-3b (two independent codes will be used) and comparison with previous studies (i.e. O’Donovan et al. 2007; Gibson et al. 2009; Sozzetti et al. 2009; Colón et al. 2010; Southworth 2010, 2011; Christiansen et al. 2011; Lee et al. 2011; Sada et al. 2012; Turner et al. 2013; Kundurthy et al. 2013). (ii) Based on the obtained transits, we will determine the mid-transit times (T_C) and with following analysis of TTV we will discuss possible presence of a hypothetical additional planet (perturber). We will try to estimate its upper-mass limit as a function of orbital periods ratio of transiting planet and the hypothetical perturber. (iii) Finally, using the long-term integration and applying the method of maximum eccentricity we will search for stability of regions inside the TrES-3b planet in context of additional planet(s).

The remainder of this paper is organized as follows. In the Section 2, we describe observations and data reduction pipelines used to produce the light curves. Section 3 presents the methods for analysis of transit light curves as well as discussion and comparison of the parameters of TrES-3 system. Sections 4 and 5 are devoted to TTV and long-term stability of the system, respectively. Finally, in Section 6 we summarize and discuss our results.

OBSERVATIONS AND DATA REDUCTION

We obtained our data using several telescopes with different instruments. This allowed us to obtain many light curves since this strategy can effectively cope with the weather problems. On the other hand, this approach results in rather heterogeneous data. We used most of the data in average quality for the TTV analysis which is not very demanding on homogeneity of the data. Only the best quality light curve was used for the planet parameter determination.

Observations used in this paper were carried out at the several observatories in Slovakia (Stará Lesná Observatory; 49° 09^′ 10^′′ N, 20° 17^′ 28^′′ E), Poland (Piwnice Observatory; 53° 05^′ 43^′′ N, 18° 13^′ 46^′′ E), Germany (Grossschwabhausen Observatory; 50° 55^′ 44^′′ N, 11° 29^′ 03^′′ E; Volkssternwarte Kirchheim Observatory, 50° 55^′ 44^′′ N, 11° 29^′ 03^′′ E and Michael Adrian Observatory, 49° 55^′ 27^′′ N, 08° 24^′ 33^′′ E) and Spain (Calar Alto Observatory; 37° 13^′ 25^′′ N, 02° 32^′ 46^′′ E). We collected 14 transit light curves obtained between 2009 May and 2011 September. The transits on 2009 May 12 and 2010 August 20 were observed simultaneously at two different observatories. The telescope diameters of 0.5 to 2.2 m allowed us to obtain photometry with 1.2–7.8 mmag precision, depending on observing conditions. Observations generally started ∼1 h before the expected beginning of a transit and ended ∼1 h after the event. Unfortunately, weather conditions and schedule constraints meant that we were not able to fit this scheme in all cases.

All instruments are equipped with CCD cameras with the Johnson-Cousins (UBVR_CI_C) standard filter system. The information from individual observatories and instruments as well as the summary of observing runs are given in Tables 1 and 2. The standard correction procedure (bias, dark and flat-field correction) and subsequently aperture photometry was performed by iraf² and task chphot (Raetz et al. 2009) (GSH and VK), c-munipack package³ (G1) and mira_pro_7⁴ (MA). Data from remaining telescopes (P and CA) were reduced with the software pipeline developed for the Semi-Automatic Variability Search sky survey (Niedzielski, Maciejewski & Czart 2003). To generate an artificial comparison star, at least 20–30 per cent of stars with the lowest light-curve scatter were selected iteratively from the field stars brighter than 2.5–3 mag below the saturation level (e.g. Broeg, Fernández & Neuhäuser 2005; Raetz et al. 2009). To measure instrumental magnitudes, various aperture radii were used. The aperture which was found to produce light curve with the smallest overall scatter was applied to generate final light curve. The linear trend in the out-of-transit parts was also removed.

LIGHT-CURVE ANALYSIS

Our photometric observations of the TrES-3 system consist of data from different instruments and are of different photometric quality. For the purpose of radius determination, we decided to analyse only a light curve with the lowest scatter and the best rms = 1.2 mmag. We choose data obtained at Calar Alto observatory on 2010 September 06 (see Figs 1 and 2). We first refined the light-curve system parameters and subsequently determined the individual times of transit as described below in order to improve ephemeris.

For the calculation of synthetic light curves we used two independent approaches: the first one is based on the routines from Mandel & Agol (2002) and the Monte Carlo (MC) simulations (described in Section 3.1), the second one uses jktebop code (Southworth, Maxted & Smalley 2004) (see Section 3.2).

SOLUTION 1

Figure 1.

Open in new tabDownload slide

The light curves of the system TrES-3 obtained at individual observatories between 2009 May and 2011 September. The best light curve (obtained at Calar Alto, 2010 September) was used to determine system parameters. The data from other observatories were used for T_C determination and TTV analysis. The typical error bar is plotted next to each light curve.

Figure 2.

Open in new tabDownload slide

Top: the light curve obtained at Calar Alto on 2010 September 06, and the best fit corresponding by SOLUTION 1. Middle: residuals from the best fit mentioned above. Bottom: difference from the fits corresponding by SOLUTION 1 and SOLUTION 2.

Figure 3.

Open in new tabDownload slide

Top: reduced |$\chi ^2_{\mathrm{r}}$| representing global minimum solution for calculated parameter R_p/R_*. Bottom: reduced |$\chi ^2_{\mathrm{r}}$| representing global minimum solution for calculated parameter i.

Figure 4.

Open in new tabDownload slide

Confidence region depicted as a projection of the 4-dimensional region into 2-dimensional parameter space (inclination and ratio of the planet to star radius). Regions of 1σ, 2σ and 3σ corresponding to Δχ = 4.72, 9.7 and 16.3, respectively, are marked with different grey colours.

Finally, we took into account the uncertainty of the stellar mass and semimajor axis according to the simple error propagation rule. The results from the first analysis are shown in the Table 3 as SOLUTION 1.

SOLUTION 2 (jktebop code)

We have modelled the light curve using the jktebop⁵ code (Southworth 2011) as well. jktebop grew out of the original ebop program written for eclipsing binary systems (Etzel 1981; Popper & Etzel 1981) via implementing the NDE (Nelson & Davis 1972) model. jktebop uses biaxial spheroids to model the component stars (or star and planet) and performs a numerical integration in concentric annuli over the surface of each body to obtain the flux coming from the system. This feature of the code allows us to avoid small and spherical planet approximations which are used in analytic light-curve generators based on Mandel & Agol (2002), and hence to derive planet's oblateness. A model is fitted to the data by the Levenberg–Marquardt least-square procedure. The code converges rapidly towards a reliable solution and diminishes the correlation between fitted parameters (Southworth 2008).

The errors of derived parameters were determined in two ways for each combination of data set and adopted limb darkening (LD) law (we have used the same coefficients like in the case of SOLUTION 1). First, we ran 1000 MC simulations, a spread range of a given parameter within 68.3 per cent was taken as its error estimate. Secondly, the prayer-bead method (e.g. Winn et al. 2009; Désert et al. 2011) was used to check whether red noise was present in our data. MC errors were found to be 2–3 times smaller than the values returned by the prayer-bead method. This indicates that the light curve is affected not only by Poisson noise, but also by additional correlated noise. Therefore, our prayer-bead error estimates were taken as our final errors.

Light-curve analysis results

The resulting values of the parameters together with their uncertainties are given in Table 3. For comparison, the parameters from previous studies are added there as well. Fig. 2 shows the resulting best fit obtained by SOLUTION 1. In order to compare both solutions, we also plotted differences from the fit obtained by routines resulting from SOLUTION 1 and SOLUTION 2.

We used MC simulations to produce of 10 000 synthetic data sets with the same probability distribution as the residuals of the fit. From each synthetic data set we estimated the synthetic transit parameters. Using the results of the simulation, the dependence of the inclination versus the planet to star radius ratio is plotted in Fig. 4. The figure demonstrates that there is a correlation between the parameters and that solution for the system TrES-3 is not unique and can be located in a relatively wide range (degeneracy of the parameters). This fact is also caused by grazing transit and by subsequent sensitivity of the solution to the determination of LD coefficients.

For determination of our mid-transit times, the best-fitting model was used as a template and fitted to other light curves for which only the mid-transit time was allowed to vary. In order to have homogeneous analysis of TTV, we collected the data points of transit light curves published in previous studies (see Fig. 5). Transit light curve data of Sozzetti et al. (2009) and Colón et al. (2010) were available online and the remaining tabulated data were obtained from mentioned authors by e-mail. We re-analysed all transit light curves and derived mid-transit times by the same procedure (the code jktebop) to get a homogeneous data set for TTV analysis (see Section 4).

We also re-analysed all light curves under consideration with k = r_b/r_A as a free parameter and tried to search for any variation. All determinations of k were consistent within error bars with a mean value, so we did not detect any transit depth variation. We also saw no significant signal in periodograms.

TRANSIT TIMING VARIATION

Our new 14 mid-transit times and 42 redetermined literature values were used to refine the transit ephemeris. The mid-transit times were transformed from JD or HJD (based on UTC) into BJD (based on Barycentric Dynamical Time – TDB) using the online converter⁶ by Eastman, Siverd & Gaudi (2010). As a result of fitting a linear function of the epoch, we obtained the mid-transit time for the initial epoch T₀ = 245 4538.581 44 ± 0.000 07 BJD_TDB and the orbital period P_b = 1.306 185 99 ± 0.000 000 23 d. The individual mid-transit errors were taken as weights. The linear fit yields reduced χ² of 2.5, that is similar to a value of 2.3 reported by Gibson et al. (2009). These values, noticeable greater than 1 might suggest the existence of an additional planet which perturbs the orbital motion of TrES-3b (Gibson et al. 2009; Sozzetti et al. 2009).

Results for new mid-transit times are shown in Table 4. The observed minus calculated (O−C) diagram, plotted in Fig. 5, shows no significant deviation from the linear ephemeris. All our data points deviated by less than 3σ. We also searched for a periodicity that could be a sign of an additional body in the system. We generated a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) for the residuals in a frequency domain limited by the Nyquist frequency and found the highest peak at ∼30 d and peak-to-peak amplitude of 70 ± 20 s. This period could coincide with a stellar rotation period which is roughly estimated to be ∼28 d. Examples of such TTV signals induced by the stellar activity are observed in the Kepler data (e.g. Mazeh et al. 2013). However, the false alarm probability (calculated empirically by a bootstrap resampling method with 10⁵ trials) of the putative signal is disqualifying with a value of 18.2 per cent. In addition, the amplitude is close to the mean 1σ timing error of our observations. These findings allow us to conclude that a strictly periodic TTV signal with the amplitude greater than ∼1 min over a 4 yr time span seems to be unlikely.

Figure 5.

Open in new tabDownload slide

Observation minus calculation O − C diagram for transit timing of TrES-3b, plotted according to the new linear ephemeris. Open circles denote re-analysed mid-transit times from the literature. Filled symbols mark new mid-transit times reported in this paper.

Following Gibson et al. (2009), we put upper constraints on a mass of a potential perturbing planet in the system with refined assumptions. We adopted a value of 1 min for the maximal amplitude of the possible TTV signal. We also increased the sampling resolution to probe resonances other than inner and outer 1:2 commensurabilities studied in detail by Gibson et al. (2009). We simplified the three-body problem by assuming that the planetary system is coplanar and initial orbits of both planets are circular. The masses of both the star and the transiting planet, as well as its semimajor axis were taken from Sozzetti et al. (2009). The orbital period of the hypothetical perturber was varied in a range between 0.2 and 5 orbital periods of TrES-3b. We produced 2800 synthetic O−C diagrams based on calculations done with the mercury code (Chambers 1999). We applied the Bulirsch–Stoer algorithm to integrate the equations of motion.

The most important feature of the Bulirsch–Stoer algorithm for N-body simulations is that it is capable of keeping an upper limit on the local errors introduced due to taking finite time-steps by adaptively reducing the step size when interactions between the particles increase in strength. Calculations covered 1500 d, i.e. a total time span of available transit observations. The results of simulations are presented in Fig. 6. Our analysis allows us to exclude an additional Earth-mass planet close to inner 1:3, 1:2 and 3:5 and outer 5:3, 2:1 and 3:1 MMRs.

Figure 6.

Open in new tabDownload slide

The upper-mass limit of a hypothetical additional planet that could perturb the orbital motion of TrES-3b as a function of ratio of orbital periods of transiting planet, P_b, and the hypothetical perturber, P_p. Orbits located in a grey area were found to be unstable due to close encounters of both planets.

LONG-TERM STABILITY OF THE SYSTEM

In this section, we investigated the long-term gravitational influence of TrES-3b planet on a potential second planet in the system. Thus, we performed numerical simulation for studying the stability of orbits and checking their chaotic behaviour using the method of maximum eccentricity (e.g. Dvorak et al. 2003).

We used long-term integration of small-mass (Earth-mass) planet orbits for inspecting the stability regions in the TrES-3b system. The time span of the integration was around 140 000 revolutions of the planet around the star. For the integration of this system, we again applied an efficient variable-time-step algorithm (Bulirsch–Stoer integration method). The parameter ϵ which controls the accuracy of the integration was set to 10⁻⁸, in this case.

We have generated 10⁵ massless particles for representing small planets in this system. We assumed that semimajor axis ranges from 0.006 25 au to 0.1 au. The inner border of the generated system of small planets is around 1.5 times star radius. This value of inner border is approximately at the location of Roche limit for an Earth-like planet and includes also the region in the vicinity of 3:1 MMR (0.01 au for this case) in our analysis. The step size in semimajor axis was Δ a = 0.000 25 au, in eccentricity Δ e = 0.02 and Δ i = 5° in inclination. The upper limit of the grid in eccentricity was 0.5 and 50° in inclination. We integrated the orbits of the small planets for the time span of 500 yr (about 140 000 revolutions of TrES-3b around the parent star). We obtained orbital evolution of each small planet in the system and also were able to study the stability regions at the end of the integration using the method of maximum eccentricity, mentioned above. The maximum eccentricity that the potential perturber of TrES-3b reached over the time of the integration is plotted on the a − e (a − i) stability map in the Fig. 7 (Fig. 8).

Figure 7.

Open in new tabDownload slide

Stability plot in the a − e plane showing the maximum eccentricity. From top left to bottom right: i = 5°, i = 20°, i = 35° and i = 50°. The MMRs with TrES-3b planet are also marked. The minimal value of semimajor axis is 0.006 25 au.

Figure 8.

Open in new tabDownload slide

Stability plot in the a − i plane showing the maximum eccentricity. From top left to bottom right: e = 0.001, e = 0.1, e = 0.3 and e = 0.5. The MMRs with TrES-3b planet are also marked. The minimal value of semimajor axis is 0.006 25 au.

Fig. 7 shows the stability maps in a − e plane for selected values of inclinations: i = 5°, i = 20°, i = 35° and i = 50° (from top left to bottom right). Also the MMRs with TrES-3b planet are marked in these plots. We can see a stable region inside the 2:1 MMR for each value of inclination. Outside the 2:1 MMR, the gravitational influence of the TrES-3b planet is very strong and leads to depletion of these regions. Only the MMRs 2:1, 3:2 and 1:1 are moderately populated but the population decreases with the increase of the inclination. For completeness, we note that 5:3 (a ≈ 0.016 au) and 3:5 (a ≈ 0.032 au) MMRs are depleted and thus not stable for additional Earth-mass planet in their vicinity. In Fig. 8, we present the stability maps in a − i plane for several values of eccentricities: e = 0.001, e = 0.1, e = 0.3 and e = 0.5. One can see the stable region inside the 2:1 MMR with TrES-3b planet, which we found stable also in the a − e plane. Considering the region beyond the 1:3 MMR, the depletion of the planet population is not so strong than in the region between 2:1 and 1:3 MMR. This feature is in good agreement with the weaker gravitational influence of the TrES-3b planet.

Based on our results, we showed that the region inside the TrES-3b planet orbit (especially inside the 2:1 up to 3:1 MMR) can be stable on longer time-scales (hundreds of years or hundred thousands of revolutions of TrES-3b around the star) from the dynamical point of view. The region from 0.015 au (near 2:1 MMR) to 0.05 au (near 1:3 MMR) was found unstable apart from moderately populated MMRs located in this area (see Figs 7 and 8). The relatively small increase of population in the mentioned MMRs depends on the initial values of semimajor axis, eccentricity and inclination. The region beyond the 0.05 au was found to have a chaotic behaviour and the depletion of the planet population increases with increasing values of initial eccentricity as well as inclination.

DISCUSSION AND CONCLUSION

Based on the transit light curves obtained at several observatories between 2009 May and 2011 September, we redetermined orbital parameters and the radius of the transiting planet TrES-3b. The best light curve (obtained at Calar Alto, 2010 September) was used for light-curve analysis, and the data from the other observatories were used for T_C determination and TTV investigation. We used two independent solutions for parameters determination and finally, we concluded that our values are consistent with previous results of Sozzetti et al. (2009) and Christiansen et al. (2011).

The aim of this paper was also to discuss possible presence of a second planet in the extrasolar system TrES-3. For this purpose, we used our new 14 mid-transit times and the individual determinations from Sozzetti et al. (2009) and Gibson et al. (2009). The resulting O−C diagram showed no significant deviation of data points from the linear ephemeris. In addition, we tried to search for a periodicity that could be caused by an additional body in the system. We can conclude that a strictly periodic TTV signal with the amplitude greater than 1 min over a 4yr time span seems to be unlikely. This result, together with refined assumptions of Gibson et al. (2009) allow us to put upper constraints on the mass of a potential perturbing planet. The additional Earth-mass planet close to inner 3:1, 2:1 and 5:3 and outer 3:5, 1:2 and 1:3 MMRs can be excluded.

Finally, we used the long-term integration of the theoretical set of massless particles generated in TrES-3 system for studying the dynamical stability of potential second planet in the system (influenced by TrES-3b gravitation). From our analysis, we found that the region inside the TrES-3b planet orbit (especially inside the 2:1 MMR) up to 3:1 MMR can be stable on longer time-scales (hundreds of years or hundred thousands of revolutions of TrES-3b around the star) from the dynamical point of view. The region from 0.015 au (near 2:1 MMR) to 0.05 au (near 1:3 MMR) was found to be unstable apart from moderately populated MMRs located in this area. The relatively small increase of population in these MMRs depends on the initial values of semimajor axis, eccentricity and inclination. The region beyond 0.05 au was found to have a chaotic behaviour and depletion of the planet population increases with increasing values of initial eccentricity as well as inclination.

This work has been supported by the VEGA grants No. 2/0078/10, 2/0094/11, 2/0011/10. MV, MJ, JB, TP and ŠP would also like to thank the project APVV-0158-11. TK thanks the Student Project Grant at MU MUNI/A/0968/2009 and the National scholarship programme of the Slovak Republic. GM acknowledges Iuventus Plus grants IP2010 023070 and IP2011 031971. SR thanks the German National Science Foundation (DFG) for support in project NE 515/33-1.

REFERENCES