Nauyaca: a New Tool to Determine Planetary Masses and Orbital Elements through Transit Timing Analysis

Nauyaca incorporates several modules with techniques adapted specifically for the TTVs inversion problem. We describe the functionalities of these modules below.

2.1.1. Setup

Here, we describe the preparation of the data before running the simulations. In the context of object-oriented programming, we define a Planetary System object by specifying the stellar mass and radius and the harbored planets. Transit times fitting will be operated over Planetary Systems. Then, we define the planets by establishing the information about their transit times ephemeris and the allowed parameter space. The transit ephemeris per planet should include the epoch number, time of transit, and lower and upper timing uncertainties. In the simulations, the transit epochs (integer transit numbers) of all of the planets are counted starting from 0 after t₀. Thus, the user must be aware that the epoch numbers are properly labeled and referenced to the same reference time t₀.

The time span of the simulations is automatically chosen to encompass the full time of observations in the ephemeris data. If no time step is defined for the simulations, it is automatically set to be 30 steps per orbit of the internal planet (∼3.33% of the internal planet period). With this time step Deck et al. (2014) demonstrated the effectiveness of determining transit times with an accuracy within 10 s for a ∼22.3 days period planet. A time step <5% of the internal planet period is recommendable to reach that accuracy.

Regarding the parameter space, the tool requires specific boundaries for each parameter in order to perform an effective sampling, avoiding nonphysical regions or redundant information. Each planet in the Planetary System must define its own parameter space. If there is no constraint for any individual or none of the parameters in Θ_j , a set of default boundaries is established. By default, masses are restricted to be in the range from 1 lunar mass to 80 Jupiter masses. However, when information about the stellar mass is supplied, the upper planetary mass limit is recalculated to be at most 1% of the stellar mass. This is done to keep valid the N-body Hamiltonian internally solved by TTVFast. In Figure 1, we show the currently measured planet masses M_p or ${M}_{p}\sin i$ from Exoplanet Archive ³ . It shows the planet mass limit corresponding to 1% of the stellar mass. We find that ∼95% of the currently known measured planet masses have mass ratios with the parent star lower than 1%, and thus the selected cutoff of 1% is statistically valid for most of the currently known planets. The default period space is limited to be between 0.1 and 1000 days. The eccentricity is limited between 1e-6 and 0.9, where the lower limit is different from exactly 0 to avoid undefined orbital angles. Inclination is defined between 0° and 180°. In the case of periodic boundaries for the argument of periastron, ω, and mean anomaly, M, fixing boundaries in the range 0°–360° could lead to an improper sampling process when the solution is near to these borders. This problem is solved internally by parameterizing the angles by means of ω + M and ω – M. The boundaries in the parameterized space encompass 720°. At the end of the sampling process, the results in the parameterized space are mapped to be between 0° and 360° for the individual ω and M. It is an internal process, so the user only defines these angles between 0° and 360°. The ascending node also must be defined between 0° and 360°. These default ranges can be modified by the user to better constrain any parameter or to keep it fixed. In the case of assuming constant parameters, these are not part of the sampling process.

Zoom In Zoom Out Reset image size

After setting the parameter space, Nauyaca normalizes the boundaries of all of the planets to dimensionless boundaries between 0 and 1. The whole sampling process (both for optimization and the MCMC) is performed with the normalized boundaries and is returned to physical values at the end of the runs. This normalization is done to remove the differences in orders of magnitude between parameters, which enhances the sampling performance.

2.1.2. Optimization Module

It combines minimization algorithms ordered sequentially to reach solutions Θ_j that best explain the transit times. These results can be used as initial guesses for the MCMC. We tested many algorithms and their calling order and found that the sequence Differential Evolution (DE; Storn & Price 1997), Powell (PW; Powell 1964), and Nelder−Mead (NM; Gao & Han 2012) progressively minimizes the total χ² (Equation (2)). This sequence is not arbitrary but reflects the nature of exploration of each algorithm. First, DE is capable of exploring a large parameter space without the necessity of requiring an initial guess. The algorithm is fully stochastic and no information about the smoothness of the space is required. The only mandatory information needed is the parameter space to explore, which has been set in the Setup module. A drawback of this algorithm is the slow convergence rate, and thus the outputs of this method could correspond to unconverged solutions. Even so, these solutions are better than any random proposal in the parameter space. Therefore, these are used as an initial guess to run PW, which is suitable to perform a minimum searching, assuming that the space around the starting point is continuous although complex. Finally, NM takes the solution previously found by PW and performs a downhill simplex method assuming that locally the parameter space is smooth and unimodal.

We show in Figure 2 an example of the progressive χ² minimization for the transit times fitting of a two-planet system with 190 transit times in total. For this example we performed 320 realizations. Background lines connect solutions of the same run, where a gradual descending of the χ² throughout the sequence is observed. The box plots show the intervals of the χ² achieved with the different algorithms, where typically the χ² is reduced by around 4–5 orders of magnitude between the first and last methods.

Zoom In Zoom Out Reset image size

Performing multiple realizations of this process enhances the chance of finding a global minimum but also provides us with a global view of the most probable regions in the parameter space (this will be addressed in detail in Section 4.3). It is a fast way of finding a starting point region with valuable information that can help to delimit the searching radius that the MCMC routine will explore.

2.1.3. MCMC Module

In order to set up planetary systems as realistically as possible, we selected masses and radii for the stellar hosts and masses, radii, and orbital periods for planets from the available data of confirmed planetary systems from the Exoplanet Archive. ⁴ This was done to create synthetic planetary systems with stellar and planetary parameters based on observations. The catalog includes different planetary multiplicities (number of planets in the system), ranging from two up to five planets. Below we describe the procedure followed to create the catalog.

First, we selected systems with planets discovered by the transit method. Second, we selected systems with reported stellar masses and radii. For those with unavailable data but with reported effective temperature, surface gravity, and metallicity, we derived the stellar mass or radius from the work of Torres et al. (2010). Third, we selected those systems for which all of their planets have reported masses, radii, and periods (systems whose planets have these parameters unreported were discarded). This procedure reduces significantly the number of planets, dropping to 2.5% (105 planets) of the original observed planets. In order to increase the number of synthetic planetary systems we applied an over-sampling technique (SMOTE; Chawla et al. 2011), which makes new samples by interpolation of the k-nearest neighbors. Thus, we tripled the number of planetary systems and their parameters, namely, stellar mass and radius, and planetary masses, radii, and periods. The remaining five orbital elements to complete a unique planetary configuration were made from random uniform distributions in the intervals: ecc [0.05, 0.2], inc [895, 905], ω [0°, 360°], M [0°, 360°], and Ω [88°, 92°]. We restricted the intervals of the inclination inc and ascending node Ω to get near coplanar and prograde orbits. These restrictions in the construction of the catalog are independent of the tool test. In practice, we can expand the boundaries of any parameter to be explored, as long as they have a physical meaning (see Section 4.1).

Once the complete set of planetary parameters is established, an N-body integration was performed over 10⁶ orbits of the internal planet using REBOUND (Rein & Liu 2012; Rein & Tamayo 2015). We used the chaos indicator MEGNO (Cincotta et al. 2003; Maffione et al. 2011; Rein & Tamayo 2015) to test the dynamical stability of the proposed planetary system. A MEGNO value of around ∼2 indicates a quasi-periodic motion (regular stable orbits). If the set of parameters of the proposed planetary system turned out to be stable (without a planetary ejection or with 1.7 < MEGNO < 2.3), then we appended the final state of the N-body simulation as an entry in the catalog. It should be pointed out that the orbital elements of these entries correspond to the osculating elements at the end of the stable N-body runs and they are not the same as the random values taken from the intervals of parameters indicated above. Finally, we used these entries as initial conditions to run TTVs simulations encompassing 130 transits of the internal planet using TTVFast. We selected that number to mimic the number of transits for planets with periods <10 days, during ∼3 yr of observation. We used a fixed number of transits rather than a fixed time span since it is more relevant for the inversion problem (Nesvorný & Morbidelli 2008). As a result, we got 130 transit times ephemeris for all of the internal planets, and a few less for the remaining planets according to their periods and orbital configurations.

We assigned synthetic errors to these transit times according to the analytical approximation of the timing precision (Agol & Fabrycky 2018),

$$\begin{eqnarray}&&{\sigma }_{t}=\displaystyle \frac{1}{\sqrt{2}}\ {\tau }^{1/2}{\dot{N}}^{-1/2}{\delta }^{-1}\end{eqnarray} \tag{ 3 }$$

where τ is the approximate duration of the transit ingress $\tau \approx 2.2\mathrm{minutes}\left({R}_{p}/{R}_{\oplus }\right){\left({M}_{* }/{M}_{\odot }\right)}^{-1/3}{\left(P/10\,{day}\right)}^{1/3}$, which is a function of planetary radius R_p , orbital period P, and stellar mass M_*, $\delta ={({R}_{p}/{R}_{* })}^{2}$ is the depth of the transit, and $\dot{N}$ is related to the Poisson noise due to the count rate of the star. We assumed $\dot{N}$ to be constant and equal to 1 × 10⁷ e⁻ minutes⁻¹ to mimic the value in the Kepler CCDs for a star of magnitude ∼12 in the Kepler band (Gilliland et al. 2011). Uncertainties from Equation (3) take planetary parameters derived from the Kepler photometry, and, therefore, our study is centered on the characterization of Kepler-like transits. Finally, we added white noise to the transit time models assuming a normal distribution with the mean equal to the true transit times and with standard deviation equal to the typical uncertainties given by Equation (3).

The full mock catalog is composed of twelve two-planet systems, eight three-planet systems, eight four-planet systems, and four five-planet systems, giving a total of 32 systems harboring 100 planets in total. We remark that all of the planets in the catalog transit, and thus our current study does not include the characterization of nontransiting planets. Tables in the Appendix include the parameters of the simulated mock catalog. Figure 3 shows the distributions of the parameters in the catalog according to the planetary multiplicity.

Zoom In Zoom Out Reset image size

The mock catalog of planetary parameters produces a variety of TTVs with different properties, namely, amplitudes, periodicity, and time span. TTVs signals in our catalog have amplitudes that range from almost zero minutes to 180 minutes, with an overall mean of ∼18 minutes. Mean TTVs amplitudes grouped by multiplicity reach 21, 21, 10, and 23 minutes for two-, three-, four-, and five-planet systems, respectively. It translates to a wide range of signal-to-noise ratios (S/N), ranging from ∼1 up to ∼400 with a mean of ∼23 (using the definition of the ratio of the TTVs amplitude and the timing uncertainty; Nesvorný 2019). There is also a variety of observation time spans, which range from ∼230 to ∼5800 days, with a mean of 1200 days.

As described in Section 2.1.1, Nauyaca requires specific boundaries for each planetary parameter in Θ. These boundaries delimit the parameter space that will be explored. Here, we discuss how to delimit the parameter space (including the establishment of constant parameters) before carrying out the recovery test over the mock catalog.

From fitting real transit observations it is possible to determine the planetary radius ratio, the orbital period from consecutive transits, and the impact parameter. In the case of the last two, we find from current data of observed exoplanets (from the Exoplanet Archive) that orbital periods are determined, on average, with uncertainties of ∼10⁻⁴ days. However, orbital inclinations are on average determined at a level better than 1° (also according to the Exoplanet Archive). Considering nearly coplanar orbits, the line of nodes of each planet are nearly aligned and mutual inclinations are close to 0°, which can help to reduce the dimensionality of the problem. Observed transiting exoplanets have been shown to have small mutual inclinations, typically below 3° (Fang & Margot 2012; Fabrycky et al. 2014), which has been demonstrated to have a negligible effects on TTVs (Nesvorný 2009; Nesvorný & Vokrouhlický 2014; Hadden & Lithwick 2016). Furthermore, because the orbital inclination is a well-restricted observational parameter (with a mean error of 078 with a dispersion of 22, according to data from the Exoplanet Archive), it can be kept as a fixed parameter. Additionally, in test runs we let the inclination vary between ∼80°–100° (which is a comprehensible range of inclinations that allows the transits) and confirmed that this angle has limited effects on the results for the transit time models. Nearly aligned ascending nodes represent prograde orbits and anti-aligned ones represent retrograde orbits. Thus, in order to model nearly coplanar orbits we kept fixed the orbital inclinations to their true values, inc_j,True. We also kept fixed the ascending node of the internal planets Ω_1,True (≈90°), which by construction of the coordinate system can be fixed.

For our recovery test, we delimited the planetary parameter space to these ranges, bracketing the lower and upper limits: mass [0.0123, 10 mass_j,True] M_⊕, P [P_True − δ t, P_True + δ t] days, ecc [10⁻⁵, 0.3], inc (fixed) [inc_True] deg, ω [0, 360] deg, M [0, 360] deg, and Ω [70, 110] deg. Here, the lower and upper boundaries for masses correspond to approximately a Moon mass and 10 times the true planetary mass, respectively. The boundaries in period P are around the true period of the planet with a width δ t corresponding to a typical observed period error given the orbital period itself. We estimated δ t by doing a linear regression between observed period errors as a function of the orbital period (data from the Exoplanet Archive), such that δ t = mP + b [days], where we determined m = 2.11 × 10⁻⁵ and b = 4.4019 × 10⁻⁵. The lower limit in eccentricity (ecc) was set to be small but different from 0 in order to avoid an undefined argument of periastron, while the upper limit was allowed to have values up to 0.3. This chosen range of eccentricities is compatible with 80% of the currently observed planet eccentricities. Argument of periastron (ω) and mean anomaly (M) take the usual definition between 0° and 360°. The boundaries of the ascending nodes (Ω) for all of the planets, except the internals, were restricted to a search radius around ≈ Ω_1,True ± 20 deg, and thus we consider only prograde orbits for simplicity.

We inspected the dependence of the parameters that govern the MCMC. The chosen parallel-tempering MCMC method (Vousden et al. 2015) is fine-tuned by two main parameters, namely, the number of temperatures (N_temps) and the maximum temperature (${T}_{\max }$) to build the temperature ladder. The sampling performance also depends on the number of walkers per temperature (N_w) and finally on the number of iterations (N_iter) per chain. We carried out many tests with combinations of these parameters and found that N_temps ≲ 15 with N_w ≲ 150 in a run over N_iter ∼ 5 × 10⁵ are enough in most cases to recover the true parameters Θ_j,True from the mock catalog within 1σ. We also note that ${T}_{\max }\sim {10}^{2}\mbox{--}{10}^{3}$ is a good choice for the maximum temperature. ${T}_{\max }=\infty$ (as suggested in Vousden et al. 2015) is not adequate for our purposes since walkers belonging to this temperature would propose steps outside our predefined boundaries. We confirmed this behavior on several occasions during our tests. Other parameters to control the dynamics of the temperature adjustment are internally defined within Nauyaca following the suggestions by Vousden et al. (2015). ⁵

For the recovery test presented in this work, we imposed uniform log-priors for simplicity with the functional form

$$\begin{eqnarray}\mathrm{log}({{ \mathcal P }}_{k})=\left\{\begin{array}{ll}0 & {b}_{\mathrm{low}}\leqslant k\leqslant {b}_{\mathrm{upp}}\\ -\infty & \mathrm{otherwise}\ \end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where b_low and b_upp correspond to the lower and upper boundaries for the kth planetary parameter. Table 1 shows a list of the model parameters for each planet as well as the adopted range of uniform priors. The choice of these validity ranges has been described previously in Section 4.1. From these parameters, inc was set as a constant to the true inclination, inc_True, taken from the mock catalog. For each system, Ω considers uniform priors except for the innermost planet (Planet 1 in the catalog entries), which is set to the true value Ω_1,True. Thereby, the posterior probability (up to a constant) is given by the sum of the log-likelihood (Equation (1)) and the log-prior (4) functions. In practice, any other prior functions can be supplied to Nauyaca to calculate the posterior probability.

Table 1. Model Parameters and the Selected Priors for the Transit Timing Fit

Parameter	Prior	Units
mass	${ \mathcal U }$(0.0123, 10 × massTrue )	M⊕
P	${ \mathcal U }$(PTrue − δ t, PTrue + δ t )	days
ecc	${ \mathcal U }$(1e-06, 0.3)
inc	incTrue	deg
ω	${ \mathcal U }$(0, 360)	deg
M	${ \mathcal U }$(0, 360)	deg
Ω	${ \mathcal U }$(70, 110) or Ω1,True	deg

Note. ${ \mathcal U }$[b_low, b_upp] denotes the uniform ranges between lower and upper boundaries. Single values are the fixed parameters in the simulations. Parameters with label True take the data values from the mock catalog (see the Appendix). See text for details.

Download table as:  ASCIITypeset image

Solving the inversion problem of a number of interacting planets, N_pla, implies the exploration of a parameter space of dimensions ∼7N_pla. The nature of the problem is computationally demanding since an N-body integration should be done at each iteration per walker. The wall-clock time also increases with the observational time span. Thereby, choosing a strategic initial guess for walkers adapted for the parallel-tempering MCMC is crucial to minimize these side effects.

Starting at random points from the prior function would be a reasonable choice assuming that we have informed previous knowledge about the parameters. In general, it would not be the case for planets characterized for the first time. This motivated us to use optimizer results to make an educated initial guess about the planetary parameters. Optimizers take advantage of having both a low computing time consumption in comparison with a full MCMC run (between 1% and 3% of the total time) and an ability to find many modes in the parameter space since realizations are independent among themselves. Although these solutions could be just rough approximations of more detailed solutions, they could help to identify high probability regions suitable for initialization. Initializing walkers near an optimum sensible place is better than using any random point in the parameter space (Hogg & Foreman-Mackey 2018). Even more, initialization using multistart local optimization results is found to enhance the exploration quality and be more suitable for multichain methods (e.g., Hug et al. 2013; Ballnus et al. 2017).

We present three strategies implemented in Nauyaca to initialize walkers using the information from the optimizers, namely, Gaussian, picked, and ladder. In order to set initial values from these strategies, we first made a filter over the total number of optimizer solutions (N_opt) by sorting them according to their χ². Then, we took a fraction (f_best; defined between 0 and 1) of that ordered list that includes the uppermost solution. That subset of solutions Θ_opt was then used to initialize walkers.

In Figure 4, we show examples of initialization using these strategies for many values of f_best. For the current example, we focus on the mass space of two planets identified with ID = pl2_id52 in Table 4. We performed 320 realizations of the optimizers to draw an initial walker population for a ladder of N_temps = 10 temperatures, considering different values of the parameter f_best. The individual solutions from Θ_opt are shown with squares in the first column, and the cyan star shows the position of the true mass. Colored dots are the initial walkers belonging to different temperatures. As f_best is reduced, solutions with high χ² are discarded. We detail these initialization strategies adapted to the parallel-tempering MCMC:

• 1.  
  Gaussian. Walkers are drawn from a Gaussian centered at the mean of Θ_opt with a 1σ value corresponding to the data dispersion, for each dimension. This is the simplest initialization method and possibly the most frequently used. The main difference is that here the mean and standard deviation are based on the independent random realizations and not on previous knowledge of the planetary parameters (for example, masses from radial velocity measurements). Thus, if the optimizers find a well-restricted solution for any dimension, the MCMC will be able to find the global high probability region faster in contrast to a uniform random initialization. From Figure 4, it is seen that using this strategy while reducing f_best could help to identify the zone in the parameter space where the global minimum could exist. Thus, by using this strategy most of the local modes are covered but there is not special attention given to the individual modes or to the possible correlations between parameters.
• 2.  
  Picked. Solutions from Θ_opt are randomly picked and the initial walkers are drawn from the vicinity of those solutions. If the optimizers find many modes, this strategy ensures that walkers will be drawn from around all of the modes which could correspond to local minima or the global minimum. From Figure 4, it is seen that this strategy confines the initial population of walkers as f_best is reduced while keeping the dependency between parameters. Here, walkers are equally distributed around any mode independent of their temperature and therefore all of the modes are initially sampled with the same frequency. This could be suitable for solving problems where apparently there is not a unique region where optimizer results agglomerate.
• 3.  
  Ladder. Solutions from Θ_opt are divided into an integer number of chunks equal to the number of temperatures. Walkers belonging to temperature 1 (the main temperature; blue dots) are drawn from the first chunk, which includes the uppermost solution (i.e., with the lower χ²). Walkers belonging to the second temperature are drawn from the first and second chunks. The same rule is followed for the rest of the temperatures until, finally, walkers for the hottest temperature (yellow dots) are drawn from around all of the solutions in Θ_opt. From Figure 4 it is seen that using this strategy, the modes with the lower χ² are highlighted as f_best is reduced. Unlike the picked strategy, ladder assigns the outstanding modes to colder walkers while hotter walkers sample the more disperse ones. It allows the exploration of other modes but avoids getting stuck in these local minima.

Zoom In Zoom Out Reset image size

The choice of the best strategy and parameter f_best can vary according to the problem. Ideally, a visual inspection of the optimizer solutions (as in Figure 4) could help to identify modes in the parameter space that allows us to decide how to initialize walkers. In practice, a statistical indicator (e.g., standard deviation) or a clustering method could be helpful to choose the initialization: Gaussian for high dispersed data, picked for highly multimodal parameter spaces, and ladder for a multimodal space with a main mode (as in the example of Figure 4).

Note however that the usage of the optimization module and the proposed initialization strategies is not a mandatory step in Nauyaca prior to the implementation of the MCMC. Any initial walker population can be provided by the user as long as these proposals are inside the physical boundaries described in Section 2.1.1. Nonetheless, we empirically find that following this heuristic procedure notably enhances the MCMC performance at a low computational cost.

We ran the optimizers to explore the parameter space with the aim of minimizing the differences between the observed and modeled transit times (Equation (2)). Since these runs are independent of each other, increasing the number of realizations enhances the chance of finding more minima that could correspond to local minima or the global minimum. We took a fraction (namely f_best) of the whole set of solutions with the best χ² to build a subset of solutions Θ_opt, which were used to initialize walkers.

Setting an optimum composition of Θ_opt is an interplay among N_opt, f_best, and the number of temperatures. Thus, there is not a unique way of defining these parameters. Although, we noticed that in most cases f_best < 10% is an appropriate fraction using N_opt ∼ 300–500. We used f_best shown in Table 2, according to planet multiplicity. That selection translates to 21, 25, 26, and 48 solutions that comprise Θ_opt for systems with two, three, four, and five planets, respectively. Taking higher values of f_best means taking more solutions from optimizers with increasing χ². Usually solutions with high χ² are located far from the meaningful modes, resembling random proposals (as can be noticed from Figure 4). Thus, selecting almost all of the optimizer solutions (f_best ∼ 1) could reduce the contribution of the optimizers to locate optimal initial regions to draw walkers.

In Figure 5, we show the results of the optimizers by measuring the differences between the solutions found and the true values. The solutions considered in this figure are part of the subset denoted previously as Θ_opt. Therefore, most of the solutions from the optimizers with higher χ² values were discarded. The color map indicates the percentage of solutions falling inside each bin. The darkest bins near to zero depict planets for which ≳50% of the solutions better approximate the true solutions. Subpanels with almost uniform yellow/white colors are those parameters less constrained by the optimizers.

Zoom In Zoom Out Reset image size

For the case of the planet masses, 90% of the solutions in Θ_opt are typically determined within a factor ∼2–6 with respect to the true masses. Even though the dispersion tends to increase with increasing mass and number of planets, the solutions approach the correct mass with a relative low dispersion. In the case of the periods, the solutions tend to span over all of the allowed space. We found this behavior is due to the narrowness between our selected lower and upper boundaries, which scales according to the period (see Section 4.1). We noticed that when the boundaries are enlarged (for example, with a width of ±0.01 days) the solutions agglomerate around the true periods. Eccentricities are in general well restricted only for two-planet systems, where the dispersion reaches up to ∼0.06 around the true values. For higher multiplicities the solutions scatter out, with a typical dispersion of ∼0.08. We found that optimizers tend to overestimate eccentricities as the number of planets increases. For the argument of periastron (ω) and mean anomaly (M), a similar behavior is found. They are better restricted only for two-planet systems, since for higher multiplicities the dispersion increases, tending to be evenly distributed over the parameter space. Ascending nodes are the less restricted angles given our established boundaries for this parameter. Note however that for the recovery test, we just considered limited prograde solutions (as described in Section 4.1).

Although the initial search region explored by the MCMC would be narrowed down as a result of the optimization process, the chains were initialized in these high density regions and explored all of the allowed parameter space since the adopted uniform priors were not modified.

In the work carried out by Nesvorný & Beaugé (2010) the authors test their tool TTVIM (which combines a minimization algorithm with an analytic approximation for inverting TTVs) using synthetic transit observations with the aim of recovering planet masses and orbital elements of a nontransiting planet. Nesvorný & Beaugé (2010) used upper limits in relative and absolute errors in planet parameters to decide whether the solutions were correctly recovered or not. Unlike Nesvorný & Beaugé (2010), we considered as recovered those solutions from the posterior distributions consistent with the true values within the 68% (1σ) and 95% (2σ) credible intervals of the total posteriors, assuming that they follow a normal distribution. These solutions are marked in Figure 6 as circles and squares. The comparison is made for planet parameters: mass, period (P), eccentricity (ecc), periastron longitude (ϖ; defined as the sum of ascending node and argument of periastron, ϖ = Ω + ω), and the true longitude (ℓ; defined as ℓ = ν + ϖ, where ν is the true anomaly, which is obtained from a Fourier expansion with the terms of the mean anomaly and eccentricity). Some of our posteriors exhibit more than one peak, manifesting the nature of the selected sampler. Median values and errors are usually centered on the prominent peaks so the data in Figure 6 are representative of our results. In practice, a dynamical stability test could help to distinguish the physically possible solutions. This study is underway, but it is out of the scope of the present work.

Zoom In Zoom Out Reset image size

From Figure 6, we take the discrete values of the differences between recovered and true data (excluding the unrecovered data marked as gray crosses), assuming that discrete medians are representative of the MCMC results. This was done for each planet multiplicity and planet parameter. From the resulting distributions, we measured the standard deviation that represents, on average, the typical precision achieved in the recovery test for different parameters and planet multiplicities. For masses, we found the minimum dispersion of 0.7 M_⊕ for the two-planet systems, and a maximum of 14 M_⊕ for five-planet systems. For P, a minimum dispersion of 9 s for two-planet systems, and a maximum of 110 s for five-planet systems. For ecc, a minimum dispersion of 0.007 for two-planet systems, and a maximum of 0.03 for three-planet systems. For ϖ, a minimum of 35° for five-planet systems, and a maximum of 50° for three-planet systems. For ℓ, a minimum dispersion of 3° for two-planet systems, and a maximum of 75 for four-planet systems. Intermediate dispersions were found for multiplicities not mentioned in these ranges.

We also calculated the Kolmogorov–Smirnov (KS) statistic over the whole distribution of true values and the posterior medians for all of the parameters for planets in the full catalog. We found, for masses, a KS statistic = 0.07 (p-value = 0.96); for P, KS statistic = 0.01 (p-value = 1.0); for ecc, KS statistic = 0.25 (p-value = 3×10⁻³); for ϖ KS statistic =0.17 (p-value = 0.11), and for ℓ, KS statistic = 0.03 (p-value = 0.99). Thus, we cannot reject the hypothesis that the distribution of masses, periods (P), and true longitudes (ℓ) are statistically drawn from the same input distributions.

Figure 7 shows the percentage of planets consistent with the true parameters within 1σ and 2σ, when grouping planets with the same multiplicity. In Table 3 we summarize these results also considering the global percentages for the full catalog (independent of multiplicity). For masses, periods, and true longitudes, the global recovery percentages are consistent with the expected statistical values, i.e., around ∼68% of the planet parameters are recovered within 1σ and ∼95% within 2σ. These parameters also exhibit consistency with the input catalog according to the KS test. However, we note that ecc and ϖ have similar recovery percentages but are far below these statistically expected values. For these two parameters we found, respectively, 52% and 55% within 1σ and 80% and 76% within 2σ (see Table 3).

Zoom In Zoom Out Reset image size

We investigate the possible causes of both the similarity and low recovery rates by inspecting the dependence between both parameters. We found that ∼45% of planets with ecc ≲ 0.05 recover both ecc and ϖ at the same time. The remaining planets do not recover at least one of these parameters. By contrast, planets with ecc ≳ 0.05 recover both parameters ∼90% of the times, showing that the low recovery rates occur mainly for low eccentricity orbits. In our catalog, most of the planets belonging to systems with three or more planets have eccentricities below ∼0.05, as shown in Figure 3. Even more, from Figure 7, it is seen that eccentricity is the only parameter that diminishes its recovery rate as the number of planets increases (teal lines with stars).

We find that the cause of the similarity and the global low recovery rate for ecc and ϖ is a result of the difficulty of finding a well-defined argument of periastron as the orbits tend to be more circular, which in our case occurs mainly for planets belonging to high multiplicity systems. We will show in Section 5.4 that the trend of the eccentricity shown in Figure 7 and the low recovery rate of ecc (and consequently ϖ) can be partially explained by the diminishing of the number of transits available for external planets belonging to high multiplicity systems.

In Figure 8 we illustrate the TTVs of two systems with simulated data of different qualities and their orbits. The top figure corresponds to the system with ID = pl2_id52 in Table 4 for which we estimate an S/N of 8.4 and 5.0 for planets 1 and 2, respectively. The bottom figure corresponds to the system with ID = pl3_id1 in Table 5 for which we estimate an S/N of 1.6, 2.08, and 2.06 for planets 1, 2, and 3, respectively. Left panels show 100 orbital configurations (thin colored orbits) reconstructed from the planetary parameters taken randomly from our MCMC posteriors. For comparison, the true orbits are plotted as solid black lines. Transits are detected each time the planets cross the +Z-axis, which corresponds to the line of sight. Right panels show the 100 TTV signals (colored solid lines) using the same planet parameters of the orbits. Synthetic observations are shown by empty circles with error bars. For these two systems, the TTV models are consistent with the synthetic data, but with a better planetary determination for the top system, which has a better S/N. Hence, the quality of the transit times directly affects the determination of the orbital configurations. For example, we identify that in the case of the three-planet system, the argument of periastron shows a wide range of possible solutions with respect to the true position (dashed black lines).

Zoom In Zoom Out Reset image size

We statistically measured the fitted residuals from the whole catalog by taking 100 random samples per planet from the MCMC posteriors and then calculating their corresponding TTVs. We took the differences between data and the fitted transit times, grouping these residuals according to their multiplicity. The resulting distributions are shown in Figure 9. We found a typical standard deviation for these residuals of 2.5 minutes for two-planet, 2.3 minutes for three-planet, 3 minutes for four-planet, and 5.3 minutes for five-planet systems.

Zoom In Zoom Out Reset image size

The TTVs inversion problem challenges not only the methods or algorithms, but also the observational data requirements. In this subsection we address this issue by considering the quality (measured by the S/N) and quantity (N_tran) of the data required in order to correctly determine planetary parameters. We focus on the determination of planetary mass and eccentricity.

Previous works have proposed different answers regarding the required number of transits and data quality to invert TTVs. Saad-Olivera et al. (2019) suggested that a large number of TTVs combined with a high S/N can be enough to robustly estimate planetary parameters without the necessity of radial velocity measurements. A detailed analysis regarding the characterization of nontransiting planets was conducted by Veras et al. (2011), finding that at least 50 transits are needed to invert the problem. Nesvorný & Morbidelli (2008) developed a method based on perturbation theory to characterize two-planet systems, and found that an S/N ∼ 15–30 and about N_tran ≳20 are typically required to uniquely characterize these systems.

We carried out an analysis of the parameters involved in the fitting procedure, looking for possible correlations between the data quality and quantity, the intrinsic planetary parameters (including the number of planets, N_pla), and the results we found using Nauyaca. A first inspection suggested a trend with the number of transits and thus it is convenient to define the total signal-to-noise ratio (S/N)_T as the product of ${\rm{S}}/{\rm{N}}\times \sqrt{{N}_{\mathrm{tran}}}$. Figure 10 shows the fractional error in the determination of the planet mass and eccentricity, and the dependency with the (S/N)_T , N_tran, and N_pla. Here, 2σ_mass and 2σ_ecc correspond to our derived uncertainties taken from the posteriors encompassing the 95% credible interval. These uncertainties are divided by the true mass and eccentricity, and thus, Figure 10 can work as a guide to put some constraints on the determination of these parameters given the quality and quantity of the data.

Zoom In Zoom Out Reset image size

In general, it is seen that the fractional error (related with the adequate determination of the parameters) has a correlation with the quality of the data and in a second term with the quantity. Moreover, the deficiency of one of these quantities is in some cases compensated for by the other one. We observed that planets with individual high (S/N)_T can constrain the planetary mass and eccentricity even with a relative low N_tran. Complementary, a large N_tran can compensate for a low (S/N)_T . Planets with a low (S/N)_T exhibit a larger dispersion of their fractional errors for their masses and eccentricities, although, those with a large N_tran usually have a better determination (i.e., a lower fractional error). Comparing the top and bottom panels, it can be discerned that the majority of the unrecovered planets are those with a number of transits ≲60–80 (white/yellow crosses), corresponding to systems with three, four, and five planets. However, planets of these multiplicities with low (S/N)_T can still find the correct solution but with a poor constraint (white/yellow circles in the upper left part of the top diagrams). Also notice from the bottom panels that the trend between the fractional errors and the total (S/N)_T seems to hold when grouped by planet multiplicities.

From these results it is seen that (S/N)_T and N_tran have a significant impact on the determination of these parameters. The low recovery rate for ecc discussed in Section 5.2 is then explained by the low (S/N)_T of our synthetic data and the low number of transits per planet (≲20–40) for systems with many planets. Low (S/N)_T can be due to noisy measurements of the transit timing or due to low amplitudes of the TTVs signals. In the first case, the limitation is in part caused by the instruments and observational strategy. In the second case, low amplitudes correspond to planets in less perturbed (almost circular) orbits. In this situation, inverting TTV signals results in less constrained orbits.