FAST INVERSION METHOD FOR DETERMINATION OF PLANETARY PARAMETERS FROM TRANSIT TIMING VARIATIONS

In Nesvorný & Morbidelli (2008) and Nesvorný (2009) (hereafter NM08 and N09), we developed and tested a fast inversion method that can be used to characterize planetary systems from the observed transit timing variations (TTVs; Agol et al. 2005; Holman & Murray 2005). See NM08 and N09 for a technical description of the method. Here we use this new method to solve the TTV inversion problem for an arbitrary planetary system. The results provide a baseline for studies of real exoplanetary systems for which TTVs will be detected. Examples of past work that would greatly benefit from the application of the fast inversion algorithm discussed here include Steffen & Agol (2005), Agol & Steffen (2007), Miller-Ricci et al. (2008), and Gibson et al. (2009).

In Section 2, we briefly describe the TTV inversion method. In Section 3, we apply it to a case with coplanar planetary orbits. Inclined planetary orbits are discussed in Section 4. We show, for example, that the mutual inclination of planetary orbits can be determined from TTVs. This important parameter, which may be used to test planet-migration theories (e.g., Rasio & Ford 1996; Goldreich & Sari 2003), is not typically available from other existing planet-detection methods.

Our TTV inversion method, hereafter TTVIM, has two parts. The first part is a fast algorithm for the computation of transit times, (δt_j)_trial, 1 ⩽ j ⩽ N, for specified planetary parameters. This algorithm is based on perturbation theory (NM08, N09). It calculates the short-period TTVs as these have been shown to be the most diagnostic (NM08). The long-term effects such as the apsidal precession produced by the perturbing planet are more difficult to detect if transit observations span only a few years (Miralda-Escudé 2002; Heyl & Gladman 2007).

The second part of TTVIM is an adaptation of the downhill simplex method (DSM; Press et al. 1992). The DSM is used to search for the minima of

$$\begin{equation} \chi ^2 = \sum _{j=1}^N \left[ (\delta t_j)_{\rm trial}-(\delta t_j)_{\rm obs} \over \sigma _j \right]^2\,, \end{equation} \tag{ 1 }$$

where (δt_j)_trial are the transit times produced by the first part of the algorithm for a trial planetary system, (δt_j)_obs are the observed mid-transit times, and σ_j are the measurement errors. It is assumed here (as indicated by δ's) that the period of the transiting planet, P, has been removed from transit observations. Thus, (δt_j)_obs = (t_j)_obs − C_jP, where (t_j)_obs are the actual transit times and integer C_j denotes the transit cycle.

The best-fit planetary parameters correspond to the global minimum of χ², denoted by χ²_min in the following. A large number of initial trials must be used to assure that the DSM method finds χ²_min. The confidence levels for the normally distributed data can be defined as Δχ² = χ² − χ²_min < (Δχ²)_cut, where the (Δχ²)_cut values are properly chosen for N and the required confidence level (NM08).

Here we assume that the mass and orbit of the transiting planet are known from transit and radial velocity (RV) measurements as this should be the most common case in practice. If so, χ² is a function of seven unknown parameters of the perturbing planet, χ² = χ²(m₂, a₂, e₂, i₂, Ω₂, ϖ₂, λ₂), where m₂ is the mass, a₂ is the semimajor axis, e₂ is the eccentricity, i₂ is the inclination, Ω₂ is the nodal longitude, ϖ₂ is the periapse longitude, and λ₂ is the mean orbital phase at t = 0 (arbitrarily defined here to correspond to cycle C₀).³ The parameters of the transiting planet will be denoted by index 1. DSM must therefore search in seven-dimensional space for the global minimum of χ². This is not a trivial task because χ² often has many deep and narrow local minima (Steffen & Agol 2007). Fortunately, several simplifications can be made.

First, as the amplitude of the short-period TTVs scales linearly with m₂, we can calculate the TTV profile for the selected a₂, e₂, i₂, Ω₂, ϖ₂, λ₂ values and obtain m₂ by the linear least-square fit. Second, the determination of χ² for a new set of the Ω₂, ϖ₂, and λ₂ values is computationally cheap in the perturbation algorithm, if χ² was determined previously for the required a₂, e₂, and i₂ values.⁴ The code can thus efficiently search for the minimum of χ²(a₂, e₂, i₂) for any value of Ω₂, ϖ₂, and λ₂. In practice, we use 5–20 values between 0° and 360° to resolve each of these parameters. This effectively reduces the number of dimensions to three. Once the interval of estimated a₂, e₂, and i₂ values is narrowed down, the solution can be refined by using the full DSM search in seven dimensions.

The tricky part of TTVIM is the choice of the initial guess in the (a₂, e₂, i₂) space. By trials and errors we found that fine (and nonlinear) sampling of a₂ is generally needed for a successful convergence of the algorithm. The best results were obtained with uniform sampling in 1/α², where α = a₁/a₂ < 1. Parameters e₂ and i₂ require less care since DSM usually finds the right minimum even in the high-e₂ and/or high-i₂ case if at least one corner of the initial simplex is stretched to e₂ > 0.2 and i₂ > 30°.

With the nominal setup, our TTVIM code (ttvim.f) requires about 2 minutes of CPU time⁵ for a coplanar fit with i₂ = 0 and about 50 minutes for the full seven-dimensional fit. In the absence of measurement errors, the success rate in finding χ²_min is better than 95%. Thus, ttvim.f is a robust code that can reliably solve the TTV inverse problem at a low computational cost.

We used a random number generator to define different sets of parameters m₂, a₂, e₂, i₂, Ω₂, ϖ₂, and λ₂. Typically, 1000–2000 different planet parameter sets were used in tests. In each case, the orbital evolution of the two planets was followed for a fixed timespan, 0 < t < T_int, with the Bulirsch–Stoer integrator (Press et al. 1992). During this timespan we interpolated for and recorded all transit times of the inner planet. These data mimic the real observations, (δt_j)_obs. They were used in a blind test where we applied the TTVIM code to each of these cases in an attempt to recover the original mass and the orbital elements of the non-transiting planet.

We start by discussing the case with star's mass m₀ = M_☉, where M_☉ = 2 × 10³³ g is the mass of the Sun, m₁ = 10⁻³ m₀, a₁ = 0.1 AU, e₁ = i₁ = 0, and N = 100 consecutive transits. Since NM08 and N09 showed that the behavior of the inversion method is insensitive to m₁, we will not test different m₁ values in this work. To distinguish between the issues related to the intrinsic limitations of ttvim.f and those arising from the finite precision of the real measurements, we first discuss an idealized case with zero measurement noise.

For coplanar orbits and σ_j = 0, the TTVIM code finds the correct planetary parameters with a high rate of success (Figure 1). The typical precision in the successful cases is |m₂ − m*₂|/m₂ < 0.2, |a₂ − a*₂|/a₂ < 0.02, |e₂ − e*₂| < 0.02, |ϖ₂ − ϖ*₂| < 10° and |λ₂ − λ*₂| < 10°, where the asterisk denotes the values determined by the TTVIM code.⁶ This is very satisfactory. In the absence of measurement errors, the result of the TTVIM code illustrated in Figure 1 with m₂ = 10⁻⁴m₀ is insensitive to the actual value of m₂.

Zoom In Zoom Out Reset image size

The main failure mode of the TTVIM code occurs near mean motion resonances between planets, because resonant perturbations are not (yet) taken into account in TTVIM. While the resonant signal can improve our chances of the TTV detection for (near-)resonant planets, it seems less useful in helping us estimate the mass and orbit of the planetary companion. Specifically, the amplitude and period of the resonant signal can be fit by a number of different planetary setups corresponding to different resonances. Thus, without an a priori knowledge of the mean motion resonance that is responsible for the observed behavior, the inversion problem from resonant frequencies alone is strongly degenerate.

Fortunately, the short-period TTVs underlying the resonant signal can still be used to determine the planetary parameters without much ambiguity. As shown in NM08, probably the best strategy is to isolate short-period frequencies in the signal by Fourier filter and apply the inversion method to the filtered signal. The application of this procedure is straightforward in individual cases (see NM08), where the resonant period, and thus the appropriate frequency cutoff, can be estimated from (t_j)_obs. We verified that this procedure works quite well in >75% of cases shown in Figure 1 in which the resonant variations are an issue.

The remaining <25% unsuccessful cases (representing <5% overall) correspond to the very large values of e₂ for which the Laplacian expansion of the perturbing function in TTVIM is not convergent (NM08), and/or planetary configurations that are not Hill stable. Direct N-body integrations can be used to address the TTV inversion problem in the very high eccentricity domain, but the CPU cost of these tests is likely to be substantial and lies beyond the scope of this Letter.

The measurement errors have a profound effect on the uniqueness of the inverse problem (Figure 2). For m₂ = 10⁻⁴ m₀, N = 100, and σ_j = σ = 3 s, corresponding to the Kepler-like precision of timing measurements for a Sun-like star with a 2 Neptune-mass planet, unique determination of planetary parameters can be achieved for most stable systems with q₂ = a₂(1 − e₂) < 3.3a₁, while for σ_j = σ = 10 s (Corot-like precision), it is required that q₂ < 2.6a₁. These limits approximately correspond to the planetary parameters for which the amplitude of the short-period TTVs is comparable to σ.

Zoom In Zoom Out Reset image size

Figure 3 shows the result of the TTVIM code for an Earth-mass planet. The region of parameter space in which unique determination can be achieved from TTVs is relatively small even with σ = 1 s. Thus, an Earth-mass planet detection and characterization of its orbit will require a rather fortuitous setup of the planetary system, in which (q₂ − a₁)/a₁ ≲ 2 (for external perturber).

Zoom In Zoom Out Reset image size

We applied the TTVIM code to mock planetary systems with 0° < i₂ < 50°. As in Section 3, we assumed that m₀ = M_☉, m₁ = 10⁻³ m₀, a₁ = 0.1 AU, e₁ = 0 and used N = 100 consecutive transits. Figure 4 shows the result for σ_j = σ = 0. The (a, e) plot does not differ much from the coplanar case although it may be noted that the quality of fits slightly degraded for the large e₂ values. This is probably related to the convergence problems of the perturbation algorithm in TTVIM. A precise N-body integrator should perform better for high e₂, although it has yet to be shown that an N-body integrator can be applied to the inclined inverse problem in practice due to the large CPU cost.

Zoom In Zoom Out Reset image size

Probably the most exciting result obtained in this work is that it was possible to determine the mutual inclination of planets for most planetary systems (Figure 4(b)). Unlike the RV technique, the TTV method can therefore be used to characterize the inclination distribution of multi-planet systems. Figure 5 shows the detailed statistic of TTVIM errors in i₂. In most cases, |i₂ − i*₂| < 2°. The tail of larger |i₂ − i*₂| values corresponds to the high-eccentricity cases. If the statistic is limited to q₂ = a₂(1 − e₂)>0.25 AU (Figure 5(a); dashed line), the fraction of successful cases with |i₂ − i*₂| < 2° increases to >90%. In the successful cases, orbital angles Ω₂, ϖ₂, and λ₂ are generally correctly determined to within a better than 5° precision.

Zoom In Zoom Out Reset image size

We also studied how the uniqueness of the inclined inverse problem is affected observational errors. The trends seen in these tests are very similar to those described in Section 3. Namely, the instrumental noise sets an upper limit on q₂ beyond which the determination of planetary parameters from TTVs is ambiguous. Again, we see that these limits approximately correspond to the planetary parameters for which the amplitude of the short-period TTVs is comparable to σ. The results in N08 can therefore be used to estimate whether (or not) a unique characterization of the specific inclined planetary system may be achieved from TTV observations with given σ.

We find that TTVs obtained in the coplanar case represent a good approximation of the TTVs for planetary orbits with i₂ < 20°. The planar version of the TTVIM algorithm can therefore be used in these cases to estimate the a₂ and e₂ values of the perturbing planet. This helps to narrow the range of initial guesses for the seven-dimensional fit and represents a factor of ∼20 speed up of the inversion. The full seven-dimensional algorithm needs to be used for i₂ > 20°.

The method developed here can be used to analyze TTVs found for any of the potentially hundreds of planets expected to be discovered by Kepler (Beatty & Gaudi 2008). Kepler should be able to detect TTVs of only a few seconds (Holman & Murray 2005), which should easily exist in many systems, extrapolating from the RV planets (Agol et al. 2005; Fabrycky 2009).

Perhaps the most interesting result that comes out of this work is that the shape of the TTV signal is generally sensitive to the orbital inclination of the non-transiting planetary companion. Thus, the TTV method can provide means of determining mutual inclinations in systems in which at least one planet is transiting. This parameter cannot be determined by other planet-detection methods.

TTVIM algorithm can be easily extended to incorporate uncertainties in the transiting planet's parameters. This can be done by sampling dimensions that correspond to the additional parameters. For example, in N09 we extended the NM08 method to the case with e₁ ≠ 0. This may be especially relevant to the transiting planets that will be discovered by Kepler because these planets are expected to have wider orbits, which are less susceptible to the circularizing effects of tides. The low CPU cost of the TTVIM algorithm is the key element which will make such studies possible.

This work was supported by the NSF AAG program.