First order TTVs

# Transit timing to first order in eccentricity

## Abstract

Characterization of transiting planets with transit timing variations (TTVs) requires understanding how to translate the observed TTVs into masses and orbital elements of the planets. This can be challenging in multi-planet transiting systems, but fortunately these systems tend to be nearly plane-parallel and low eccentricity. Here we present a novel derivation of analytic formulae for TTVs that are accurate to first order in the planet-star mass ratios and in the orbital eccentricities. These formulae are accurate in proximity to first order resonances, as well as away from resonance, and compare well with more computationally expensive N-body integrations in the low eccentricity, low mass-ratio regime when applied to simulated and to actual multi-transiting Kepler planet systems. We make code available for implementing these formulae.

planets and satellites: detection
1

## 1. Introduction

No planet orbits on a precisely Keplerian orbit: post-Newtonian corrections, stellar oblateness, and, most importantly, planetary perturbations cause deviations from a periodic ephemeris for transiting exoplanets (Miralda-Escudé, 2002; Schneider, 2003, 2004; Holman & Murray, 2005; Agol et al., 2005; Heyl & Gladman, 2007; Nesvorný & Morbidelli, 2008; Fabrycky, 2010). Transit-timing variations (TTVs) have been used to confirm that transit signals are in fact due to planets (Holman et al., 2010; Ford et al., 2012b; Steffen et al., 2012, 2013; Fabrycky et al., 2012; Ford et al., 2012a; Xie, 2013; Xie et al., 2014), to detect and characterize non-transiting planets (Ballard et al., 2011; Nesvorný et al., 2012, 2013), and to make precise measurements of the masses and dynamical states of multi-transiting exoplanet systems (e.g. Carter et al., 2012).

For the latter two applications, fast computation of TTVs is required for rapid searching through parameter space for perturbing companions, and for rapid computation of the posterior distributions of the masses and orbital elements of transiting planet systems. Numerical computation of TTVs can be sped up through symplectic integration, through a more efficient numerical solution of Kepler’s equation, and through transit time interpolation (Deck et al., 2014); however, this approach can still be too computationally intensive for high multiplicity systems, and does not pinpoint the physical origin of constraints upon planetary system properties. Analytic formulae based on perturbation theory can greatly speed computation, but the perturbation theory to high order in eccentricity and inclination becomes complicated quickly, and numerical codes that implement the analytic formulae have not been released or widely used (Nesvorný & Morbidelli, 2008; Nesvorný, 2009; Nesvorný & Beaugé, 2010). Much can be accomplished with first order (in eccentricity) perturbation theory because orbital eccentricities of many planets exhibiting TTVs are small. TTVs are most easily observed for pairs of planets near a mean motion resonance; thanks to the low eccentricity of the systems in consideration (Fabrycky et al., 2014; Hadden & Lithwick, 2014; Limbach & Turner, 2014; Van Eylen & Albrecht, 2015), the first order resonances are most represented among TTV pairs, and it is for first order resonances that a first order theory is adequate. TTVs caused by these first order resonant interactions are primarily sinusoidal and are subject to a degeneracy between mass and eccentricity (Boué et al., 2012), caused by mixing of two frequencies of perturbation which are aliased at the frequency of the transiting planet, as explained in an elegant analysis by Lithwick et al. (2012, hereafter LXW12). To break this degeneracy requires the measurement of additional modes, such as the short-timescale TTVs known as ‘chopping’ variations (Holman et al., 2010; Deck & Agol, 2015), or statistical analysis of many systems (Wu & Lithwick, 2013; Hadden & Lithwick, 2014; Xie, 2014). In addition, the LXW12 analysis is only approximate, and breaks down for pairs further from resonance (Deck & Agol, 2015). These issues motivate the current paper in which we derive an explicit formula for TTVs accurate to first order in eccentricity and planet-star mass ratio, valid for (nearly) plane-parallel transiting planets (although Nesvorný & Vokrouhlický 2014 showed that mutual inclinations of planets can be large and still be well described by coplanar TTVs).

We expect that these results will be useful a) for determining how different frequencies within the TTV signal constrain the planetary masses and orbital elements (Deck & Agol, 2015); b) for analyzing systems with a large number of interacting transiting planets by making linear additions of the analytic formula for pairs of planets (Lissauer et al., 2013; Jontof-Hutter et al., 2014); and c) for rapid search through parameter space of perturbing planets.

We first summarize the TTV solution to first order in eccentricity and mass in §2 (the full derivation is given in appendix A new approach to TTVs). We then compare these with prior results, both analytic and numeric (§3), including comparison of analytic and numeric analyses of specific systems. We discuss the numerical implementation and speed in §4. We end with a discussion of the possible applications and future directions (§5).

## 2. First-order solution:

Here we give a complete summary of the assumptions and variables used, and the solution to the first-order equations for readers that wish to simply use the results of this computation. The details of the derivation are given in Appendix A new approach to TTVs. Since multi-planet transiting systems typically have nearly edge-on orbits and hence small mutual inclinations, it is usually sufficient in analytic approximations to treat the problem in the plane-parallel approximation. This leaves four orbital elements for each planet (semi-major axis, , mean longitude, , eccentricity, , and longitude of periastron, ), plus the mass ratio of each planet to the star, , where is the mass of the inner planet, is the mass of the outer planet, and is the mass of the star. For nearly circular planetary orbits, there are two small dimensionless parameters in the problem: and . The usual procedure for computing transit timing variations is to: 1) to write down a Hamiltonian (or disturbing function) for perturbations due to another planet; 2) expand the Hamiltonian as a function of the orbital elements to the order in eccentricity desired plus one (e.g. if a transit timing solution is needed to first order in eccentricity, then the Hamiltonian must be expanded to second order in eccentricity), including the linear combinations of mean-longitudes leading to the important resonant terms necessary for sufficient accuracy; 3) compute the variation in the orbital elements using Hamilton’s equations, which are four first-order partial differential equations for each planet, and involves differentiating the Hamiltonian with respect to the orbital elements (which can be a rather complex operation); 4) integrate the resulting equations as a function of time; 5) compute the true longitudes, , as a function of time, where is the unperturbed Keplerian orbit, and is the perturbation of the th planet caused by its planet companion(s); 6) compute the transit timing variations:

 δti=−˙θ−1i,Kδθi. (1)

This is the approach taken by Agol et al. (2005), Nesvorný & Vokrouhlický (2014), and LXW12. A different approach employing Hamiltonian perturbation theory (Nesvorný & Morbidelli, 2008) was used in Deck & Agol (2015). This involves determining the canonical transformation between the full canonical orbital element set and the average set; the TTVs, which are deviations from an average “Keplerian” orbit, can be derived from this transformation.

The standard procedure outlined in detail above (based on Hamilton’s equations) has the advantages of requiring only first-order differential equations for the computation, and the advantage of using standard methods in celestial mechanics for the computation. However, there are two possible drawbacks: 1) the expansion of the Hamiltonian in orbital elements can be rather complex; 2) the main quantity of interest for transit-timing variations is , rather than the perturbed orbital elements. The derivation based on canonical transformations (Nesvorný & Morbidelli, 2008), though elegant, has the disadvantage of requiring the extra machinery and knowledge of Hamiltonian perturbation theory.

In our new derivation we forgo computing the orbital elements, and simply treat the problem in polar coordinates . We then use Newton’s equations in terms of a disturbing function which can be expressed as a function of polar coordinates, with the added advantage that Newton’s equations make clearer which forces are causing the perturbations. This approach has some possible advantages: 1) only two differential equations are necessary (albeit second-order rather than first-order); 2) the derivatives of the disturbing function with respect to the polar coordinates are easy to compute; 3) the perturbed polar coordinates directly yield the transit timing variations; 4) the resulting expression is more compact than in the Hamiltonian formulation. The second-order differential equation may seem like a drawback, but it can be solved using complex notation (as in LXW12) and by expanding the derivatives of the disturbing function in terms of orbital elements, which yields harmonic functions which are easy to integrate. The final answer is expressed as a sum over harmonics of the perturbing planet’s orbital frequency (Deck & Agol, 2015). Each coefficient for each planet in the harmonic series solution can be solved for by inverting three two-by-two matrices, which have a standard format, resulting directly in the transit timing variations at a particular frequency.

The unperturbed orbital frequencies, , are defined by . As usual, . We define where is the Laplace coefficient,

 b(j)1/2(α)=1π∫2π0dθcos(jθ)(1+α2−2αcosθ)1/2. (2)

The difference in mean longitude of the planets is . Auxiliary dimensionless quantities are:

 βj =j(n1−n2)/n1=j(1−α3/2), (3) κj =j(n1−n2)/n2=j(α−3/2−1). (4)

The functions we use below are given by:

 ~Aj00 = b(j)1/2(α), (5) ~Aj10 = α∂b(j)1/2/∂α, (6) ~Aj20 = α2∂2b(j)1/2/∂α2, (7) ~Aj01 = −(~Aj10+~Aj00)=−(α∂b(j)1/2/∂α+b(j)1/2), (8) ~Aj02 = 2~Aj00+4~Aj10+~Aj20=2b(j)1/2+4α∂b(j)1/2/∂α+α2∂2b(j)1/2/∂α2, (9) ~Aj11 = −(2~Aj10+~Aj20)=−2α∂b(j)1/2/∂α−α2∂2b(j)1/2/∂α2. (10)

To use this solution in computing TTVs, the longitudes need to be computed from the observed transit times; the mean ephemeris, , may be used in computing the (unperturbed) orbital ephemeris. Now, the mean longitudes are given to first order in eccentricity by

 λi=2π(t−t0,iPi)+2eisinϖi, (11)

if we assume that the orbital reference is along the line of sight, and thus at the times of transit.

The solutions for the inner planet () and outer planet () are given by:

 δt1 = P12πμ2∑j≥1[f(0)1,jsin(jψ)+f(−1)1,je1sin[jψ−(λ1−ϖ1)]+f(+1)1,je1sin[jψ+(λ1−ϖ1)] (12) + f(−2)1,j−1e2sin[jψ−(λ1−ϖ2)]+f(+2)1,j+1e2sin[jψ+(λ1−ϖ2)]], (13) δt2 = P22πμ1∑j≥1[f(0)2,jsin(jψ)+f(−2)2,je2sin[jψ−(λ2−ϖ2)]+f(+2)2,je2sin[jψ+(λ2−ϖ2)] (14) + f(−1)2,j+1e1sin[jψ−(λ2−ϖ1)]+f(+1)2,j−1e1sin[jψ+(λ2−ϖ1)]], (15)

where the functions are given by:

 f(±k)i,j(α) = u(γ,c1,c2)+δikv±(ζ,d1,d2), (16) u(γ,c1,c2) = (3+γ2)c1+2γc2γ2(1−γ2) (17) v±(ζ,d1,d2) = (±(1−ζ2)+6ζ)d1+(2+ζ2)d2ζ(1−ζ2)(ζ±1)(ζ±2). (18)

where , , and and , , and are given in Table 1, and is the Kronecker delta function. Note that the top signs in correspond to values, while the bottom correspond to . The functions are solely a function of , , and .

In practice the sum over from 1 to must be truncated at a finite value of . Typically does not need to be chosen to be too large since the Laplace coefficients decline in amplitude with (Deck & Agol, 2015). We recommend choosing a large enough such that the resulting computation is converged.

As an example of using Table 1, the coefficient has , , and , , , and . Then, the coefficient is given by:

 f(−1)1,j(α)=u(βj−1,c1,c2)+v−(βj,d1,d2). (19)

## 3. Comparison with other formulae

The zeroth-order solution (in the limit ) compares exactly with the results given in Agol et al. (2005), Nesvorný & Vokrouhlický (2014), and Deck & Agol (2015) which were derived with Hamilton’s equations and the approach based on canonical transformations; this is reassuring given the very different approach used in this derivation. We have also rederived the first-order eccentricity equations using the approach based on canonical perturbation theory employed in Deck & Agol (2015), and found exact agreement with the results presented here to first order in eccentricity.

In Figures (1) and (2) we plot the eccentricity-dependent coefficients, , as a function of period ratio, . The zeroth-order eccentricity coefficients are plotted in Deck & Agol (2015). The first-order coefficients can show three singularities for the terms with superscripts and near first-order resonance, second-order resonance, and .

### 3.1. Comparison with first order resonant equations

LXW12 present a formula valid near (but not in) first order mean-motion resonances that captures the behavior of resonant terms in an elegant, but approximate, manner. Here we compare the complete formulae given here to their near-resonant formulae.

The expressions for and do not show the same dependence in the denominator as the expressions in LXW12; their expression just contains the resonant frequency, , while ours contains additional frequencies. We carried out the partial fraction expansion of to isolate the denominator which matches LXW12’s expression, and we find that the expressions agree exactly with their expressions (A28) and (A29).

To compare our full expression with LXW12’s, we have computed the eccentricity-dependent (near 2:3) expression for the inner and outer planets (as this term is unaffected by the indirect terms). We re-write the TTV formulae derived here and in LXW12 for the th planet perturbed by planet as

 δti =μkni∑j≥1[Aji,ieisin(jλk+ϕji,i)+Aji,keksin(jλk+ϕji,k)] (20)

where we have set at transit (this incurs some error, at order , but that is a second order effect since we are comparing the TTV term linear in ). When written in this way, the amplitude and phase depend only on , and .

For the 3:2 resonance, the LXW12 resonant term depends on and , which both decline quickly away from resonance, and thus other terms that depend on and make a more significant contribution further from resonance; hence our formulae agree close to resonance but diverge away from exact commensurability. Figure 3 shows the fractional error in the amplitude and phase of the terms that are proportional to the eccentricities of the planets for radians. The error in the LXW12 expression (A28) and (A29) reaches 20% at a 5% separation from exact resonance in this case; if we use the further approximate expression given in their main paper in lieu of (A28), the discrepancy for the inner planet increases to 30%. This is similar to the error in the zeroth-order component of their expression (Deck & Agol, 2015). The zeroth- and first-order eccentricity terms have a different dependence on longitudes: for example, for the inner planet the zeroth order term scales as , while the first order term scales as , where . Since the mean longitude of the inner planet is nearly identical at each transit of the inner planet, the term in the exponent is approximately constant, while the dependence is identical in both terms, leading to aliasing of these coefficients. Thus the error incurred in their approximation can lead to a different phase dependence and amplitude for this aliased term away from resonance.

### 3.2. Comparison with N-body integrations

We have carried out extensive integrations of three-body systems using TTVFast (Deck et al., 2014), and compared the results with the first-order analytic formulae (equation 12). Note that the TTVFast code uses the convention of the longitude of periastron being measured from the sky plane to match the convention of radial velocity surveys, while here we use the observer’s line of sight as the reference direction, as done in LXW12 and Deck & Agol (2015). The longitudes computed from the TTVFast code need to have subtracted to make the plots shown below. In addition, the orbital elements accepted by TTVFast are the instantaneous/osculating orbital elements (initial conditions) at the specified initial time, while the orbital elements used in these formulae are the mean orbital elements of the planets over the timescale of the observations.

#### Eccentricity and period ratio dependence

Figure 4 compares the precision of the analytic formula as a function of and eccentricity of both planets, which are set to be equal, . We have set , which we found (approximately) maximizes the discrepancy of the analytic model compared with the N-body model, and which (approximately) minimizes the discrepancy; hence the figures bracket the precision of the analytic model. This is due to the fact that the anti-aligned longitude geometry causes the planets to be closer at conjunctions that occur when the inner planet is at apoapse and the outer is at periapse; their proximity at these conjunctions causes their gravitational interactions to be more sensitive to deviations from the epicyclic approximation, which are second order in eccentricity, and thus missing from our computation. We have assumed that the period of the inner planet is days, and we have integrated the system with TTVFast for 1600 days, about the duration of the initial Kepler mission, assuming plane-parallel orbits. For these tests we assume , and we selected random values for the longitudes of the planets at the initial time. For each set of initial conditions, we output the orbital elements at regular intervals during the N-body integration, from which we computed the average orbital elements over the duration of the integration. These averaged orbital elements were used for computing the amplitudes of the analytic model, which we summed up to . We optimized the fit of the analytic formula to the numerical TTVs by allowing the ephemerides of the planets to vary in the formula, but holding the eccentricity vectors and mass ratios fixed at the values computed from the time-averaged N-body simulation, while we computed in the analytic formula from the ratio of the periods derived from the best-fit ephemerides.

The fractional precision was computed from the RMS of the residuals of the best analytic model fit to the TTVs, divided by the RMS of the TTVs computed from the N-body integration. Figure 4 shows that the formula works to better than 10% precision for a wide range of parameter space. However, it fails near resonances, most significantly for the : resonances indicated in green, and : in blue. For the outer planet, there are narrow regions near : period ratios for which the formula does poorly. The disagreement grows in breadth for larger eccentricities. This diagram can be used to pinpoint the relevance of the analytic formula for a particular system, and we suggest that the analytic formulae should be used with caution in the regions where the formula disagrees by more than 10% precision.

Most of the regions where the formula fails are near resonance. In these cases, the residuals can frequently be fit by sinusoidal variations at the relevant resonant frequencies of the higher order resonant terms that are not captured in the first-order model; when including these sinusoidal terms in the fit, the residuals drop dramatically near the resonances. Thus, the analytic first-order solution plus a sinusoid with arbitrary amplitude and phase can be used for systems in which only the shape of the transit timing variations plus the specific variations of the non-resonant terms is necessary (although this approach breaks down for large enough eccentricity).

#### Mass dependence

We have carried out simulations for a range of masses, keeping . We find that the fractional error of the analytic formula grows near resonances and near as the mass increases with weaker dependence upon eccentricity. For small , the formula works well up to , while near : period ratio, for example, the error broadens around the resonance.

Figure 5 shows the fractional error in the formula (computed as in the eccentricity dependent case) for (the mass dependence of the precision is nearly independent of eccentricity) with (the results look very similar for ). The formula is accurate for a broad range of masses, but for some systems, such as Planet Hunters 3c/d, indicated in Cyan in Fig. 5, the discrepancy becomes large, % (the masses and eccentricity vectors for this plot differ from PH3c/d, but a plot made for the parameters of that pair of planets looks very similar).

### 3.3. Comparison with two-planet systems

We have carried out fits to systems with two interacting planets described in LXW12, and we have re-fit Planet Hunters (PH3) c/d. We have carried out N-body dynamical analyses using TTVFast in addition to fits with the analytic first-order formula in order to assess its utility in analyzing multi-planet systems.

Our first case is Kepler-18c/d, which was published by Cochran et al. (2011) and also analyzed by LXW12. We used the same transit times and uncertainties from the Cochran et al. (2011) paper to allow for direct comparison to their results; these transit times were also used in LXW12. We carried out a markov chain monte carlo (MCMC) analysis using an affine-invariant population approach (Goodman & Weare, 2010). We allow a multiplicative factor for the timing uncertainties on each planet and we place a uniform prior on the eccentricity of each planet (Ford, 2006). Figure 6 shows a comparison of the results of the MCMC analyses with N-body integration versus the analytic formulae with . For this system the mass ratios are and the eccentricities are of order (1-), while , which is a regime in which the first-order formula is accurate to 9% compared with N-body integration (see Figure 4; the cyan dot indicates the approximate location of the upper end of the eccentricities of Kepler-18c/d). For conversion of the mass ratios to planet masses, we assume that (we ignore the uncertainty on this stellar mass). The masses of the planets derived from N-body are and , while from the analytic formula are and . These are well within of one another, the markov chain posteriors show very similar distributions (Fig. 6). We note that these results differ from those reported by Cochran et al. (2011) who used N-body to estimate the transit times, found the best fit using Levenberg-Marquardt optimization, and estimated the uncertainties from the Hessian matrix at the best fit parameters rather than a full posterior analysis. Our results also differ from the estimates in LXW12, which only solve for a ‘nominal’ mass assuming zero eccentricity using the approximate near-resonant formula. We feel that our results should be superior to these prior results, and warrant a more extensive analysis with the full Kepler dataset as well as radial velocity measurements.

We next compared analyses of the Kepler-28 system, which was originally studied by Steffen et al. (2012), and included in the analysis of LXW12. We used the transit times published by Steffen et al. (2012), and followed the same procedure as Kepler-18. Figure 7 shows a comparison of the masses and eccentricity vectors for this system, which has a mean period ratio of , just wide of 2:3 period ratio, corresponding to . The masses are poorly constrained due to the degeneracy with eccentricity, which allows the eccentricity to wander to larger values. For conversion of the mass ratios to planet masses, we assume that (we ignore the uncertainty on this stellar mass). A comparison between the mass constraints from N-body and the analytic formula gives: versus , and versus . The 85% confidence value of is 0.098, while for is 0.076, while the longitudes of periastron within the posterior distribution are primarily anti-aligned; the location of these points is indicated with a magenta datapoint in Figure 4. Note that in the anti-aligned case the analytic formula is valid to larger eccentricities, and thus is adequate to describe this system.

As Figure 5 indicates, the Planet Hunters 3 (PH3) system, the outer two planets (c/d) analyzed in Deck & Agol (2015), has a large discrepancy due to the large mass of the outer planet. The excellent agreement with the chopping formula given in Deck & Agol (2015) is still imperfect; the figure in that paper was mistakenly produced with larger timing error bars than used in Schmitt et al. (2014) which caused the agreement to appear slightly better than the first-order formula indicates. Figure 8 shows the results of a comparison of N-body and analytic fits to the PH3 transit times given in Schmitt et al. (2014); the planet masses assume . The analytic formula gives a significant discrepancy due to the large mass of the outer planet and due to the proximity to the 1:2 period ratio.

To confirm that the large mass of PH3d led to this discrepancy, we took the best-fit parameters resulting from our N-body analysis of the real PH3 data, and reduced the masses of the outer two planets by a factor of 10. Using the outer two planets alone, we simulated transit times and added Gaussian noise at a level of that of the noise of the real data to maintain the same signal-to-noise ratio as the actual data. We then modeled this simulated data using TTVFast and the analytic formulae. We found that the agreement between the N-body and analytic analyses becomes excellent (Fig. 8).

### 3.4. Comparison with multi-planet systems

The two-body solution can be used for more than two bodies by addition of two-body TTV solutions for each pair of two planets (Lissauer et al., 2011):

 δti1=∑i2≠i1δti1,i2, (21)

where are the solutions from equation (12) for the th planet due to the th planet. The sum over for each pair of planets can be carried up to to give sufficient precision for that pair of planets that is smaller than the measured timing precision.

Our first system of study is Kepler-51 (Masuda, 2014), consisting of planets with period ratios close to 1:2:3. We used the transit times reported in Masuda (2014) to carry out dynamical models with N-body/TTVFast and with an analytic TTV signal given as the sum of the TTVs of the three adjacent pairs of planets. We included up to in the TTV signals. The results show excellent agreement; Figure 9 shows the measured transit-timing variations, as well as the best-fit N-body and analytic TTVs. The results of the MCMC analyses are compared in Figure 10 which shows the posterior distribution of masses and eccentricities measured with both analyses, assuming . For two of the planets the eccentricities are consistent with zero; in these cases the tail of the eccentricity histogram is heavier for the N-body than the for the analytic formula. The masses of the inner two planets from our N-body and analytic MCMC analyses agree well with the masses from the analysis in Masuda (2014), while the mass of the outer planet is small by .

## 4. Numerical implementation and speed

The primary computational burden of equation (12) lies in computing the Laplace coefficients. We use a series solution for these coefficients, which gives both speed and accuracy, using code shared by Jack Wisdom. The secondary computational burden is in computing the sine and cosine terms, which involves four angles, and thus requires eight evaluations. We carry out the computation of higher sines and cosines using trigonometric addition formulae, which means we only need to compute eight trigonometric functions at each transit time once; the rest are gotten from addition and multiplication of these.

In the cases that we run a Markov chain for a set of planets, the initial value of is known fairly well from the period ratio of the planets. In this case the Laplace coefficients and their derivatives needed for the solution can be Taylor expanded at the given by an initial fit to the transit times, and these coefficients can be stored for evaluation of the coefficients at slightly different values of encountered during the MCMC simulation. This approach would not work if is being varied over a grid (for example, in the case of searching for a perturbing planet with unknown period); however, computational efficiency can still be achieved by reusing the Laplace coefficients at different eccentricities (Nesvorný & Morbidelli, 2008).

We have coded the first-order formula, equation (12), in C, IDL, Python, and Julia (Bezanzon et al., 2012). We carried out a benchmark comparison of the Julia implementation of the formula with the C implementation of TTVFast, and we find that it is faster when the Laplace coefficients are approximated from a Taylor expansion. As TTVFast is about 20 times faster than TTVs computed with standard N-body integrators, this represents nearly four orders of magnitude in speed up, similar to that found by Nesvorný & Morbidelli (2008). Note that if integer period ratios are chosen, sometimes the denominators of and can become infinite, causing divergence; we expect that this will not be encountered in practice as the formulae only apply to non-resonant planets.

The code implementing these equations may be accessed at https://github.com/ericagol/TTVFaster.

## 5. Discussion and conclusions

In modeling transit timing variations, degeneracies and computational speed can each prohibit the accurate measurement of transiting planet masses and orbital properties, with their attendant uncertainties. The degeneracy due to aliasing near first order resonances (LXW12) can be broken with very high signal-to-noise due to the slight difference in the eccentricity dependence as a function of period ratio, as well the presence of perturbations at other frequencies (Deck & Agol, 2015). Here we have tried to improve the modeling of the terms which are linearly dependent upon eccentricity to provide a higher-fidelity analytic model to address both the degeneracy and the computation barriers.

To this end, we have presented a first-order solution in eccentricity and mass ratio to the plane-parallel, near-circular 3-body problem on timescales shorter than the secular timescale. This improves to first order in eccentricity the original solution given in Agol et al. (2005) which was derived to zeroth order in eccentricity and first order in mass ratio (this solution has also been given in different forms in Nesvorný & Vokrouhlický 2014 and Deck & Agol 2015). The expressions are accurate compared with numerical integration over a wide range of parameter space relevant to the hundreds of multi-transiting planetary systems being found at short orbital periods with Kepler (Rowe et al., 2014; Lissauer et al., 2014). We find that this expression is more accurate than the stripped-down near-resonant formula given in LXW12, although their formula has a simpler form which clearly highlights the mass-eccentricity degeneracy. The first-order eccentricity formulae also can be used to model more than two planets with linear combinations of 2-planet formulae, and works well for the system we tested here, Kepler-51.

We used an approach starting with the Newtonian equations of motion rather than the Hamiltonian, and compute the perturbed polar coordinates of the planets’ orbits; this approach is akin to solving dispersion relations of differential equations for mode and stability analysis (for example, the magneto-rotational instability is derived with this approach, Chandrasekhar 1961; Balbus & Hawley 1998). The unperturbed solution to these differential equations represents Keplerian motion expanded in eccentricity. The terms with various frequencies in the disturbing function give an inhomogeneous component to the solution, which cause TTVs to vary at frequencies which depend upon integer combinations of the orbital frequencies of the two planets. Since the answer obtained in the end is the same as in the Hamiltonian (for ) and canonical transformation approaches, this approach might be useful pedagogically for those more familiar with stability analyses. In addition, this approach might be useful for other problems, such as a stability analysis of a two-planet system or for carrying out the TTV computation to second-order in mass ratios . The latter is interesting as it would reveal how TTVs of a transiting planet may be used to measure the mass of that planet (and not just the mass of the perturbing planet).

We expect that these formulae will be used in carrying out initial fits to multi-transiting planets which show TTVs (Mazeh et al., 2013), in searching for companion perturbing planets to isolated planets showing TTVs, in characterizing multi-planet systems with TTVs to confirm and check for convergence of N-body MCMC analyses, in forecasting TTV amplitude for follow-up measurement, in estimating the optimum times for transit observation, and in making predictions for transit times to plan observations. It should be useful for estimating the densities, masses, and radii of the host stars and their exoplanets (Agol et al., 2005; Montet & Johnson, 2012; Hadden & Lithwick, 2014; Kipping et al., 2014; Jontof-Hutter et al., 2015), and in comparing the TTV solutions to radial velocity solutions for the masses and orbits of exoplanets. The analytic nature of our solution should be amenable to automatic differentiation (Fournier et al., 2012), which could speed up optimization based on gradient computation, and could also enable Hamiltonian/Hybrid Markov Chain Monte Carlo (Neal, 2011). When a large number of planets transit a star and each show evidence for dynamical interactions, the number of free parameters describing the system becomes large, and thus MCMC becomes prohibitively computationally expensive. The first-order analytic formula developed here can be used for modeling these systems if their masses and eccentricities are in the allowable range. As the formula is about 400 times faster than TTVFast, which is already about 20 times faster than Bulirsch-Stoer based integrators, the total speedup of about 8000 should make running chains long enough to converge more feasible, especially in tandem with parallel computation which can be easily adapted for population MCMC (Foreman-Mackey et al., 2013). The analytic formula also has the advantage of being able to pinpoint which features in the TTVs constrain the parameters of the system (LWX12). Since the TTVs of a planet display harmonics of the perturbing planet (Deck & Agol, 2015), the amplitudes and phases of each of these harmonics can be measured directly from the TTVs, and then these can be used to place individual constraints upon the masses and eccentricity vectors of the planets. The regions where the constraints overlap may reveal the consistency and uniqueness of the solution for the system parameters in some cases (Deck & Agol, 2015).

Although in principle TTVs allow for unique measurements of planet mass and eccentricities, degeneracies between these parameters are often found for systems with low signal to noise. In these cases, the eccentricities can become extremely large, indicating unstable orbits, as long as the masses are adjusted in a corresponding manner. We applied Hill stability to avoid this problem when using the analytic formulae; the full N-body computation avoids this issue naturally since large eccentricities introduce second-order (and higher) variations (that our calculation ignores) which prevent the high-eccentricity cases from fitting the data well. It may be possible to break some degeneracies with transit duration variations (Pál & Kocsis, 2008; Nesvorný et al., 2013), which can be computed with the same formalism we have described here, albeit in the plane-parallel limit.

The first-order formulae described here could be extended to higher order in eccentricity and/or mass ratio, albeit with much more computational effort. A slightly more accurate formula might be obtained by computing the longitudes from Kepler’s equation at the times of transit of each planet rather than using the first-order eccentricity formula (11), as well as using the exact formula for at the transit times in equation (1); this requires very little additional computational effort, but will still be only accurate to first order in eccentricity. The solutions for the perturbed polar coordinates, , can be derived in the same manner that we have derived the transit timing variations. These are needed for carrying out the higher-order perturbation solutions, and in turn could be used for modeling astrometric variations, radial velocity varations, and pulsar timing variations of host stars to account for the interactions of planets to first order in eccentricity.

EA acknowledges support from NASA grants NNX13AF20G, NNX13AF62G, and NASA Astrobiology Instituteâs Virtual Planetary Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. EA thanks the Kavli Institute for Theoretical Physics and the organizers of the “Dynamics and Evolution of Earth-like Planets” workshop where a portion of this work was completed; this manuscript is preprint number NSF-KITP-15-132. KD acknowledges support from the Joint Center for Planetary Astronomy fellowship. We thank Jack Wisdom for sharing laplace.c which computes Laplace coefficients and their derivatives with series summation, we thank Eric Ford for advice on implementation of the formula in Julia, and we thank Brett Morris and Ethan Kruse for advice on implementation of the formula in Python (requested by the referee).

APPENDIX A

## A new approach to TTVs

Here we give the detailed derivation of the first-order solution presented in Section 2.

### .1. TTVs from angular and radial variations

We start with the equations of motion in Murray & Dermott (MD), 6.10-6.11, which are in heliocentric coordinates. We use the index to label the planets, and denote the inner planet with and outer with . TTVs can be computed from the true longitudes, and , with respect to the star; hence heliocentric coordinates are ideal for transit timing computation. We convert the equations of motion to polar coordinates, and isolate the radial and longitudinal equations:

 r2i¨θi+2ri˙ri˙θi = ∂∂θi(Ui+Ri)=∂Ri∂θi=˙li, (1) ¨ri−ri˙θ2i = ∂∂ri(Ui+Ri). (2)

The first equation can be rewritten as the time derivative of specific angular momentum, , and without the disturbing force it expresses conservation of angular momentum, while the second equation includes centripetal acceleration in the radial direction as the second term on the left hand side. The term is the standard Keplerian potential, while is the disturbing function which reflects the gravitational potential energy of planet-planet interactions. Because the planet-planet interactions lead to only small perturbations of the base Keplerian orbits, we seek a solution to equations (1) of the form: and , where is the unperturbed Keplerian polar coordinates of the orbits, and are the small perturbations. We can then plug these solutions into the equations of motion (1) and expand in powers of , mass ratio, and eccentricity.

We normalize by , which is the semi-major axis of the unperturbed Keplerian orbit (we do not perturb , , or , so they are fixed at the mean, unperturbed values). Then , where is a dimensionless radial coordinate, so that and . The solution to the unperturbed orbits to first order in eccentricity is and where , is the complex eccentricity vector (LXW12), is the complex conjugate of , and is the real part.

### .2. Perturbed equations of motion

The perturbed equations become:

 r2i,Kδ¨θi = −2ri,K¨θi,Kδri−2˙ri,K˙θi,Kδri−2ri,K˙θi,Kδ˙ri−2ri,K˙ri,Kδ˙θi+∂Ri∂θi, (3) δ¨ri = 2ri,K˙θi,Kδ˙θi+˙θ2i,Kδri+2G(m⋆+mi)r3i,Kδri+∂Ri∂ri. (4)

In these equations, we have cancelled the terms for the unperturbed Keplerian on both sides of the equation. From hereon we will drop the ’’ subscript from the unperturbed Keplerian orbital elements. To first order in eccentricity, the TTVs of the th planet are:

 δti=−n−1i(δθi+2ϵiδθ(0)i), (5)

where is the perturbed solution to zeroth order in eccentricity. To lowest order in eccentricity, , so:

 ϵi = −12(eie\mathbbmi(λi−ϖi)+eie−\mathbbmi(λi−ϖi)) (6) ˙ϵi = −\mathbbmini2(eie\mathbbmi(λi−ϖi)−eie−\mathbbmi(λi−ϖi)). (7)

### .3. First order in eccentricity equations

Let the specific angular momentum equal . Since is a constant, . To first order in eccentricity, , where . Substituting these into the above equations, and expanding to first order in eccentricity, we find

 δ¨ϵi−2ni(1−ϵi)δ˙θi−n2i(3−10ϵi)δϵi = 1a2i∂Ri∂ϵi (8) (1+2ϵi)δ¨θi−2ni˙ϵiδϵi+2ni(1−ϵi)δ˙ϵi+2˙ϵiδ˙θi = 1a2i∂Ri∂θi. (9)

Note that the terms with or are first order in eccentricity; hence, the other quantities in these terms need to only be expanded to zeroth order in eccentricity. Denoting the zeroth-order solutions as and , we find the differential equations governing the transit timing solution:

 δ¨θi+2niδ˙ϵi = −2(ϵiδ¨θ(0)i+˙ϵiδ˙θ(0)i)+2ni(˙ϵiδϵ(0)i+ϵiδ˙ϵ(0)i)+1a2i∂Ri∂θi, (10) δ¨ϵi−3n2iδϵi−2niδ˙θi = −2niϵi(δ˙θ(0)i+5niδϵ(0)i)+1a2i∂Ri∂ϵi. (11)

We assume that both sides of these equations are complex, and the final solution is found from taking their real parts.

The quantity is the disturbing function, which for the inner planet can be broken into two pieces: (MD 6.44), where and . For the outer planet, there are also two pieces, , where (MD 6.45). As usual, .

### .4. Expansion of the disturbing functions

The expansion for is given in MD 6.66. We are considering the plane-parallel case, so , in which case we only need to include the term. Also, we would like a solution that is first order in eccentricity (of the unperturbed Keplerian orbit), so the term in brackets in MD 6.66 needs to be expanded to second order in (noting again that is first order in eccentricity):

 2∑l=01l!l∑k=0(lk)ϵk1ϵl−k2A0,j,k,l−k=A0,j,0,0+A0,j,1,0ϵ1+A0,j,0,1ϵ2+12A0,j,0,2ϵ22+A0,j,1,1ϵ1ϵ2+12A0,j,2,0ϵ21, (12)

where (MD 6.63) and is a Laplace coefficient (MD 6.67). We define to simplify the expressions below.

We rewrite in complex notation (in the plane-parallel limit, expanded to second order in ), giving:

 RD = RD,0+Re[∑j≥1(~Aj00+~Aj10ϵ1+~Aj01ϵ2+12~Aj02ϵ22+~Aj11ϵ1ϵ2+12~Aj20ϵ21)e\mathbbmij(θ1−θ2)], (13) RD,0 = 12(~A000+~A010ϵ1+~A001ϵ2+12~A002ϵ22+~A011ϵ1ϵ2+12~A020ϵ21). (14)

where we have used the fact that since . Likewise,

 RE=−Re[(1+ϵ1)(1+ϵ2)−2e\mathbbmi(θ1−θ2)]. (15)

Taking the derivative of and with respect to gives to first order in :

 1a21∂R1∂θ1=1a21δ˙l1=n21μ2αRe[∑j≥1(~Aj00+~Aj10ϵ1+~Aj01ϵ2−α(1+ϵ1−2ϵ2)δj1)\mathbbmije\mathbbmij(θ1−θ2)], (16)

where we have used and . The derivative of with respect to is:

 1a21∂R1∂ϵ1 = n21μ2αRe[12(~A010+~A011ϵ2+~A020ϵ1) (17) + ∑j≥1(~Aj10+~Aj11ϵ2+~Aj20ϵ1−α(1−2ϵ2)δj1)e\mathbbmij(θ1−θ2)]. (18)

For the outer planet,

 1a22∂R2∂θ2=1a22δ˙l2=