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/00046361/201424080  
Published online  04 November 2014 
TRADES: A new software to derive orbital parameters from observed transit times and radial velocities
Revisiting Kepler11 and Kepler9
^{1}
Department of Physics and AstronomyUniversità degli Studi di
Padova,
Via Marzolo 8,
35131
Padova,
Italy
email:
luca.borsato.2@studenti.unipd.it
^{2}
INAF – Osservatorio Astronomico di Padova, vicolo
dell’Osservatorio 5, 35122
Padova,
Italy
Received: 28 April 2014
Accepted: 7 August 2014
Aims. With the purpose of determining the orbital parameters of exoplanetary systems from observational data, we have developed a software, named TRADES (TRAnsits and Dynamics of Exoplanetary Systems), to simultaneously fit observed radial velocities and transit times data.
Methods. We implemented a dynamical simulator for Nbody systems, which also fits the available data during the orbital integration and determines the best combination of the orbital parameters using grid search, χ^{2} minimization, genetic algorithms, particle swarm optimization, and bootstrap analysis.
Results. To validate TRADES, we tested the code on a synthetic threebody system and on two real systems discovered by the Kepler mission: Kepler9 and Kepler11. These systems are good benchmarks to test multiple exoplanet systems showing transit time variations (TTVs) due to the gravitational interaction among planets. We have found that orbital parameters of Kepler11 planets agree well with the values proposed in the discovery paper and with a a recent work from the same authors. We analyzed the first three quarters of Kepler9 system and found parameters in partial agreement with discovery paper. Analyzing transit times (T_{0}s), covering 12 quarters of Kepler data, that we have found a new bestfit solution. This solution outputs masses that are about 55% of the values proposed in the discovery paper; this leads to a reduced semiamplitude of the radial velocities of about 12.80 ms^{1}.
Key words: methods: numerical / celestial mechanics / stars: individual: Kepler11 / stars: individual: Kepler9
© ESO, 2014
1. Introduction
Now, more than 1779 planets^{1} 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 transitlike 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 (T_{0}s) of a planet may be not strictly periodic (see Agol et al. 2005; Holman & Murray 2005; MiraldaEscudé 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 T_{0}s are computed as a Fourier series. A drawback of this method is that it does not take into account for the socalled 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 groundbase facilities.
The method described in this paper is based on a direct numerical Nbody 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 KOI142c 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. 3–5, we run TRADES on a synthetic 3body system, Kepler11, and Kepler9, 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 extrasolar 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 (T_{0}, or midtransit times) and radial velocities (RVs). These T_{0}s and RVs are computed during the integration of the planetary orbits. We have developed TRADES from zero because we want to avoid blackbox 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) LevenbergMarquardt^{2} (LM) algorithm; 3) genetic algorithm (GA, we used the implementation named PIKAIA^{3}, Charbonneau 1995); 4) and particle swarm optimization (PSO^{4}, Tada 2007). In each mode, TRADES compares observed transit times (T_{0,obs}s) and radial velocities (RV_{obs}s) with the simulated ones (T_{0,sim}s and RV_{sim}s).

1.
In the grid search method, TRADES samples the orbital elementsof a perturbing body in a fourdimensional grid: the mass, M, the period, P (or the semimajor 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 T_{0}s 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 LevenbergMarquardt 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 T_{0}s 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 forwarddifference 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 bestfit parameter set. We generate a set of T_{0}s and RVs from the fitted parameters, and we add a Gaussian noise having the calculated value (of T_{0}s 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 Zaxis pointing to the observer, while the XY plane is the sky plane. At a given reference epoch, the Keplerian orbital elements of each planet are: period P (or semimajor 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 = E − esinE, with the NewtonRaphson method (e.g., see Danby 1988; Murray & Dermott 2000; Murray & Correia 2011).
Then, we rotate the state vector by applying three consecutive rotation matrices, R_{l}(φ) (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 XY plane is the sky plane with the Zaxis 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 skyplane. Indeed, a planet with inclination of 0° has an orbit that lies on the skyplane (XY plane), that is, it is seen faceon. The orbit of a planet with i = 90° is seen edgeon (it transits exactly through the center of the star), and it lies on the XZ 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 M_{1} 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 RungeKuttaCashKarp integrator (RKCK, Cash & Karp 1990; Press et al. 1996). It is not a symplectic integrator^{5}, and it is not well suited for long–term time integrations. Instead, it uses small and variable steps (it self adjusts the stepsize 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 skyprojected separation r_{s,k} ≡ (X_{k},Y_{k}) with the NewtonRaphson method by solving (5)then moving and iterating by the quantity, (6)In this way, we can determine with high precision the midtransit time and the corresponding state vector (r_{mid}, ) with an accuracy equal to the selected δt. We decided to set this accuracy in TRADES at the machine precision, which can be finetuned 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 skyprojected separation at the transit time,  r_{s,mid} , with the radius of the star, R_{⋆}, and of the planets, R_{k}, as in Fabrycky (2011). If the transit (or the occultation) does exist, we move about from the t_{mid} 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 bisectionNewtonRaphson 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(t_{i},t_{i − 1}) and G(t_{i},t_{i − 1}) functions (called f(t,t_{0}) and g(t,t_{0}) 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 t_{i} = t_{i − 1} + δt_{i − 1} at the ith iterations; E_{i} and E_{i − 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 lighttime 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 zcomponent of the barycentric velocity of the star ( in the right unit of measurement). The observed RV is defined as RV_{obs} = γ + rv_{obs}, where γ is the motion of the barycenter of the system and rv_{obs} is the reflex motion of the star induced by the planets. The program TRADES calculates γ_{sim} as the weighted mean of the difference Δrv_{j} = RV_{j,obs} − rv_{j,sim} with j from one to the number of RVs. The final simulated RV is RV_{j,sim} = γ_{sim} + rv_{j,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 semimajoraxis for the system, a_{min} and a_{max}, 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 semimajor 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 (M_{Jup} and R_{Jup}), and a second planet named c with mass and radius of Saturn (M_{Sat} and R_{Sat}). We assumed a coplanar system with inclination of 90° (perfectly edgeon). The input orbital elements of the system are summarized in Table 1.
We simulated the system with TRADES for 500 days. We computed all T_{0}s for each body and we call these times the “true” transit times (T_{0,true}s) of the system. Then, we created sets of synthetic transit times (T_{0,synth}s), 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 T_{0,synth}s, 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.
Fig. 1
Mean (redopen circles with 1σ error bars) and median (bluefilled circles) variation (Δ) of fitted parameters of planet c for each value of the measurement errors on T_{0}, which are here parametrized as s × P/ 3 (see text for details). 
To further test the robustness of the algorithm, we took the T_{0,true}s (without added noise) and varied the initial semimajor 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 semimajor 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 3body 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 Intel^{r} Xeon^{r} CPU E52680 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 semimajor axis of planet c, for which we used the same limit of the grid search (a_{c} = [ 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.
Fig. 2
Variations (Δ) of the fitted parameters for different initial values of the semimajor 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 colorcoded 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 semimajor axis a. 
4. Test case: Kepler11 system
The system Kepler11 (KOI157, 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 highresolution 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 Kepler11 star to be 0.95 ± 0.10 M_{⊙} and 1.1 ± 0.1 R_{⊙}.
Parameters of the Kepler11 system. Epoch of reference: 2 455 190.0 (BJD_{UTC}).
We first performed an analysis of the Kepler11 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, M_{0} at the reference epoch t_{epoch} = 2 455 190.0 (BJD_{UTC}^{6}), setting the transit time (T_{0}, 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 Kepler11g, 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 Kepler11 planets, according to the Lis2011 model.
Fig. 3
Orbits of the Kepler11 system with initial parameters from Lis2011 model (circular model, see Table 2). The planet marker size is scaled with the mass of the planet. Topright: “Sky Plane”, Kepler11 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. Bottomright: projection of the system as seen face on. The big markers are the initial points of the integration. Bottomleft: RV model from the simulation. 
We fitted a linear ephemeris (Table 3) to the observed transit times of each planet for the first three quarters and computed the ObservedCalculated (OC) diagrams, where O is T_{0,obs}s 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 Kepler11 and K9 for Kepler9) and a Roman number (K11I, K11II, and so on). See solution K11I 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 T_{0}s. The residuals after the TRADESLM 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 K11I, 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 nonzero value of 0.0001. This small change was able to let the LM algorithm to properly return reasonable parameter values; see solution K11II in Table 2 for a summary of the parameters.
The resulting masses of the solutions K11I and K11II all agree within 2σ with the discovery paper (Lissauer et al. 2011a) and with all the bestfit 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 Nbody simulation with genetic algorithm, LevenbergMarquardt, and bootstrap), but they directly fit the flux of Kepler light curves (socalled dynamical photometric model) without fitting the transit times.
Ephemeris fitted to the first three quarters of data of the Kepler11 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 T_{0}s 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 T_{0}, 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 OCs for the T_{0}s from Mazeh et al. (2013). We found a remarkable mismatch with the OCs plotted in the paper by Lissauer et al. (2013). We stress that the T_{0}s by Mazeh et al. (2013) are calculated with an automated algorithm. It would be advisable to more carefully analyze the light curves determining the T_{0}s 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 T_{0}s for the first twelve quarters of Kepler11 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 K11II; 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: e_{b} = 0.05, e_{c} = 0.05, e_{d} = 0.001, e_{e} = 0.005, e_{f} = 0.005, and e_{g} = 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. 4–6). 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 T_{0} we used, so it would be interesting to reestimate the T_{0} from the light curves and reanalyze the system.
Fig. 4
Masses (upperpanel) and eccentricities (lowerpanel) for the Kepler11 planets, calculated with TRADES in grid+LM mode (whiteblue 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 blueyellow circle (darkgreen error bars) is the best simulation (number 10, , calculated mass, M_{g,out}, and input mass, M_{g,out}, are reported in the legend at the top of the plots). The different simulations (different initial mass of planet g, m_{g,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 (lightgray filled bars). Red lines at 1 M_{⊕} (solid) and at 10 M_{⊕} (dashed). 
However, all the simulations have a final of about 2 or lower; the ninth simulation of the third set is our best solution (hereafter K11III, see Table 2) with a (Fig. 6). The resulting OC 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 K11III 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 Kepler11 system.
Fig. 6
Same as Fig. 4 but for simulations with different initial eccentricities: e_{b} = 0.05, e_{c} = 0.05, e_{d} = 0.001, e_{e} = 0.005, e_{f} = 0.005, and e_{g} = 0.1. The best solution of this plot is the socalled K11III solution; see Table 2 for the summary of the final parameters. 
Fig. 7
OC diagrams for planets b, c, and d (from top to bottom) of the Kepler11 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 K11III (best simulation in Fig. 6) in Table 2. Residual plots, as the difference between observed and simulated central time (T_{0,obs} − T_{0,sim}), are in the lower panel of each OC plot. The unit of measurement of the left OC yaxis 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). 
Main differences in the Kepler11 analysis for each solution.
Fig. 9
OC diagrams (with residuals) from linear ephemeris for planet Kepler9b (top panel) and Kepler9c (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 BJD_{UTC} as open blue circles, and the dotted blue line is the RV model for the whole simulation. 
5. Test case: Kepler9 system
Another ideal benchmark for testing TRADES is the multiple planet system Kepler9 (KOI377). The star Kepler9 is a Solarlike 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 Saturnsized planet candidates (Kepler9 b and c with radii of about ~0.8 R_{Jup}) near the 2:1 mean motion resonance (MMR). They detected an additional signal related to a third, smaller planet (KOI377.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 followup so far (the expected RV semiamplitude 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).
Linear and quadratic ephemeris (in BJD_{UTC}) fitted to data of Kepler9 system.
From the analysis of the dP/ dt of the parabolic fit (quadratic ephemeris) for the T_{0}s of each planet Holman et al. (2010) inferred the masses of Kepler9b and Kepler9c to be 0.252 ± 0.013 M_{Jup}, and 0.171 ± 0.013 M_{Jup}, 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 KOI377.03 of about 7 M_{⊕}, but they could not fix the lower mass limit. The authors proposed 1 M_{⊕} for a volatilerich planet with a hot extended atmosphere. Torres et al. (2011) could not determine a mass value for KOI377.03 but estimated a radius of .
Parameters of the Kepler9 system at epoch t_{epoch} = 2 455 088.212 BJD_{UTC}.
We assumed the orbital parameters of the two planets at t_{epoch} = 2 455 088.212 BJD_{UTC} from Table S6 in Holman et al. (2010,supporting online material, SOM) and set (Holman 2012, priv. comm.; the value of reported in the supporting online 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 OC diagrams with those from the simulations (Fig. 9). With the parameters from Holman et al. (2010), we obtained a simulated OC for Kepler9c, 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 Kepler9c.
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 midtransit 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 Kepler9) 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, P_{c}, agrees with the discovery paper within 2σ. The small changes in the parameter values are enough to explain the OC 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 K9I) and in Fig. 10 the OCs and the RV diagrams are plotted (notations and colors as in Fig. 9).
Fig. 10
Kepler9 system: same plots as in Fig. 9 but with the parameters determined with TRADESLM (K9I of Table 6). 
5.1. Transit time analysis of the twelve quarters
As for the Kepler11 system, we extended the analysis of Kepler9 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 Kepler11. First of all, we extended the integration of the orbits of the planets from the solution K9I to the twelve quarters (we did not fit any parameters in this simulation) and we compare the observed T_{0}s 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 OC 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 T_{0}s.
Fig. 11
Kepler9 system: OC diagrams from the parameters obtained (solution K9I 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 T_{0}s from Mazeh et al. (2013), while the RVs are from the discovery paper. The epoch of the transits (N in xaxis) are calculated from the linear ephemeris from Mazeh et al. (2013). 
Fig. 12
Kepler9 system: OC diagrams from the fit with TRADES for the data from the twelve quarters of Kepler. T_{0}s from Mazeh et al. (2013), while the RV are from the discovery paper. Initial parameters from the solution K9I. for 62 d.o.f. Given the high value of the , we did not report the parameters of this solution. 
We considered the orbital solution K9I in Table 6, and we ran a simulation on the T_{0}s 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 K9I.
Then, we fit all the 12 quarters with initial condition from solution K9I. The final is ~ 33 (for 62 d.o.f). The OC 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 T_{0}s 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 K9I 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 K9I 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 reanalyzed the Kepler9 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^{6}M_{Jup} and 1 M_{Jup} and the eccentricities between 0 and 0.3. The best solution (K9II) 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 K9II 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 semiamplitude of about 16.11 ms^{1}, which is smaller than the one from the solution K9I (~28.91 ms^{1}) that is extended to all 12 quarters.
Fig. 13
OC diagrams (top and middle plot) for the solution K9II in Table 6 from the fit of the T_{0}s 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 K9II with over plotted RV observations (black dots with errorbars). We do not have the RV residuals from TRADES, because we have not fitted the RVs. The RV model has a lower RV semiamplitude of about 12.80 ms^{1} with respect to the RVs from the discovery paper. 
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 Kepler11. 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 (K11III) 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 Kepler9 system (without KOI377.03), and we found that the parameters from the SOM of the discovery paper (Holman et al. 2010) cannot properly reproduce the OC diagram of Kepler9c. We tested the orbits and OC 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 OC diagram. We found the same results after analyzing the T_{0}s 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 OC compared to the observed one (Fig. 11). The solution K9I, as obtained with the same temporal baseline of the discovery paper cannot explain the observed transit time for the whole set of T_{0}s.
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 (K9II) 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 semiamplitude 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). Followup 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 Kepler9 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 MonteCarloMarkovchain 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 Sclass mission CHEOPS (Broeg et al. 2013). The launch of this satellite is foreseen by December 2017. CHEOPS will characterize the structure of Neptune and Earthlike exoplanets monitoring their transits. The CHEOPS targets will be stars hosting planets that are known via accurate Doppler and groundbased 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 Kepler19. We plan to make a selfconsistent 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.
http://exoplanet.eu/catalog, 2014 March 21st.
lmdif converted to Fortran 90 by Alan Miller (http://jblevins.org/mirror/amiller/) from MINPACK.
PIKAIA (http://www.hao.ucar.edu/modeling/pikaia/pikaia.php) converted to Fortran 90 by Alan Miller.
based on the public Fortran 90 code at http://www.nda.ac.jp/cc/users/tada/
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 ASIINAF n. 2013016R.0 del 9 luglio 2013.
References
 Agol, E., & Steffen, J. H. 2007, MNRAS, 374, 941 [NASA ADS] [CrossRef] [Google Scholar]
 Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Barros, S. C. C., Diaz, R. F., Santerne, A., et al. 2014, A&A, 561, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, EPJ Web Conf., 47, 3005 [Google Scholar]
 Cash, J. R., & Karp, A. H. 1990, ACM Trans. Math. Softw., 16, 201 [CrossRef] [Google Scholar]
 Chambers, J. E., & Migliorini, F. 1997, in AAS/Division for Planetary Sciences Meeting Abstracts #29 (Richmond: WillamannBell), BAAS, 29, 1024 [Google Scholar]
 Charbonneau, P. 1995, ApJS, 101, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Danby, J. M. A. 1988, Fundamentals of celestial mechanics (Richmond, Va., USA: WillmannBell) [Google Scholar]
 Deck, K. M., Agol, E., Holman, M. J., & Nesvorný, D. 2014, ApJ, 787, 132 [Google Scholar]
 Deprit, A. 1969, Celest. Mech., 1, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Eberhart, R. C. 2007, Computational Intelligence: Concepts to Implementations (San Francisco: Morgan Kaufmann Publishers Inc.) [Google Scholar]
 Fabrycky, D. C. 2011, NonKeplerian Dynamics of Exoplanets (University of Arizona Press), ed. S. Seager, 217 [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., Rowe, J. F., Fabrycky, D. C., et al. 2011, ApJS, 197, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Goldberg, D. E. 1989, Genetic Algorithms in Search, Optimization and Machine Learning, 1st edn. (Boston: AddisonWesley Longman Publishing Co., Inc.) [Google Scholar]
 Holland, J. 1975, Adaptation in Natural and Artificial Systems (The University of Michigan Press) [Google Scholar]
 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Holman, M. J., Fabrycky, D. C., Ragozzine, D., et al. 2010, Science, 330, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hori, G. 1966, PASJ, 18, 287 [NASA ADS] [Google Scholar]
 Irwin, J. B. 1952, ApJ, 116, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Kennedy, J., & Eberhart, R. 1995, in Neural Networks, 1995, Proc. IEEE Inter. Conf., 4, 1942 [Google Scholar]
 Latham, D. W., Rowe, J. F., Quinn, S. N., et al. 2011, ApJ, 732, L24 [Google Scholar]
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011a, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011b, ApJS, 197, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., JontofHutter, D., Rowe, J. F., et al. 2013, ApJ, 770, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Masuda, K., Hirano, T., Taruya, A., Nagasawa, M., & Suto, Y. 2013, ApJ, 778, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Mazeh, T., Nachmani, G., Holczer, T., et al. 2013, ApJS, 208, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Migaszewski, C., Słonina, M., & Goździewski, K. 2012, MNRAS, 427, 770 [NASA ADS] [CrossRef] [Google Scholar]
 MiraldaEscudé, J. 2002, ApJ, 564, 1019 [NASA ADS] [CrossRef] [Google Scholar]
 Moré, J. J., Garbow, B. S., & Hillstrom, K. E. 1980, User guide for MINPACK1, Tech. Rep. ANL8074, Argonne Nat. Lab., Argonne, IL [Google Scholar]
 Murray, C. D., & Correia, A. C. M. 2011, in Keplerian Orbits and Dynamics of Exoplanets (Tucson: University of Arizona Press), ed. S. Piper, 15 [Google Scholar]
 Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Nascimbeni, V., Piotto, G., Bedin, L. R., & Damasso, M. 2011, A&A, 527, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nesvorný, D. 2009, ApJ, 701, 1116 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., & Morbidelli, A. 2008, ApJ, 688, 636 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., Kipping, D., Terrell, D., et al. 2013, ApJ, 777, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1996, Numerical recipes in Fortran 90: the art of parallel scientific computing, 2nd edn. (New York: Cambridge University Press) [Google Scholar]
 Steffen, J. H., & Agol, E. 2005, MNRAS, 364, L96 [Google Scholar]
 Tada, T. 2007, Journal of Japan Society of Hydrology & Water Resources, 20, 450 [Google Scholar]
 Torres, G., Fressin, F., Batalha, N. M., et al. 2011, ApJ, 727, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in SPIE Conf. Ser. 2198, eds. D. L. Crawford, & E. R. Craine, 362 [Google Scholar]
All Tables
Parameters of the Kepler11 system. Epoch of reference: 2 455 190.0 (BJD_{UTC}).
Linear and quadratic ephemeris (in BJD_{UTC}) fitted to data of Kepler9 system.
All Figures
Fig. 1
Mean (redopen circles with 1σ error bars) and median (bluefilled circles) variation (Δ) of fitted parameters of planet c for each value of the measurement errors on T_{0}, which are here parametrized as s × P/ 3 (see text for details). 

In the text 
Fig. 2
Variations (Δ) of the fitted parameters for different initial values of the semimajor 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 colorcoded 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 semimajor axis a. 

In the text 
Fig. 3
Orbits of the Kepler11 system with initial parameters from Lis2011 model (circular model, see Table 2). The planet marker size is scaled with the mass of the planet. Topright: “Sky Plane”, Kepler11 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. Bottomright: projection of the system as seen face on. The big markers are the initial points of the integration. Bottomleft: RV model from the simulation. 

In the text 
Fig. 4
Masses (upperpanel) and eccentricities (lowerpanel) for the Kepler11 planets, calculated with TRADES in grid+LM mode (whiteblue 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 blueyellow circle (darkgreen error bars) is the best simulation (number 10, , calculated mass, M_{g,out}, and input mass, M_{g,out}, are reported in the legend at the top of the plots). The different simulations (different initial mass of planet g, m_{g,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 (lightgray filled bars). Red lines at 1 M_{⊕} (solid) and at 10 M_{⊕} (dashed). 

In the text 
Fig. 5
Same as Fig. 4 but for simulations with initial eccentricities of 0.1. 

In the text 
Fig. 6
Same as Fig. 4 but for simulations with different initial eccentricities: e_{b} = 0.05, e_{c} = 0.05, e_{d} = 0.001, e_{e} = 0.005, e_{f} = 0.005, and e_{g} = 0.1. The best solution of this plot is the socalled K11III solution; see Table 2 for the summary of the final parameters. 

In the text 
Fig. 7
OC diagrams for planets b, c, and d (from top to bottom) of the Kepler11 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 K11III (best simulation in Fig. 6) in Table 2. Residual plots, as the difference between observed and simulated central time (T_{0,obs} − T_{0,sim}), are in the lower panel of each OC plot. The unit of measurement of the left OC yaxis 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). 

In the text 
Fig. 8
Same as Fig. 7 but for Kepler11 planets e, f, and g (top to bottom). 

In the text 
Fig. 9
OC diagrams (with residuals) from linear ephemeris for planet Kepler9b (top panel) and Kepler9c (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 BJD_{UTC} as open blue circles, and the dotted blue line is the RV model for the whole simulation. 

In the text 
Fig. 10
Kepler9 system: same plots as in Fig. 9 but with the parameters determined with TRADESLM (K9I of Table 6). 

In the text 
Fig. 11
Kepler9 system: OC diagrams from the parameters obtained (solution K9I 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 T_{0}s from Mazeh et al. (2013), while the RVs are from the discovery paper. The epoch of the transits (N in xaxis) are calculated from the linear ephemeris from Mazeh et al. (2013). 

In the text 
Fig. 12
Kepler9 system: OC diagrams from the fit with TRADES for the data from the twelve quarters of Kepler. T_{0}s from Mazeh et al. (2013), while the RV are from the discovery paper. Initial parameters from the solution K9I. for 62 d.o.f. Given the high value of the , we did not report the parameters of this solution. 

In the text 
Fig. 13
OC diagrams (top and middle plot) for the solution K9II in Table 6 from the fit of the T_{0}s 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 K9II with over plotted RV observations (black dots with errorbars). We do not have the RV residuals from TRADES, because we have not fitted the RVs. The RV model has a lower RV semiamplitude of about 12.80 ms^{1} with respect to the RVs from the discovery paper. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.