Free Access
Issue
A&A
Volume 571, November 2014
Article Number A38
Number of page(s) 13
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201424080
Published online 04 November 2014

© ESO, 2014

1. Introduction

Now, more than 1779 planets1 have been discovered and confirmed in about 1102 planetary systems. Around 460 planetary systems are known to be multiple planet systems. Hundreds of Kepler planetary candidates with multiple transit-like signals are still waiting confirmation (see Latham et al. 2011; Lissauer et al. 2011b). The usual way to characterize multiple planet systems is by combining information from both transits and radial velocities. An effect due to the presence of multiple planets is the transit time variation (TTV): the gravitational interaction between two planets causes a deviation from the Keplerian orbit and, as a consequence, the transit times (T0s) of a planet may be not strictly periodic (see Agol et al. 2005; Holman & Murray 2005; Miralda-Escudé 2002). This effect can be also exploited to infer the presence of an unknown planet, even if it does not transit the host star (Agol et al. 2005; Holman & Murray 2005). For example, the Kepler Transit Timing Observations series (TTO, Ford et al. 2011, and references therein) and TASTE project (Nascimbeni et al. 2011, and references therein) demonstrate the use of this technique.

The problem of determining the masses and the orbital parameters of the planets in a multiple system is a difficult inverse problem. In some works, the authors adopted an analytic approach to the problem (e.g., Nesvorný & Morbidelli 2008; Nesvorný 2009), developing a method from the perturbation theory (Hori 1966; Deprit 1969), where the T0s are computed as a Fourier series. A drawback of this method is that it does not take into account for the so-called mean motion resonance (MMR) or cases that are just outside the MMR. We have a MMR when the ratio of the periods of two planets is a multiple of a small integer number, such as 2:1. Two planets reciprocally in MMR, or just outside it, show a strong TTV signal that is easily detectable even with ground-base facilities.

The method described in this paper is based on a direct numerical N-body approach (Steffen & Agol 2005; Agol & Steffen 2007), which is conceptually simpler, but computationally intensive. Very recently, Deck et al. (2014) have developed TTVFast, a symplectic integrator that computes transit times and radial velocities of an exoplanetary system.

An example of the application based on the TTV technique can be found in Nesvorný et al. (2013), where the authors have predicted the presence of the planet KOI-142c in the system, which has been recently confirmed by Barros et al. (2014).

In Sect. 2 we introduce the TRADES program, the basic formulas, and methods to calculate radial velocities and transit times; in Sects. 35, we run TRADES on a synthetic 3-body system, Kepler-11, and Kepler-9, respectively. We summarize and discuss the results in the Sect. 6.

2. TRADES

We have developed a computer program (in Fortran 90 and openMP) for determining the possible physical and dynamical configurations of extra-solar planetary systems from observational data, known as TRADES, which stands for TRAnsits and Dynamics of Exoplanetary Systems. The program TRADES models the dynamics of multiple planet systems and reproduces the observed transit times (T0, or mid-transit times) and radial velocities (RVs). These T0s and RVs are computed during the integration of the planetary orbits. We have developed TRADES from zero because we want to avoid black-box programs, it would be easier to parallelize it with openMP, and include additional algorithms.

To solve the inverse problem, TRADES can be run in four different modes: 1) “grid” search; 2) Levenberg-Marquardt2 (LM) algorithm; 3) genetic algorithm (GA, we used the implementation named PIKAIA3, Charbonneau 1995); 4) and particle swarm optimization (PSO4, Tada 2007). In each mode, TRADES compares observed transit times (T0,obss) and radial velocities (RVobss) with the simulated ones (T0,sims and RVsims).

  • 1.

    In the grid search method, TRADES samples the orbital elementsof a perturbing body in a four-dimensional grid: the mass, M, the period, P (or the semi-major axis, a), the eccentricity, e, and the argument of the pericenter, ω. The grid parameters can be evenly sampled on a fixed grid by setting the number of steps, or the step size, or by a number of points chosen randomly within the parameter bounds. For any given set of values, the orbits are integrated, and the residuals between the observed and computed T0s and RVs are computed. For each combination of the parameters, the LM algorithm can be called and the best case is the one with the lowest residuals (lowest χ2). We have selected these four parameters for the grid search because they represent the minimal set of parameters required to model a coplanar system. In the future, we intend to add the possibility of making the grid search over all the set of parameters for each body.

  • 2.

    After an initial guess on the orbital parameters of the perturber, which could be provided by the previously described grid approach, the LM algorithm exploits the Levenberg-Marquardt minimization method to find the solution with the lowest residuals. The LM algorithm requires the analytic derivative of the model with respect to the parameters to be fitted. Since the T0s are determined by an iterative method and the radial velocities are computed using the numerical integrator, we cannot express these as analytic functions of fitting parameters. We have adopted the method described in Moré et al. (1980) to compute the Jacobian matrix, which is determined by a forward-difference approximation. The epsfcn parameter, which is the parameter that determines the first Jacobian matrix, is automatically selected in a logarithmic range from the machine precision up to 10-6; the best value is the one that returns the lowest χ2. This method has the advantage to be scale invariant, but it assumes that each parameter is varied by the same epsfcn value (e.g., a variation of 10% of the period has a different effect than a variation of the same percentage of the argument of pericenter).

  • 3.

    The GA mode searches for the best orbit by performing a genetic optimization (e.g. Holland 1975; Goldberg 1989), where the fitness parameter is set to the inverse of the χ2. This algorithm is inspired by natural selection which is the biological process of evolution. Each generation is a new population of “offspring” orbital configurations, that are the result of “parent” pairs of orbital configurations that are ranked following the fitness parameter. A drawback of the GA is the slowness of the algorithm, when compared to other optimizers. However, the GA should converge to a global solution (if it exists) after the appropriate number of iterations.

  • 4.

    The PSO is another optimization algorithm that searches for the global solution of the problem; this approach is inspired by the social behavior of bird flock and fish school (e.g., Kennedy & Eberhart 1995; Eberhart 2007). The fitness parameter used is the same as the GA, the inverse of the χ2. For each “particle”, the next step (or iteration) in the space of the fitted parameters is mainly given by the combination of three terms: random walk, best “particle” position (combination of parameters), and best “global” position (best orbital configuration of the all particles and all iterations).

The grid search is a good approach in case that we want to explore a limited subset of the parameter space or if we want to analyze the behavior of the system by varying some parameters, for example to test the effects of a growing mass for the perturbing planet. GA and PSO are good methods to be used in case of a wider space of parameters. The orbital solution determined with the GA or the PSO method is eventually refined with the LM mode.

For each mode, TRADES can perform a bootstrap analysis to calculate the interval of confidence of the best-fit parameter set. We generate a set of T0s and RVs from the fitted parameters, and we add a Gaussian noise having the calculated value (of T0s and RVs) as the mean and the corresponding measurement error as variance. We fit each new set of observables with the LM. We iterate the whole process thousands of times to analyze the distribution for each fitted parameter.

2.1. Orientation of the reference frame

For the transit time determination, the propagation of the trajectories of all planets in the system is performed in a reference frame with the Z-axis pointing to the observer, while the X-Y plane is the sky plane. At a given reference epoch, the Keplerian orbital elements of each planet are: period P (or semi-major axis a), inclination i, eccentricity e, argument of the pericenter ω, longitude of the ascending node Ω, and time of the passage at the pericenter τ (or the mean anomaly M). Given the orbital elements we first compute the initial radius r and velocity vectors in the orbital plane (e.g., see Murray & Dermott 2000):

where n = 2π/P is the mean motion, E is the eccentric anomaly obtained from the solution of the Kepler’s equation, M = EesinE, with the Newton-Raphson method (e.g., see Danby 1988; Murray & Dermott 2000; Murray & Correia 2011).

Then, we rotate the state vector by applying three consecutive rotation matrices, Rl(φ) (e.g., see Danby 1988; Murray & Dermott 2000; Murray & Correia 2011) where φ is the rotation angle and l is the rotation axis (where l is { 1,2,3 } for { x′,y′,z′ }). To rotate the initial state vector from the orbital plane to the observer reference frame, we have to use the transpose of the rotation matrix, with angles ω, i, and Ω. After this rotation, the X-Y plane is the sky plane with the Z-axis pointing to the observer, and we determine the initial state vector of each kth planet: (3)The same rotations have to be applied to the initial velocity vector. The inclinations are measured from the sky-plane. Indeed, a planet with inclination of has an orbit that lies on the sky-plane (X-Y plane), that is, it is seen face-on. The orbit of a planet with i = 90° is seen edge-on (it transits exactly through the center of the star), and it lies on the X-Z plane. From the initial state vector, TRADES integrates the astrocentric equation of motion (e.g. Murray & Dermott 2000; Fabrycky 2011) of planet k: (4)where M1 is the mass of the star and N the number of bodies; the first term is the direct gravitational force, and the second term is the indirect force due to mutual interaction of the planets. The orbits are computed with the Runge-Kutta-Cash-Karp integrator (RKCK, Cash & Karp 1990; Press et al. 1996). It is not a symplectic integrator5, and it is not well suited for long–term time integrations. Instead, it uses small and variable steps (it self adjusts the step-size to maintain the numerical precision during the computation of the orbits), is fast, and preserves the total energy and the total angular momentum during the time scales of our simulations.

2.2. Transit determination

We chose the change of sign of the X or Y coordinates between two consecutive steps of each planet trajectory as first condition of an eclipse. When this condition is met, following Fabrycky (2011,, we have to seek roots of the sky-projected separation rs,k ≡ (Xk,Yk) with the Newton-Raphson method by solving (5)then moving and iterating by the quantity, (6)In this way, we can determine with high precision the mid-transit time and the corresponding state vector (rmid, ) with an accuracy equal to the selected δt. We decided to set this accuracy in TRADES at the machine precision, which can be fine-tuned in the source of the code and defines the type of the chosen variables.

Then, we determine if we have four contact times, or just two (in the case of a grazing eclipse), or no transit, comparing the module of the sky-projected separation at the transit time, | rs,mid |, with the radius of the star, R, and of the planets, Rk, as in Fabrycky (2011). If the transit (or the occultation) does exist, we move about from the tmid backward () for first and second contact, and forward (+) for third and fourth contact. Then, we solve (7)and move of (8)until δt is less than the accuracy found (the same adopted in finding the transit time).

We based our approach on Fabrycky (2011), but we used a bisection-Newton-Raphson hybrid method, which is guaranteed to be bound near the solution, and we assume that the orbital elements of the bodies are almost constant around the center of the transit. Because of the latter assumption, we use F(ti,ti − 1) and G(ti,ti − 1) functions (called f(t,t0) and g(t,t0) in Danby 1988; Murray & Dermott 2000) to compute the planetary state vectors instead of the integrator while seeking the transit times: (9)where (10)with (11)and (12)where ti = ti − 1 + δti − 1 at the ith iterations; Ei and Ei − 1 are the eccentric anomalies at the ith and ith− 1 iterations. This allows the code to run faster than using the integrator, which has a lower number of function calls, to seek the transit and contact times.

The light coming from the star is delayed due to the motion of the star around the barycenter of the system (Irwin 1952), and so TRADES corrects for the light-time travel effect (LTE = , see Fabrycky 2011), each contact, and center transit time.

2.3. RV calculation and other constraints

For each observed RV (when available), TRADES integrates the orbits of the planets to the instant of the RV point and calculates the RV as the opposite of the z-component of the barycentric velocity of the star ( in the right unit of measurement). The observed RV is defined as RVobs = γ + rvobs, where γ is the motion of the barycenter of the system and rvobs is the reflex motion of the star induced by the planets. The program TRADES calculates γsim as the weighted mean of the difference Δrvj = RVj,obs − rvj,sim with j from one to the number of RVs. The final simulated RV is RVj,sim = γsim + rvj,sim. We are planning to implement the γsim fitting rather than the described weighted mean method.

Furthermore, we added some constraints on the orbit during the integration, setting a minimum and a maximum semi-major-axis for the system, amin and amax, respectively. The lower limit has been set equal to the star radius, while the maximum limit has been set equal to five times the largest semi-major axis of the system calculated from the periods of the planets. In the GA, we have used the largest period boundary. We use the definition of the Hill’s sphere to obtain minimum distance allowed between two planets (Murray & Dermott 2000). In this case, when these constraints are not be respected, the integration is stopped, and the χ2 returned is set to the maximum value allowed by the compiler, so the combination of the parameters is rejected.

3. Validation with simulated system

To validate TRADES, we simulated a synthetic system with two planets having known orbital parameters. We chose a star with a mass and radius equal to the Sun, a first planet named b with Jupiter mass and radius (MJup and RJup), and a second planet named c with mass and radius of Saturn (MSat and RSat). We assumed a co-planar system with inclination of 90° (perfectly edge-on). The input orbital elements of the system are summarized in Table 1.

Table 1

Parameters of the simulated system in Sect. 3.

We simulated the system with TRADES for 500 days. We computed all T0s for each body and we call these times the “true” transit times (T0,trues) of the system. Then, we created sets of synthetic transit times (T0,synths), such as (13)where s is a scaling factor varying from 0.01 to 1.5 (on twenty logarithmic steps) needed to simulate good to very bad measurement cases; N(0,1) is Gaussian noise with 0 as mean and 1 as variance. The P/ 3 factor is needed to scale the Gaussian noise in the right unit of time, and avoids confusion between transits and occultations at the same time. Furthermore, for each set of T0,synths, we selected a random number of transits (at least N/ 3 with N being the total number of transits of each planet) to simulate observed transits.

We fixed the orbital parameters of planet b, and we fitted M, P, e, ω, M, and Ω of planet c. We ran TRADES in LM mode for each scaling factor, and we calculated the difference of the parameters (Δ) as the determined parameters minus the input parameters. We repeated the simulation ten times (we calculated new Gaussian noise and the number of observed times every time). In Fig. 1, we plotted then mean and median of ten simulations for each s scaling factor value. The parameters of the system derived by TRADES depart from the “true” values only for extremely large measurement errors.

thumbnail Fig. 1

Mean (red-open circles with 1σ error bars) and median (blue-filled circles) variation (Δ) of fitted parameters of planet c for each value of the measurement errors on T0, which are here parametrized as s × P/ 3 (see text for details).

Open with DEXTER

To further test the robustness of the algorithm, we took the T0,trues (without added noise) and varied the initial semi-major axis of planet c from 0.19 AU to 0.21 AU with the TRADES grid+LM mode by fitting the same parameters of the previous test. The algorithm nicely converged to the values from which the synthetic data were generated, except for initial parameters that are too far from the right solution. This is due to the known limitation of the LM algorithm, which converges to close local minima from an initial set of parameters. Figure 2 shows the variation of the parameter differences (Δ) as function of the initial semi-major axis. This test shows how well TRADES recovers the parameters in the case of a bad guess of the initial parameters.

We measured the computational time required by TRADES, and we found that it can integrate (initial step size of 0.0001 day) a 3-body synthetic system for 3000 days writing the orbits, the Keplerian elements, and the constants of motion for each 0.1 days in about 2.3 seconds. An integration of 1000 days has been performed in less than 1 second and in less than half second for 500 days of integration time, but most of the time have been spent writing files. We want to stress that TRADES write these files only at the end of the simulations, so the real computation is faster than these estimates. The time required by TRADES to complete the grid search was about 51 min with 10 processors of an Intelr Xeonr CPU E5-2680 based workstation. For each combination of the initial parameters in the grid search, TRADES runs 10 times the LM to select the best value for the parameter epsfcn, that is needed to construct the initial Jacobian.

We tested the PSO+LM algorithm by fitting the same parameters of the grid search with limited boundaries except for the semi-major axis of planet c, for which we used the same limit of the grid search (ac = [ 0.19,0.21 ] AU). We ran this test four times with 200 particles for 2000 iterations, and TRADES always returned the right parameter values in less than about 1 hour and 40 min with 10 processors.

thumbnail Fig. 2

Variations (Δ) of the fitted parameters for different initial values of the semi-major axis of the planet c; each point corresponds to a different simulation. In this case, the input parameters (of planet c) are the parameters in Table 1 used to generate the exact transit times. The goodness of the fit has been color-coded so that good fits () have been plotted as blue points and bad fits () as open circles. The small gaps are due to the random sampling used to generate the grid in the semi-major axis a.

Open with DEXTER

4. Test case: Kepler-11 system

The system Kepler-11 (KOI-157, Lissauer et al. 2011a) has six transiting planets packed in less than 0.5 AU, making a complex and challenging case to be tested with TRADES.

From the spectroscopic analysis of HIRES high-resolution spectra, Lissauer et al. (2011a) derived the stellar parameters (effective temperature, surface gravity, metallicity, and projected stellar equatorial rotation) and determined the mass and the radius of Kepler-11 star to be 0.95 ± 0.10 M and 1.1 ± 0.1 R.

Table 2

Parameters of the Kepler-11 system. Epoch of reference: 2 455 190.0 (BJDUTC).

We first performed an analysis of the Kepler-11 system only on the data from the first three quarters of Kepler observations published by Lissauer et al. (2011a) and supplementary information (SI). We used the first circular model from the Lissauer et al. (2011a, SI) as an initial guess, which fixed the eccentricity and the longitude of the ascending node to zero for all the planets; hereafter we call this model Lis2011 (see first column in Table 2 for a summary of the orbital parameters). We used this model because the authors did not provide any information about the mean anomalies (or the time of the passage at the pericenter) for those planets. In this case, the argument of the pericenter, ω, is undetermined, so we fixed it to ω = 90° for each planet. We then calculated the initial mean anomaly, M0 at the reference epoch tepoch = 2 455 190.0 (BJDUTC6), setting the transit time (T0, Lissauer et al. 2011a, SI Table S4) as the time of passage at the pericenter: (14)where n is the mean motion of the planet.

Lissauer et al. (2011a) gave an upper limit of 300 M on the mass of the planet Kepler-11g, while they set it to zero in the three dynamical models of the supplementary information , and we followed the same approach. Figure 3 shows the orbits and RVs of the Kepler-11 planets, according to the Lis2011 model.

thumbnail Fig. 3

Orbits of the Kepler-11 system with initial parameters from Lis2011 model (circular model, see Table 2). The planet marker size is scaled with the mass of the planet. Top-right: “Sky Plane”, Kepler-11 system as seen from the Kepler satellite; we plotted only one orbit near a transit for each planet. Each circle is the position of a planet at a given integration time step. Bottom-right: projection of the system as seen face on. The big markers are the initial points of the integration. Bottom-left: RV model from the simulation.

Open with DEXTER

We fitted a linear ephemeris (Table 3) to the observed transit times of each planet for the first three quarters and computed the Observed-Calculated (O-C) diagrams, where O is T0,obss and C is the transit time calculated from the linear ephemeris.

We ran TRADES and fitted masses, periods, and mean anomalies of each planet. Hereafter, the orbital solution we have determined with TRADES is named with a short ID of the system (K11 for Kepler-11 and K9 for Kepler-9) and a Roman number (K11-I, K11-II, and so on). See solution K11-I in Table 2 for a summary of the parameters determined with TRADES (with 2σ confidence intervals from bootstrap analysis), which agree with those published by Lissauer et al. (2011a). For each bootstrap run, we run 1000 iterations to obtain the confidence intervals at the 97.72 percentile (2σ) of the distribution of each parameter. We calculated the residuals as the difference between the observed and simulated T0s. The residuals after the TRADES-LM fit are smaller than those with simulated transit times, that obtained with original parameters, and the final is around 1.25 for 88 degrees of freedom (dof, calculated as the difference between the number of data and the number of fitted parameters).

We also used the same initial conditions of solution K11-I, but this time we fitted the eccentricity and the argument of pericenter of all the planets. In this case, the LM did not move from the initial conditions even if it properly ended the simulation and returned reasonable errors. In the user guide of MINPACK (Moré et al. 1980), the user is warned to carefully analyze the case in which one has a null initial parameter. We set the initial eccentricities to a small but non-zero value of 0.0001. This small change was able to let the LM algorithm to properly return reasonable parameter values; see solution K11-II in Table 2 for a summary of the parameters.

The resulting masses of the solutions K11-I and K11-II all agree within 2σ with the discovery paper (Lissauer et al. 2011a) and with all the best-fit solutions determined by Migaszewski et al. (2012). In the latter work, the authors presented different sets of orbital parameters determined with an approach similar to ours (direct N-body simulation with genetic algorithm, Levenberg-Marquardt, and bootstrap), but they directly fit the flux of Kepler light curves (so-called dynamical photometric model) without fitting the transit times.

Table 3

Ephemeris fitted to the first three quarters of data of the Kepler-11 system.

4.1. Transit time analysis of the twelve quarters

Recently, Lissauer et al. (2013) analyzed the transit times covering fourteen quarters of Kepler data (in long and short cadence mode). A new independent extraction of T0s from the light curves made by Lissauer et al. (2013, hereafter we call the dynamical model from this work Lis2013, see column five of Table 2 led the authors to change the value of some parameters of the system. For example, they determined a mass of of the planet c that is lower than 15.82 ± 2.21 M published in the discovery paper (Lissauer et al. 2011a). Unfortunately the authors have not published the T0, so we used the data from Mazeh et al. (2013) that recently published the transit times for twelve quarters of the Kepler mission for 721 KOIs.

We used the linear ephemeris by Lissauer et al. (2013) to compute the O-Cs for the T0s from Mazeh et al. (2013). We found a remarkable mismatch with the O-Cs plotted in the paper by Lissauer et al. (2013). We stress that the T0s by Mazeh et al. (2013) are calculated with an automated algorithm. It would be advisable to more carefully analyze the light curves determining the T0s with higher precision, but it is not the purpose of this work. We analyzed the system with all transit times from Mazeh et al. (2013) without any selection, and they lead to unphysical results. We then decided to discard data with duration and depth of the transits that are 5σ away from the median values. This selection defines the sample of T0s for the first twelve quarters of Kepler-11 exoplanets on which we bases our next analysis.

We ran simulations with TRADES in grid+LM mode on twelve quarters with initial set of parameters as in K11-II; we fitted M, P, e, ω, and M of each planet (Ω = 0° fixed for all the planets). In particular, we varied the mass of planet g from 1 M to 100 M with a logarithmic step (ten simulations including the boundaries of 1 and 100 M) in the grid. We repeated this set of simulations for three different initial values of the eccentricity: in the first sample, we set the initial eccentricity of all planets to 0.001; in the second sample this is equal to 0.1; and in the third sample, we used a different value of the eccentricity for each planet, which closer to the Lissauer et al. (2013) ones: eb = 0.05, ec = 0.05, ed = 0.001, ee = 0.005, ef = 0.005, and eg = 0.1. With these simulations, we intended to test whether a forest of local minima are met during the search for the lowest χ2 (see Figs. 46). According to Figs. 4 and 5, this indeed seems to be the case significantly complicating the identification of the real minimum. The LM was not able to properly change the eccentricity of the planets that, which got stuck close to the initial value in the majority of the cases. Furthermore, when the initial eccentricity have been set to 0.1 the masses of planets d and f have decreased (Fig. 5) compared to those in the previous simulation (Fig. 4). Maybe, this could be an effect due to the particular sample of T0 we used, so it would be interesting to re-estimate the T0 from the light curves and re-analyze the system.

thumbnail Fig. 4

Masses (upper-panel) and eccentricities (lower-panel) for the Kepler-11 planets, calculated with TRADES in grid+LM mode (white-blue circle with blue error bars, see the legend on top of the upper plot) and 2σ confidence intervals from bootstrap analysis (red filled bars), with initial eccentricity of 0.001 for each planet. The blue-yellow circle (dark-green error bars) is the best simulation (number 10, , calculated mass, Mg,out, and input mass, Mg,out, are reported in the legend at the top of the plots). The different simulations (different initial mass of planet g, mg,in) have been plotted from left (first simulation) to right (eleventh simulation) for each planet. Masses and eccentricities by Lissauer et al. (2013) plotted as black lines (values on top of the plots) with the 2σ confidence intervals (light-gray filled bars). Red lines at 1 M (solid) and at 10 M (dashed).

Open with DEXTER

However, all the simulations have a final of about 2 or lower; the ninth simulation of the third set is our best solution (hereafter K11-III, see Table 2) with a (Fig. 6). The resulting O-C diagrams of the best simulation are shown in Figs. 7 (planets b, c, and d) and 8 (planet e, f, and g). All the masses and the eccentricities of solution K11-III agree well with the values found by Lissauer et al. (2013). Some of our simulations converged to parameter values, which are different from those proposed by Lissauer et al. (2013). Furthermore, some simulations show very narrow confidence intervals. This could be due both to the high complexity of the problem and to a strong selection effect: the distribution of the parameters in the bootstrap analysis are strongly bounded to the parameter values found by the LM algorithm.

In Table 4, we report a brief summary of the main differences of the characteristics of the analysis that led us to each solution for the Kepler-11 system.

thumbnail Fig. 5

Same as Fig. 4 but for simulations with initial eccentricities of 0.1.

Open with DEXTER

thumbnail Fig. 6

Same as Fig. 4 but for simulations with different initial eccentricities: eb = 0.05, ec = 0.05, ed = 0.001, ee = 0.005, ef = 0.005, and eg = 0.1. The best solution of this plot is the so-called K11-III solution; see Table 2 for the summary of the final parameters.

Open with DEXTER

thumbnail Fig. 7

O-C diagrams for planets b, c, and d (from top to bottom) of the Kepler-11 system; the black filled circles are the observed points fitted by a linear ephemeris, the blue open circles are the simulated points fitted by the same linear ephemeris of the observations. The simulated points are calculated from the solution K11-III (best simulation in Fig. 6) in Table 2. Residual plots, as the difference between observed and simulated central time (T0,obsT0,sim), are in the lower panel of each O-C plot. The unit of measurement of the left O-C y-axis is days (d) and minutes (m) are on the right. The N in the abscissa identifies the transit number respect to the reference transit time of the ephemeris of each body (second column of Table 3).

Open with DEXTER

thumbnail Fig. 8

Same as Fig. 7 but for Kepler-11 planets e, f, and g (top to bottom).

Open with DEXTER

Table 4

Main differences in the Kepler-11 analysis for each solution.

thumbnail Fig. 9

O-C diagrams (with residuals) from linear ephemeris for planet Kepler-9b (top panel) and Kepler-9c (middle panel) with the discovery paper’s parameters (see column two of Table 6); observations plotted as solid black circles, simulations plotted as open blue circles. The parameter N, in the abscissa, has the same meaning of Fig. 7. Bottom panel shows the RV observations as solid black circles, simulations at the same BJDUTC as open blue circles, and the dotted blue line is the RV model for the whole simulation.

Open with DEXTER

5. Test case: Kepler-9 system

Another ideal benchmark for testing TRADES is the multiple planet system Kepler-9 (KOI-377). The star Kepler-9 is a Solar-like G2 dwarf with a magnitude V = 13.9 (Holman et al. 2010), mass of 1.07 ± 0.05 M, and radius R = 1.02 ± 0.05 R (Torres et al. 2011). From the first three quarters of the Kepler data, Holman et al. (2010) identified two transiting Saturn-sized planet candidates (Kepler-9 b and c with radii of about ~0.8 RJup) near the 2:1 mean motion resonance (MMR). They detected an additional signal related to a third, smaller planet (KOI-377.03, estimated radius ~1.5 R), which is validated with BLENDER in Torres et al. (2011) but still unconfirmed. The last planet is not input in our simulations, given that there is no confirmation by spectroscopic follow-up so far (the expected RV semi-amplitude of about ~1.5 ms-1 would not increase the scatter in the RV data). Moreover, Holman et al. (2010) stated that the dynamical influences of the fourth body on other planets is undetectable on Kepler data (TTV amplitude of the order of ten seconds).

Table 5

Linear and quadratic ephemeris (in BJDUTC) fitted to data of Kepler-9 system.

From the analysis of the dP/ dt of the parabolic fit (quadratic ephemeris) for the T0s of each planet Holman et al. (2010) inferred the masses of Kepler-9b and Kepler-9c to be 0.252 ± 0.013 MJup, and 0.171 ± 0.013 MJup, respectively; they used the RV measurements from six spectra with the HIRES echelle spectrograph at Keck Observatory (Vogt et al. 1994) only to put a constraint on the masses. Holman et al. (2010) set an upper limit to the mass of the KOI-377.03 of about 7 M, but they could not fix the lower mass limit. The authors proposed 1 M for a volatile-rich planet with a hot extended atmosphere. Torres et al. (2011) could not determine a mass value for KOI-377.03 but estimated a radius of .

Table 6

Parameters of the Kepler-9 system at epoch tepoch = 2 455 088.212 BJDUTC.

We assumed the orbital parameters of the two planets at tepoch = 2 455 088.212 BJDUTC from Table S6 in Holman et al. (2010,supporting on-line material, SOM) and set (Holman 2012, priv. comm.; the value of reported in the supporting on-line material is inconsistent with the transit geometry). We simulated the system with TRADES without fitting any parameter, spanning the first three quarters of the Kepler observations. We fitted a linear ephemeris (see Table 5) to the observations and compared the resulting O-C diagrams with those from the simulations (Fig. 9). With the parameters from Holman et al. (2010), we obtained a simulated O-C for Kepler-9c, which is systematically offset from the observed data points by ~300 min (see Fig. 9, middle panel). In the bottom panel of the Fig. 9, we plot the RV model compared to the observations and the residuals.

We also reported the quadratic ephemeris for comparison with the discovery paper in Table 5. Our linear and quadratic ephemeris in Table 5 adopt the transit closest to the median epoch as time of reference for each body, while Holman et al. (2010) used the last transit time as reference for Kepler-9c.

To investigate whether this behavior is due to bugs in TRADES, we ran a second analysis with the MERCURY package (Chambers & Migliorini 1997). We simulated the same system with MERCURY and used the same technique as described in Sect. 2.2 to calculate the central time of the simulated transits. The maximum absolute difference between the mid-transit times from TRADES and MERCURY (with RADAU15 and Hybrid integrator) is ~0.16 s for an integration of 500 days. We did the calculation of the Keplerian orbital elements both for TRADES and MERCURY, and we verified the trend of the X coordinate (the coordinate used as alarm in case of eclipses for Kepler-9) of each planet as function of time (in a range of time around an observed transit): we did not find any difference or unexpected behavior between TRADES and MERCURY. These tests support our results showing that the problem is not in the integrator or in the subroutine used to calculate the transit times.

Then, we fitted M, P, e, ω, and M (mean anomaly) of both planets and Ω of planet c using the LM algorithm in TRADES. We found that the new values are consistent for all the fitted parameters with those by Holman et al. (2010, see column one and two of Table 6. Only one parameter, Pc, agrees with the discovery paper within 2σ. The small changes in the parameter values are enough to explain the O-C offset of planet c. The mean longitudes (λ = Ω + ω + M) of the two planets differs from the two solution only by few degrees, but this determines a small misalignment of the initial condition that could have a strong effect in MMR configuration. This simulation gives a χ2 ≈ 28.39, for 10 d.o.f., resulting in a . The results are summarized in Table 6 (solution K9-I) and in Fig. 10 the O-Cs and the RV diagrams are plotted (notations and colors as in Fig. 9).

thumbnail Fig. 10

Kepler-9 system: same plots as in Fig. 9 but with the parameters determined with TRADES-LM (K9-I of Table 6).

Open with DEXTER

5.1. Transit time analysis of the twelve quarters

As for the Kepler-11 system, we extended the analysis of Kepler-9 to the first twelve quarters of Kepler data using the transit times from Mazeh et al. (2013). We did not find any transit time to discard when using the same criteria used for Kepler-11. First of all, we extended the integration of the orbits of the planets from the solution K9-I to the twelve quarters (we did not fit any parameters in this simulation) and we compare the observed T0s and RVs with the simulated ones. In Fig. 11 it is clear that the simulation diverges quite soon from the observations. We run a simulation with the MERCURY package with same initial parameters of TRADES, compared the resulting O-C diagrams, and we found the same behavior. Furthermore, we calculated the transit time differences between TRADES and MERCURY, and we found that the maximum absolute difference is about 12 seconds, which is really smaller than the error bars of the T0s.

thumbnail Fig. 11

Kepler-9 system: O-C diagrams from the parameters obtained (solution K9-I in Table 6) with TRADES for the data from Holman et al. (2010) extended to the twelve quarters of Kepler. The simulations are compared with the T0s from Mazeh et al. (2013), while the RVs are from the discovery paper. The epoch of the transits (N in x-axis) are calculated from the linear ephemeris from Mazeh et al. (2013).

Open with DEXTER

thumbnail Fig. 12

Kepler-9 system: O-C diagrams from the fit with TRADES for the data from the twelve quarters of Kepler. T0s from Mazeh et al. (2013), while the RV are from the discovery paper. Initial parameters from the solution K9-I. for 62 d.o.f. Given the high value of the , we did not report the parameters of this solution.

Open with DEXTER

We considered the orbital solution K9-I in Table 6, and we ran a simulation on the T0s of Mazeh et al. (2013) for the same first three quarters of Holman et al. (2010). The six RV points are taken into account. The fitted orbital solution has all the parameters agree with the solution K9-I.

Then, we fit all the 12 quarters with initial condition from solution K9-I. The final is ~ 33 (for 62 d.o.f). The O-C diagrams (Fig. 12) are fitted better than those in Fig. 11, and the RV plot (bottom diagram in Fig 12) shows a lower amplitude.

To investigate the origin of this disagreement between observations and simulations when fitting 12 quarters (Figs. 11 and 12), we analyzed the T0s by Mazeh et al. (2013) with N simulations, with each one fitting three adjacent quarters of data (we called it “3 moving quarters”) and the six RV, in that we considered quarter 1 to 3, 2 to 4, and up to 10 to 12. We set the parameters from solution K9-I as the initial parameters of each simulation. We had good fits up to the simulation with quarters 6, 7, and 8; the following simulations showed increased (>700) that dropped to 44 only for the last three moving simulation (quarters 10, 11, and 12). In this analysis, we found that the bad fit starts when the solution K9-I diverges in Fig. 11. This could be an hint that the original solution determined by analyzing only the first three quarters of data is biased by the short time scale.

5.2. Dynamical analysis without RV points

Due to the high χ2 in Fig. 12, we re-analyzed the Kepler-9 system in a different way. We chose to run many simulations with GA+LM and PSO+LM on all 12 quarters with and without fitting the RV points. We set quite wide bounds on the parameters; in particular, we set the masses to be bound between 10-6MJup and 1 MJup and the eccentricities between 0 and 0.3. The best solution (K9-II) has been obtained with the TRADES mode PSO+LM without an RV fit. This solution has a of about 1.44 for 56 d.o.f. (summary of the final parameters in Table 6). The masses of solution K9-II are about 55% of the masses published by Holman et al. (2010). Furthermore, the eccentricities are smaller than the published ones (calculated from the SOM of Holman et al. 2010) and of the order of 0.06. These small values of the masses and the eccentricities imply a RV semi-amplitude of about 16.11 ms-1, which is smaller than the one from the solution K9-I (~28.91 ms-1) that is extended to all 12 quarters.

thumbnail Fig. 13

O-C diagrams (top and middle plot) for the solution K9-II in Table 6 from the fit of the T0s for the 12 quarters and neglecting the six RV points. Colors and markers as in Fig. 12 for 56 d.o.f. In the bottom panel, the RV model (blue dots) from the solution K9-II with over plotted RV observations (black dots with error-bars). We do not have the RV residuals from TRADES, because we have not fitted the RVs. The RV model has a lower RV semi-amplitude of about 12.80 ms-1 with respect to the RVs from the discovery paper.

Open with DEXTER

6. Summary

We have developed a program, TRADES, that simulates the dynamics of exoplanetary systems and that does a simultaneous fit of radial velocities and transit times data.

Analyzing a simulated planetary system, we have shown that TRADES can determine the parameters even from low precision data or for a rough guess of the initial orbital elements.

We validated TRADES by reproducing the packed exoplanetary system Kepler-11. The orbital parameters we determined agree with the values of the discovery paper and the recent analysis of 14 quarters by Lissauer et al. (2013). Furthermore, our analysis agrees with the results by Migaszewski et al. (2012) that used a similar approach. Our best simulation (K11-III) returned a value for the mass of the planet g of about 25 M, which agrees within the error bars and the confidence interval proposed by Lissauer et al. (2013). Furthermore, all our simulations showed a final , and the final mass of planet g agrees with the previous works.

We reproduced the Kepler-9 system (without KOI-377.03), and we found that the parameters from the SOM of the discovery paper (Holman et al. 2010) cannot properly reproduce the O-C diagram of Kepler-9c. We tested the orbits and O-C diagrams with an independent program, MERCURY, and we found the same result. A difference of a few degrees in λ for both planets is enough to explain the offset of about 300 min in the O-C diagram. We found the same results after analyzing the T0s by Mazeh et al. (2013) that cover the same quarters of Holman et al. (2010).

Extending the analysis, we found that the original solution is not compatible to the whole set of data from the 12 quarters by Mazeh et al. (2013). It shows a divergence of the simulated O-C compared to the observed one (Fig. 11). The solution K9-I, as obtained with the same temporal baseline of the discovery paper cannot explain the observed transit time for the whole set of T0s.

Only using the combination of a “quasi” global optimized search algorithm, such as genetic PIKAIA or particle swarm (PSO), with the LM algorithm, it has been possible to improve the fit on all 12 quarters but only if we neglect the six RV points. With this approach, we have found a new solution (K9-II) with a for 56 d.o.f. This solution lead to mass values that are about 55% of the mass values given in the discovery paper and smaller eccentricities. Due to these values, our RV model has a smaller semi-amplitude of about 12.80 ms-1 if compared to the observations.

We need to better study this system by, for example, obtaining more RV points because it can shed light on the issue of the different masses of the exoplanets calculated from RV data and from the TTV (for a similar case see Masuda et al. 2013). Follow-up transit observations with CHEOPS will extend the time coverage and are advisable.

In both Kepler systems we carried out the analysis of twelve quarters using the transit times calculated by Mazeh et al. (2013) using an automated algorithm. We point out that it would be advisable to analyze the light curves of the KOIs that show TTV signals and recompute the transit times of the planets in more detail.

It is known that the LM algorithm cannot return reliable errors (with physical meaning) in presence of correlated errors and complex parameter spaces. For some parameters, the bootstrap analysis returned small intervals of confidence. A possible explanation is that the parameter distributions in the bootstrap analysis are limited to values close to those found by the LM algorithm. This is probably due to a strong selection effect of the forest of minima in the χ2 space. Furthermore, for the Kepler-9 case, the measurement errors have tiny effect compared to the TTV signal that dominates the distribution of the parameters.

In the near future, we add a Monte-Carlo-Markov-chain algorithm (or an another Bayesian algorithm, e.g., multiNEST, Feroz et al. 2009) in TRADES to perform parameter estimation and model selection using a Bayesian approach. The idea is to provide the initial parameters to the Bayesian algorithm via the LM algorithm or use a grid search. Furthermore, we plan to include the transit duration in the fitting procedure to put more constraint on the parameter determination.

Recently, the European Space Agency adopted the space S-class mission CHEOPS (Broeg et al. 2013). The launch of this satellite is foreseen by December 2017. CHEOPS will characterize the structure of Neptune- and Earth-like exoplanets monitoring their transits. The CHEOPS targets will be stars hosting planets that are known via accurate Doppler and ground-based surveys. We plan to optimize TRADES for CHEOPS mission to dynamically characterize observed exoplanetary systems and to analyze possible TTV detections.

The work by Mazeh et al. (2013) provided a list of exoplanetary systems that show TTV signals. We are currently selecting a wide sample of these systems that would be suitable for an analysis with the program TRADES, such as Kepler-19. We plan to make a self-consistent analysis of the raw Kepler data, determine the transit times, and combine them with available RV data.

TRADES will be publicly release as soon as possible. Those interested can send an email to the first author of this work, and a copy of the code will be provided. At the moment the documentation is still in a preliminary form, but we are willing to provide all information needed.


1

http://exoplanet.eu/catalog, 2014 March 21st.

2

lmdif converted to Fortran 90 by Alan Miller (http://jblevins.org/mirror/amiller/) from MINPACK.

4

based on the public Fortran 90 code at http://www.nda.ac.jp/cc/users/tada/

5

A symplectic integrator has been designed to numerically solve the Hamilton’s equation by preserving the Poincaré invariants.

6

In the FITS header of Kepler data the time standard is reported as Barycentric Julian Day in Barycentric Dynamical Time (BJDTDB), but it is specified in the KSCI-19059 Subsect. 3.4 that the correct time is Barycentric Julian Day in Coordinated Universal Time (BJDUTC).

Acknowledgments

We thank the referee David Nesvorný for his careful reading and useful comments and suggestions. We acknowledge support from Italian Space Agency (ASI) regulated by Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013.

References

All Tables

Table 1

Parameters of the simulated system in Sect. 3.

Table 2

Parameters of the Kepler-11 system. Epoch of reference: 2 455 190.0 (BJDUTC).

Table 3

Ephemeris fitted to the first three quarters of data of the Kepler-11 system.

Table 4

Main differences in the Kepler-11 analysis for each solution.

Table 5

Linear and quadratic ephemeris (in BJDUTC) fitted to data of Kepler-9 system.

Table 6

Parameters of the Kepler-9 system at epoch tepoch = 2 455 088.212 BJDUTC.

All Figures

thumbnail Fig. 1

Mean (red-open circles with 1σ error bars) and median (blue-filled circles) variation (Δ) of fitted parameters of planet c for each value of the measurement errors on T0, which are here parametrized as s × P/ 3 (see text for details).

Open with DEXTER
In the text
thumbnail Fig. 2

Variations (Δ) of the fitted parameters for different initial values of the semi-major axis of the planet c; each point corresponds to a different simulation. In this case, the input parameters (of planet c) are the parameters in Table 1 used to generate the exact transit times. The goodness of the fit has been color-coded so that good fits () have been plotted as blue points and bad fits () as open circles. The small gaps are due to the random sampling used to generate the grid in the semi-major axis a.

Open with DEXTER
In the text
thumbnail Fig. 3

Orbits of the Kepler-11 system with initial parameters from Lis2011 model (circular model, see Table 2). The planet marker size is scaled with the mass of the planet. Top-right: “Sky Plane”, Kepler-11 system as seen from the Kepler satellite; we plotted only one orbit near a transit for each planet. Each circle is the position of a planet at a given integration time step. Bottom-right: projection of the system as seen face on. The big markers are the initial points of the integration. Bottom-left: RV model from the simulation.

Open with DEXTER
In the text
thumbnail Fig. 4

Masses (upper-panel) and eccentricities (lower-panel) for the Kepler-11 planets, calculated with TRADES in grid+LM mode (white-blue circle with blue error bars, see the legend on top of the upper plot) and 2σ confidence intervals from bootstrap analysis (red filled bars), with initial eccentricity of 0.001 for each planet. The blue-yellow circle (dark-green error bars) is the best simulation (number 10, , calculated mass, Mg,out, and input mass, Mg,out, are reported in the legend at the top of the plots). The different simulations (different initial mass of planet g, mg,in) have been plotted from left (first simulation) to right (eleventh simulation) for each planet. Masses and eccentricities by Lissauer et al. (2013) plotted as black lines (values on top of the plots) with the 2σ confidence intervals (light-gray filled bars). Red lines at 1 M (solid) and at 10 M (dashed).

Open with DEXTER
In the text
thumbnail Fig. 5

Same as Fig. 4 but for simulations with initial eccentricities of 0.1.

Open with DEXTER
In the text
thumbnail Fig. 6

Same as Fig. 4 but for simulations with different initial eccentricities: eb = 0.05, ec = 0.05, ed = 0.001, ee = 0.005, ef = 0.005, and eg = 0.1. The best solution of this plot is the so-called K11-III solution; see Table 2 for the summary of the final parameters.

Open with DEXTER
In the text
thumbnail Fig. 7

O-C diagrams for planets b, c, and d (from top to bottom) of the Kepler-11 system; the black filled circles are the observed points fitted by a linear ephemeris, the blue open circles are the simulated points fitted by the same linear ephemeris of the observations. The simulated points are calculated from the solution K11-III (best simulation in Fig. 6) in Table 2. Residual plots, as the difference between observed and simulated central time (T0,obsT0,sim), are in the lower panel of each O-C plot. The unit of measurement of the left O-C y-axis is days (d) and minutes (m) are on the right. The N in the abscissa identifies the transit number respect to the reference transit time of the ephemeris of each body (second column of Table 3).

Open with DEXTER
In the text
thumbnail Fig. 8

Same as Fig. 7 but for Kepler-11 planets e, f, and g (top to bottom).

Open with DEXTER
In the text
thumbnail Fig. 9

O-C diagrams (with residuals) from linear ephemeris for planet Kepler-9b (top panel) and Kepler-9c (middle panel) with the discovery paper’s parameters (see column two of Table 6); observations plotted as solid black circles, simulations plotted as open blue circles. The parameter N, in the abscissa, has the same meaning of Fig. 7. Bottom panel shows the RV observations as solid black circles, simulations at the same BJDUTC as open blue circles, and the dotted blue line is the RV model for the whole simulation.

Open with DEXTER
In the text
thumbnail Fig. 10

Kepler-9 system: same plots as in Fig. 9 but with the parameters determined with TRADES-LM (K9-I of Table 6).

Open with DEXTER
In the text
thumbnail Fig. 11

Kepler-9 system: O-C diagrams from the parameters obtained (solution K9-I in Table 6) with TRADES for the data from Holman et al. (2010) extended to the twelve quarters of Kepler. The simulations are compared with the T0s from Mazeh et al. (2013), while the RVs are from the discovery paper. The epoch of the transits (N in x-axis) are calculated from the linear ephemeris from Mazeh et al. (2013).

Open with DEXTER
In the text
thumbnail Fig. 12

Kepler-9 system: O-C diagrams from the fit with TRADES for the data from the twelve quarters of Kepler. T0s from Mazeh et al. (2013), while the RV are from the discovery paper. Initial parameters from the solution K9-I. for 62 d.o.f. Given the high value of the , we did not report the parameters of this solution.

Open with DEXTER
In the text
thumbnail Fig. 13

O-C diagrams (top and middle plot) for the solution K9-II in Table 6 from the fit of the T0s for the 12 quarters and neglecting the six RV points. Colors and markers as in Fig. 12 for 56 d.o.f. In the bottom panel, the RV model (blue dots) from the solution K9-II with over plotted RV observations (black dots with error-bars). We do not have the RV residuals from TRADES, because we have not fitted the RVs. The RV model has a lower RV semi-amplitude of about 12.80 ms-1 with respect to the RVs from the discovery paper.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.