The following article isFree article

TRANSIT TIMING TO FIRST ORDER IN ECCENTRICITY

Eric Agol andKatherine Deck

Published 2016 February 18 • © 2016. The American Astronomical Society. All rights reserved.
The Astrophysical Journal,Volume 818,Number 2Citation Eric Agol and Katherine Deck 2016ApJ818 177DOI 10.3847/0004-637X/818/2/177

DownloadArticle PDF

Eric Agol

AFFILIATIONS

Astronomy Department, University of Washington, Seattle, WA 98195, USA; agol@uw.edu

Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA

NASA Astrobiology Institutes Virtual Planetary Laboratory, Seattle, WA 98195, USA

Katherine Deck

AFFILIATIONS

Division of Geological and Planetary Sciences, Caltech, 1200 E. California Blvd., Pasadena, CA, 91125 USA

Skip to each figure in the article

Skip to each table in the article

Skip to each data item in the article

What is article data?

Eric Agol

AFFILIATIONS

Astronomy Department, University of Washington, Seattle, WA 98195, USA; agol@uw.edu

Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA

NASA Astrobiology Institutes Virtual Planetary Laboratory, Seattle, WA 98195, USA

Katherine Deck

AFFILIATIONS

Division of Geological and Planetary Sciences, Caltech, 1200 E. California Blvd., Pasadena, CA, 91125 USA

Article metrics

2336 Total downloads
0 Video abstract views

Permissions

Get permission to re-use this article

Share this article

Dates

Received2015 September 4
Accepted2015 November 24
Published2016 February 18

Keywords

planets and satellites: detection;planets and satellites: dynamical evolution and stability;planets and satellites: fundamental parameters

Journal RSS

Create or edit your corridor alerts

What are corridors?

0004-637X/818/2/177

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 expensiveN-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.

Export citation and abstractBibTeX RIS

Previous article in issue

Next article in issue

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; Schneider2003,2004; Agol et al.2005; Holman & Murray2005; Heyl & Gladman2007; Nesvorný & Morbidelli2008; Fabrycky2010). Transit timing variations (TTVs) have been used to confirm that transit signals are in fact due to planets (Holman et al.2010; Fabrycky et al.2012; Ford et al.2012a,2012b; Steffen et al.2012,2013; Xie2013; 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 on planetary system properties. Analytic formulae based on perturbation theory can greatly speedup 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ý & Morbidelli2008; 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 & Lithwick2014; Limbach & Turner2014; Van Eylen & Albrecht2015), 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 that 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 & Agol2015), or statistical analysis of many systems (Wu & Lithwick2013; Hadden & Lithwick2014; Xie et al.2014). In addition, theLXW12 analysis is only approximate and breaks down for pairs farther from resonance (Deck & Agol2015). 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 & Agol2015); (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 Section2 (the full derivation is given in AppendixA). We then compare these with prior results, both analytic and numeric (Section3), including comparison of analytic and numeric analyses of specific systems. We discuss the numerical implementation and speed in Section4. We end with a discussion of the possible applications and future directions (Section5).

2. FIRST-ORDER SOLUTION:

Here we give a complete summary of the assumptions and variables used, as well as the solution to the first-order equations for readers who wish to simply use the results of this computation. The details of the derivation are given in AppendixA. 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 the inner planet,i = 1, and four for the outer planet,i = 2 (semimajor axis,a_i, mean longitude, ${\lambda }_{i}$ , eccentricity,e_i, and longitude of periastron, ${\varpi }_{i}$ ), plus the mass ratio of each planet to the star, ${m}_{1}/{m}_{\star },{m}_{2}/{m}_{\star }$ , wherem₁ is the mass of the inner planet,m₂ is the mass of the outer planet, and ${m}_{\star }$ is the mass of the star. For nearly circular planetary orbits, there are two small dimensionless parameters in the problem: ${\mu }_{i}={m}_{i}/{m}_{\star }$ ande_i. The usual procedure for computing TTVs 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 which 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, ${\theta }_{i}={\theta }_{i,{\rm{K}}}+\delta {\theta }_{i}$ , as a function of time, where ${\theta }_{i,{\rm{K}}}$ is the unperturbed Keplerian orbit and $\delta {\theta }_{i}$ is the perturbation of theith planet caused by its planet companion(s); and (6) compute the TTVs:

This is the approach taken by Agol et al. (2005), Nesvorný & Vokrouhlický (2014), andLXW12. A different approach employing Hamiltonian perturbation theory (Nesvorný & Morbidelli2008) 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 $\delta {\theta }_{i}$ , rather than the perturbed orbital elements. The derivation based on canonical transformations (Nesvorný & Morbidelli2008), 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 $({r}_{i},{\theta }_{i})$ . We then use Newton’s equations in terms of a disturbing function that 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 per planet (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 TTVs; (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 inLXW12) and by expanding the derivatives of the disturbing function in terms of orbital elements, which yields harmonic functions that are easy to integrate. The final answer is expressed as a sum over harmonics of the perturbing planet’s orbital frequency (Deck & Agol2015). 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 TTVs at a particular frequency.

The unperturbed orbital frequencies, ${n}_{i}=2\pi /{P}_{i}$ , are defined by ${n}_{i}^{2}={{Gm}}_{\star }/{a}_{i}^{3}$ . As usual, $\alpha ={a}_{1}/{a}_{2}\approx {({P}_{1}/{P}_{2})}^{2/3}$ . We define ${\tilde{A}}_{{jmn}}={a}_{1}^{m}{a}_{2}^{n+1}\displaystyle \frac{{\partial }^{m+n}}{\partial {a}_{1}^{m}\partial {a}_{2}^{n}}({a}_{2}^{-1}{b}_{1/2}^{(j)}(\alpha ))$ , where ${b}_{1/2}^{(j)}(\alpha )$ is the Laplace coefficient,

The difference in mean longitude of the planets is $\psi ={\lambda }_{1}-{\lambda }_{2}$ . Auxiliary dimensionless quantities are

The functions ${\tilde{A}}_{{jmn}}$ we use below are given by

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

if we assume that the orbital reference is along the line of sight, and thus ${\lambda }_{i}\approx 0$ at the times of transit.

The solutions for the inner planet (i = 1) and outer planet (i = 2) are given by

where the functions ${f}_{i,j}^{(\pm k)}$ are given by

whereγ,c₁, andc₂ andζ,d₁, andd₂ are given in Table1, and ${\delta }_{{ik}}$ is the Kronecker delta function. Note that the top signs in $\pm ,\;\mp$ correspond to $+k$ values, while the bottom correspond to $-k$ . The functions ${f}_{i,j}^{(\pm k)}$ are solely a function ofj, $\pm k$ , andα.

Table 1. Coefficients for $u(\gamma ,{c}_{1},{c}_{2})$ and ${v}_{\pm }\;(\zeta ,{d}_{1},{d}_{2})$ in First-order TTV Solution

i	$\pm k$	$\gamma [\zeta ]$	${c}_{1}[{d}_{1}]$	${c}_{2}[{d}_{2}]$
1	$0[\pm 1]$	${\beta }_{j}$	$\alpha j({\tilde{A}}_{j00}-\alpha {\delta }_{j1})$	$\alpha ({\tilde{A}}_{j10}-\alpha {\delta }_{j1})$
1	±1	${\beta }_{j}\pm 1$	$\alpha j\left(\pm j{\tilde{A}}_{j00}-\displaystyle \frac{1}{2}{\tilde{A}}_{j10}+\displaystyle \frac{1}{2}(1\mp 2)\alpha {\delta }_{j1}\right)$	$\alpha \left(\pm j{\tilde{A}}_{j10}-\displaystyle \frac{1}{2}{\tilde{A}}_{j20}\mp \alpha {\delta }_{j1}\right)$
1	±2	${\beta }_{j}\pm {\alpha }^{3/2}$	$\alpha j\left(\mp j{\tilde{A}}_{j00}-\displaystyle \frac{1}{2}{\tilde{A}}_{j01}-(1\mp 1)\alpha {\delta }_{j1}\right)$	$\alpha \left(\mp j{\tilde{A}}_{j10}-\displaystyle \frac{1}{2}{\tilde{A}}_{j11}-(1\mp 1)\alpha {\delta }_{j1}\right)$
2	$0[\pm 2]$	${\kappa }_{j}$	$-j({\tilde{A}}_{j00}-{\alpha }^{-2}{\delta }_{j1})$	${\tilde{A}}_{j01}-{\alpha }^{-2}{\delta }_{j1}$
2	±1	${\kappa }_{j}\pm {\alpha }^{-3/2}$	$-j\left(\pm j{\tilde{A}}_{j00}-\displaystyle \frac{1}{2}{\tilde{A}}_{j10}-(1\pm 1){\alpha }^{-2}{\delta }_{j1}\right)$	$\pm j{\tilde{A}}_{j01}-\displaystyle \frac{1}{2}{\tilde{A}}_{j11}-(1\pm 1){\alpha }^{-2}{\delta }_{j1}$
2	±2	${\kappa }_{j}\pm 1$	$-j\left(\mp j{\tilde{A}}_{j00}-\displaystyle \frac{1}{2}{\tilde{A}}_{j01}+\displaystyle \frac{1}{2}(1\pm 2){\alpha }^{-2}{\delta }_{j1}\right)$	$\mp j{\tilde{A}}_{j01}-\displaystyle \frac{1}{2}{\tilde{A}}_{j02}\pm {\alpha }^{-2}{\delta }_{j1}$

Note. The ${v}_{\pm }$ coefficients correspond to thek values in brackets [..].

Download table as: ASCII Typeset image

In practice, the sum overj from 1 to $\infty$ must be truncated at a finite value ofj_max. Typicallyj_max does not need to be chosen to be too large since the Laplace coefficients decline in amplitude withj (Deck & Agol2015). We recommend choosing aj_max large enough such that the resulting computation is converged.

As an example of using Table1, the coefficient ${f}_{1,j}^{(-1)}$ has $i=k=1$ , $\gamma ={\beta }_{j}-1$ , ${c}_{1}=\alpha j\left(-j{\tilde{A}}_{j00}-\displaystyle \frac{1}{2}{\tilde{A}}_{j10}+\displaystyle \frac{3}{2}\alpha {\delta }_{j1}\right)$ and ${c}_{2}=\alpha \left(-j{\tilde{A}}_{j10}-\displaystyle \frac{1}{2}{\tilde{A}}_{j20}+\alpha {\delta }_{j1}\right)$ , $\zeta ={\beta }_{j}$ , ${d}_{1}=\alpha j({\tilde{A}}_{j00}-\alpha {\delta }_{j1})$ , and ${d}_{2}=\alpha ({\tilde{A}}_{j10}-\alpha {\delta }_{j1})$ . Then, the coefficient is given by

3. COMPARISON WITH OTHER FORMULAE

The zeroth-order solution (in the limit ${e}_{1}={e}_{2}=0$ ) 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, ${f}_{i,j}^{(\pm k)}(\alpha )$ , as a function of period ratio, ${P}_{2}/{P}_{1}\approx {\alpha }^{-3/2}$ . The zeroth-order eccentricity coefficients are plotted in Deck & Agol (2015). The first-order coefficients can show three singularities for the terms with superscripts (−1) and (−2) near first-order resonance, second-order resonance, and $\alpha =1$ .

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Dimensionless coefficients, ${f}_{1,j}^{(\pm k)}$ for the inner planet. The dashed lines show where the coefficients are negative.
Download figure:
Standard image High-resolution image

**Figure 1.** Dimensionless coefficients, ${f}_{1,j}^{(\pm k)}$ for the inner planet. The dashed lines show where the coefficients are negative.
Download figure:
Standard image High-resolution image

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Dimensionless coefficients, ${f}_{2,j}^{(\pm k)}$ for the outer planet. The dashed lines show where the coefficients are negative.
Download figure:
Standard image High-resolution image

**Figure 2.** Dimensionless coefficients, ${f}_{2,j}^{(\pm k)}$ for the outer planet. The dashed lines show where the coefficients are negative.
Download figure:
Standard image High-resolution image

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 foru andv_± do not show the same dependence in the denominator as the expressions inLXW12; their expression just contains the resonant frequency, ${{jn}}_{1}-(j+1){n}_{2}$ , while ours contains additional frequencies. We carried out the partial fraction expansion ofu to isolate the denominator that matchesLXW12's expression, and we find that the expressions agree exactly with their expressions (A28) and (A29).

To compare our full expression withLXW12's, we have computed the eccentricity-dependentj = 2 (near 2:3) expression for the inner and outer planets (as this term is unaffected by the indirect terms). We rewrite the TTV formulae derived here and inLXW12 for theith planet perturbed by planetk as

where we have set ${\theta }_{i}={\lambda }_{i}=0$ at transit (this incurs some error, at ordere, but that is a second-order effect since we are comparing the TTV term linear ine). When written in this way, the amplitude and phase depend only on $\alpha ,{\varpi }_{1}$ , and ${\varpi }_{2}$ .

For the 3:2 resonance, theLXW12 resonant term depends on ${e}_{1}/{(2{n}_{1}-3{n}_{2})}^{2}$ and ${e}_{2}/{(2{n}_{1}-3{n}_{2})}^{2}$ , which both decline quickly away from resonance, and thus other terms that depend one₁ ande₂ make a more significant contribution farther from resonance; hence, our formulae agree close to resonance but diverge away from exact commensurability. Figure3 shows the fractional error in the amplitudeA and phaseϕ of the terms that are proportional to the eccentricities of the planets for ${\varpi }_{1}\approx {\varpi }_{2}=0.45$ rad. The error in theLXW12 expressions (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 & Agol2015). 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 ${e}^{{\mathbb{i}}(j+1)({\lambda }_{1}-{\lambda }_{2})}$ , while the first-order term scales as ${e}^{{\mathbb{i}}(j{\lambda }_{1}-(j+1){\lambda }_{2})}$ , when ${\mathbb{i}}=\sqrt{-1}$ . Since the mean longitude of the inner planet is nearly identical at each transit of the inner planet, the ${\lambda }_{1}$ term in the exponent is approximately constant, while the ${\lambda }_{2}$ 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.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Fractional error in the coefficients of the 3:2 resonant TTV expression given inLXW12 compared with thej = 2 terms of our analytic expression. Black indicates the dependence on thee₁ term, while red indicates the dependence on thee₂ term. The top panels show the fractional errors in the amplitudes of the inner (left) and outer (right) planets. The bottom panels show the fractional error on the phases.
Download figure:
Standard image High-resolution image

3.2. Comparison with N-body Integrations

We have carried out extensive integrations of three-body systems usingTTVFast (Deck et al.2014) and compared the results with the first-order analytic formulae (Equation (6)). Note that theTTVFast 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 inLXW12 and Deck & Agol (2015). The longitudes computed from theTTVFast code need to have $\pi /2$ subtracted to make the plots shown below. In addition, the orbital elements accepted byTTVFast 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.

3.2.1. Eccentricity and Period Ratio Dependence

Figure4 compares the precision of the analytic formula as a function of $\alpha ={({P}_{1}/{P}_{2})}^{2/3}$ and eccentricity of both planets, which are set to be equal, ${e}_{1}={e}_{2}$ . We have set ${\varpi }_{1}={\varpi }_{2}+\pi$ , which we found (approximately) maximizes the discrepancy of the analytic model compared with theN-body model, and ${\varpi }_{1}={\varpi }_{2}$ , 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 ${P}_{1}=30$ days, and we have integrated the system withTTVFast for 1600 days, about the duration of the initial Kepler mission, assuming plane-parallel orbits. For these tests we assume ${\mu }_{1}={\mu }_{2}={10}^{-5}$ , 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 theN-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 ${j}_{\mathrm{max}}=10$ . 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-averagedN-body simulation, while we computedα in the analytic formula from the ratio of the periods derived from the best-fit ephemerides.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Fractional precision of the analytic formula compared withTTVFast. Left: aligned longitudes of periastron ( ${\varpi }_{1}={\varpi }_{2}$ ); right: anti-aligned longitudes of periastron ( ${\varpi }_{1}={\varpi }_{2}+\pi$ ). The dotted lines indicate the 10% precision level. Cyan dots show the approximate position of Kepler-18b/c at the 85.15% posterior eccentricity value, while magenta is used for Kepler-28. The region in the upper right is Hill unstable; these models were not computed and default to 100% uncertainty in this plot. The green dashed lines show the locations ofj:j + 1 resonances, while the blue dashed lines showj:j + 2.
Download figure:
Standard image High-resolution image

**Figure 4.** Fractional precision of the analytic formula compared withTTVFast. Left: aligned longitudes of periastron ( ${\varpi }_{1}={\varpi }_{2}$ ); right: anti-aligned longitudes of periastron ( ${\varpi }_{1}={\varpi }_{2}+\pi$ ). The dotted lines indicate the 10% precision level. Cyan dots show the approximate position of Kepler-18b/c at the 85.15% posterior eccentricity value, while magenta is used for Kepler-28. The region in the upper right is Hill unstable; these models were not computed and default to 100% uncertainty in this plot. The green dashed lines show the locations ofj:j + 1 resonances, while the blue dashed lines showj:j + 2.
Download figure:
Standard image High-resolution image

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 theN-body integration. Figure4 shows that the formula works to better than 10% precision for a wide range of $\alpha \mbox{--}e$ parameter space. However, it fails near resonances, most significantly for thej:j + 1 resonances indicated in green andj:j + 2 in blue. For the outer planet, there are narrow regions near 1:j 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 TTVs plus the specific variations of the nonresonant terms are necessary (although this approach breaks down for large enough eccentricity).

3.2.2. Mass Dependence

We have carried out simulations for a range of masses, keeping ${m}_{1}={m}_{2}({\mu }_{1}={\mu }_{2})$ . We find that the fractional error of the analytic formula grows near resonances and near $\alpha =1$ as the mass increases with weaker dependence on eccentricity. For smallα, the formula works well up to ${m}_{1}={m}_{2}={10}^{-3}{m}_{\star }$ , while near 1:2 period ratio, for example, the error broadens around the resonance.

Figure5 shows the fractional error in the formula (computed as in the eccentricity-dependent case) for ${e}_{1}={e}_{2}=0.001$ (the mass dependence of the precision is nearly independent of eccentricity) with ${\varpi }_{1}={\varpi }_{2}+\pi$ (the results look very similar for ${\varpi }_{1}={\varpi }_{2}$ ). The formula is accurate for a broad range of masses, but for some systems, such as Planet Hunters 3c/d, indicated in cyan in Figure5, the discrepancy becomes large, $\approx 10$ % (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).

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Fractional precision of the analytic formula compared withTTVFast vs.α and mass ratio of the planets to the star ( ${m}_{1}={m}_{2}$ ). Top: inner planet; bottom: outer planet. The dotted line indicates the 10% precision level. Cyan dots show the approximate position of PH3 (although note that PH3 does not have equal masses; however, the plot is similar for parameters appropriate for PH3). The region in the upper right is Hill unstable; these models were not computed and default to 100% uncertainty in this plot. The green dashed lines show the locations ofj:j + 1 resonances, while the blue dashed lines showj:j + 2.
Download figure:
Standard image High-resolution image

**Figure 5.** Fractional precision of the analytic formula compared withTTVFast vs.α and mass ratio of the planets to the star ( ${m}_{1}={m}_{2}$ ). Top: inner planet; bottom: outer planet. The dotted line indicates the 10% precision level. Cyan dots show the approximate position of PH3 (although note that PH3 does not have equal masses; however, the plot is similar for parameters appropriate for PH3). The region in the upper right is Hill unstable; these models were not computed and default to 100% uncertainty in this plot. The green dashed lines show the locations ofj:j + 1 resonances, while the blue dashed lines showj:j + 2.
Download figure:
Standard image High-resolution image

3.3. Comparison with Two-planet Systems

We have carried out fits to systems with two interacting planets described inLXW12, and we have refit Planet Hunters (PH3) c/d. We have carried outN-body dynamical analyses usingTTVFast 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 byLXW12. 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 inLXW12. We carried out a Markov Chain Monte Carlo (MCMC) analysis using an affine-invariant population approach (Goodman & Weare2010). We allow a multiplicative factor for the timing uncertainties on each planet and place a uniform prior on the eccentricity of each planet (Ford2006). Figure6 shows a comparison of the results of the MCMC analyses withN-body integration versus the analytic formulae with ${j}_{\mathrm{max}}=5$ . For this system the mass ratios are $\approx 5\times {10}^{-5}$ and the eccentricities are of order 0–0.02 (1σ), while $\alpha \approx 0.64$ , which is a regime in which the first-order formula is accurate to <9% compared withN-body integration (see Figure4; 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 ${m}_{\star }=0.972\quad {M}_{\odot }$ (we ignore the uncertainty on this stellar mass). The masses of the planets derived fromN-body are ${m}_{1}(\mathrm{nbody})={14.7}_{-6.7}^{+5.4}{M}_{\oplus }$ and ${m}_{2}(\mathrm{nbody})={14.3}_{-4.1}^{+2.3}{M}_{\oplus }$ , while from the analytic formula they are ${m}_{1}(\mathrm{analytic})={13.2}_{-7.0}^{+5.4}{M}_{\oplus }$ and ${m}_{2}(\mathrm{analytic})={13.6}_{-4.9}^{+2.4}{M}_{\oplus }$ . These are well within $1\sigma$ of one another; the Markov chain posteriors show very similar distributions (Figure6). We note that these results differ from those reported by Cochran et al. (2011), who usedN-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 inLXW12, 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 data set, as well as radial velocity measurements.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Comparison of numerical and analytic analysis of the transit times of Kepler-18. Left: histogram of the masses of each planet. Right: comparison of the 68% confidence distribution of the eccentricity vectors of the two planets (18c, solid; 18d, dashed). In each case black indicates the results from theN-body analysis, while blue indicates the results of the analytic formula.
Download figure:
Standard image High-resolution image

We next compared analyses of the Kepler-28 system, which was originally studied by Steffen et al. (2012) and included in the analysis ofLXW12. We used the transit times published by Steffen et al. (2012) and followed the same procedure as for Kepler-18. Figure7 shows a comparison of the masses and eccentricity vectors for this system, which has a mean period ratio of $\langle {P}_{2}/{P}_{1}\rangle =1.52$ , just wide of the 2:3 period ratio, corresponding to $\alpha =0.7563$ . The masses are poorly constrained owing 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 ${m}_{\star }=0.89{M}_{\odot }$ (we ignore the uncertainty on this stellar mass). A comparison between the mass constraints fromN-body and the analytic formula gives ${m}_{1}(\mathrm{nbody})={3.8}_{-2.3}^{+4.6}{M}_{\oplus }$ versus ${m}_{1}(\mathrm{analytic})={3.1}_{-1.7}^{+3.4}{M}_{\oplus }$ and ${m}_{2}(\mathrm{nbody})={5.1}_{-3.0}^{+5.9}{M}_{\oplus }$ versus ${m}_{2}(\mathrm{analytic})={4.1}_{-2.3}^{+4.6}{M}_{\oplus }$ . The 85% confidence value ofe₁ is 0.098, while fore₂ it 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 data point in Figure4. Note that in the anti-alignedϖ case the analytic formula is valid to larger eccentricities and thus is adequate to describe this system.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Comparison of numerical and analytic analysis of the transit times of Kepler-28. Left: histogram of the masses of each planet. Right: comparison of the 68% confidence distribution of the eccentricity vectors of the two planets (28b, solid; 28c, dashed). In each case black indicates the results from theN-body analysis, while blue indicates the results of the analytic formula.
Download figure:
Standard image High-resolution image

As Figure5 indicates, the Planet Hunters 3 (PH3) system, the outer two planets (c/d) analyzed in Deck & Agol (2015), has a large discrepancy owing 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. Figure8 shows the results of a comparison ofN-body and analytic fits to the PH3 transit times given in Schmitt et al. (2014); the planet masses assume ${m}_{\star }=1\quad {M}_{\odot }$ . 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.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Comparison of numerical and analytic analysis of the transit times of Planet Hunters 3c/d. Left: 68% confidence contour of the masses of both planets. Right: comparison of the 68% confidence limits with the mass of the both planets reduced by a factor of 10.
Download figure:
Standard image High-resolution image

To confirm that the large mass of PH3d led to this discrepancy, we took the best-fit parameters resulting from ourN-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 1/10 that of the noise of the real data to maintain the same signal-to-noise ratio as the actual data. We then modeled these simulated data usingTTVFast and the analytic formulae. We found that the agreement between theN-body and analytic analyses becomes excellent (Figure8).

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):

where $\delta {t}_{{i}_{1},{i}_{2}}$ are the solutions from Equation (6) for thei₁th planet due to thei₂th planet. The sum overj for each pair of planets can be carried up toj_max 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 (Masuda2014), 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 withN-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 toj = 6 in the TTV signals. The results show excellent agreement; Figure9 shows the measured transit-timing variations, as well as the best-fitN-body and analytic TTVs. The results of the MCMC analyses are compared in Figure10, which shows the posterior distribution of masses and eccentricities measured with both analyses, assuming ${m}_{\star }=1\quad {M}_{\odot }$ . For two of the planets the eccentricities are consistent with zero; in these cases the tail of the eccentricity histogram is heavier for theN-body than for the analytic formula. The masses of the inner two planets from ourN-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 $\approx 1\sigma$ .

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** TTVs from Masuda (2014) for Kepler-51 (red; 0, 1, 2 stand for 51b,c, 620.02), compared with the best-fitN-body model computed withTTVFast (black) and the best-fit two-planet, first-order eccentricity formula summed over pairs of planets (blue).
Download figure:
Standard image High-resolution image

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Comparison of numerical and analytic analyses of the transit times of Kepler-51 (0, 1, 2 stand for 51b, c, 620.02). Left: histograms of the masses of each planet. Right: histograms of the eccentricities of the planets.
Download figure:
Standard image High-resolution image

4. NUMERICAL IMPLEMENTATION AND SPEED

The primary computational burden of Equation (6) 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-j 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ý & Morbidelli2008).

We have coded the first-order formula, Equation (6), 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 ofTTVFast, and we find that it is $400\times$ faster when the Laplace coefficients are approximated from a Taylor expansion. AsTTVFast is about 20 times faster than TTVs computed with standardN-body integrators, this represents nearly four orders of magnitude in speedup, similar to that found by Nesvorný & Morbidelli (2008). Note that if integer period ratios are chosen, sometimes the denominators ofu and ${v}_{\pm }$ can become infinite, causing divergence; we expect that this will not be encountered in practice as the formulae only apply to nonresonant planets.

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

5. DISCUSSION AND CONCLUSIONS

In modeling TTVs, 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 ratio owing to the slight difference in the eccentricity dependence as a function of period ratio, as well the presence of perturbations at other frequencies (Deck & Agol2015). Here we have tried to improve the modeling of the terms that are linearly dependent on 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 three-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; Deck & Agol2015). 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 (Lissauer et al.2014; Rowe et al.2014). We find that this expression is more accurate than the stripped-down near-resonant formula given inLXW12, although their formula has a simpler form that 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 two-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 magnetorotational instability is derived with this approach; Chandrasekhar1961; Balbus & Hawley1998). 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 causes TTVs to vary at frequencies that depend on 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 ${ \mathcal O }({e}^{0})$ ) 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 ${\mu }_{i}={m}_{i}/{m}_{\star }$ . 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 that 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 ofN-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 & Johnson2012; Hadden & Lithwick2014; 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 MCMC (Neal2011). When a large number of planets transit a star and each shows 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 thanTTVFast, 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 & Agol2015), 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 on 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 & Agol2015).

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 ratio. 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 fullN-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 & Kocsis2008; 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 (5), as well as using the exact formula for ${\dot{\theta }}_{i}$ at the transit times in Equation (1); this requires very little additional computational effort, but $\delta {\theta }_{i}$ will still be only accurate to first order in eccentricity. The solutions for the perturbed polar coordinates, $(\delta {r}_{i},\delta {\theta }_{i})$ , can be derived in the same manner that we have derived the TTVs. 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.

E.A. acknowledges support from NASA grants NNX13AF20G, NNX13AF62G, and NASA Astrobiology Institutes 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. E.A. 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. K.D. 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 Section2.

A.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 indexi to label the planets, and denote the inner planet withi = 1 and outer withi = 2. TTVs can be computed from the true longitudes, ${\theta }_{1}$ and ${\theta }_{2}$ , 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:

The first equation can be rewritten as the time derivative of specific angular momentum,l_i, 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 ${U}_{i}=G({m}_{\star }+{m}_{i})/{r}_{i}$ is the standard Keplerian potential, while ${{ \mathcal R }}_{i}$ is the disturbing function that 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 (11) of the form ${r}_{i}={r}_{i,{\rm{K}}}+\delta {r}_{i}$ and ${\theta }_{i}={\theta }_{i,{\rm{K}}}+\delta {\theta }_{i}$ , where $({r}_{i,{\rm{K}}},{\theta }_{i,{\rm{K}}})$ is the unperturbed Keplerian polar coordinates of the orbits and $(\delta {r}_{i},\delta {\theta }_{i})$ are the small perturbations. We can then plug these solutions into the equations of motion (11) and expand in powers of $\delta {r}_{i},\delta {\theta }_{i}$ , mass ratio, and eccentricity.

We normalizer_i bya_i, which is the semimajor axis of the unperturbed Keplerian orbit (we do not perturba_i, ${\varpi }_{i}$ , ore_i, so they are fixed at the mean, unperturbed values). Then ${r}_{i}/{a}_{i}=1+{\epsilon }_{i}$ , where ${\epsilon }_{i}$ is a dimensionless radial coordinate, so that ${\dot{r}}_{i}/{a}_{i}={\dot{\epsilon }}_{i}$ and $\delta {r}_{i}/{a}_{i}=\delta {\epsilon }_{i}$ . The solution to the unperturbed orbits to first order in eccentricity is ${\theta }_{i,{\rm{K}}}={\lambda }_{i}-\mathrm{Re}[2{\mathbb{i}}{z}_{i}^{*}{e}^{{\mathbb{i}}{\lambda }_{i}}]={\lambda }_{i}+2{e}_{i}\mathrm{sin}({\lambda }_{i}-{\varpi }_{i})$ and ${\epsilon }_{i,{\rm{K}}}=-{\rm{R}}{\rm{e}}[{z}_{i}^{*}{e}^{{\mathbb{i}}{\lambda }_{i}}]$ , where ${\mathbb{i}}=\sqrt{-1}$ , ${z}_{i}={e}_{i}{e}^{{\mathbb{i}}{\varpi }_{i}}={e}_{i}(\mathrm{cos}{\varpi }_{i}+{\mathbb{i}}\mathrm{sin}{\varpi }_{i})$ is the complex eccentricity vector (LXW12),z* is the complex conjugate ofz, and $\mathrm{Re}[.]$ is the real part.

A.2. Perturbed Equations of Motion

The perturbed equations become

In these equations, we have cancelled the terms for the unperturbed Keplerian on both sides of the equation. From here on we will drop the “K” subscript from the unperturbed Keplerian orbital elements. To first order in eccentricity, the TTVs of theith planet are

where $\delta {\theta }_{i}^{(0)}$ is the perturbed solution to zeroth order in eccentricity. To lowest order in eccentricity, ${\epsilon }_{i}={-e}_{i}\mathrm{cos}({\lambda }_{i}-{\varpi }_{i})$ , so

A.3. First Order in Eccentricity Equations

Let the specific angular momentum equal ${l}_{i}={r}_{i}^{2}{\dot{\theta }}_{i}$ . Sincel_i is a constant, ${\ddot{\theta }}_{i}=-2{l}_{i}{\dot{r}}_{i}/{r}_{i}^{3}$ . To first order in eccentricity, ${l}_{i}={n}_{i}{a}_{i}^{2}$ , where ${n}_{i}=2\pi /{P}_{i}$ . Substituting these into the above equations, and expanding to first order in eccentricity, we find

Note that the terms with ${\epsilon }_{i}$ or ${\dot{\epsilon }}_{i}$ 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 $\delta {\epsilon }_{i}^{(0)}$ and $\delta {\theta }_{i}^{(0)}$ , we find the differential equations governing the transit timing solution:

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

The quantity ${{ \mathcal R }}_{i}$ is the disturbing function, which for the inner planet can be broken into two pieces: ${{ \mathcal R }}_{1}=\frac{{{Gm}}_{2}}{{a}_{2}}({{ \mathcal R }}_{D}+\alpha {{ \mathcal R }}_{E})$ (MD 6.44), where ${{ \mathcal R }}_{D}={a}_{2}/| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}|$ and ${{ \mathcal R }}_{E}=-({r}_{1}/{a}_{1}){({a}_{2}/{r}_{2})}^{2}\mathrm{cos}({\theta }_{1}-{\theta }_{2})$ . For the outer planet, there are also two pieces, ${{ \mathcal R }}_{2}=\left({{Gm}}_{1}/{a}_{2}\right)({{ \mathcal R }}_{D}+{\alpha }^{-2}{{ \mathcal R }}_{I})$ , where ${{ \mathcal R }}_{I}=-(1+{\epsilon }_{2}){(1+{\epsilon }_{1})}^{-2}\mathrm{cos}({\theta }_{1}-{\theta }_{2})$ (MD 6.45). As usual, $\alpha ={a}_{1}/{a}_{2}\approx {({P}_{1}/{P}_{2})}^{2/3}$ .

A.4. Expansion of the Disturbing Functions

The expansion for ${{ \mathcal R }}_{D}$ is given in MD 6.66. We are considering the plane-parallel case, so ${\boldsymbol{\Psi }}=\mathrm{cos}\psi -\mathrm{cos}({\theta }_{1}-{\theta }_{2})=0$ , in which case we only need to include the $\propto {{\boldsymbol{\Psi }}}^{0}$ 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 ${\epsilon }_{i}=\displaystyle \frac{{r}_{i}}{{a}_{i}}-1$ (noting again that ${\epsilon }_{i}$ is first order in eccentricity):

where ${A}_{0,j,m,n}={a}_{1}^{m}{a}_{2}^{n}\frac{{\partial }^{m+n}}{\partial {a}_{1}^{m}\partial {a}_{2}^{n}}({a}_{2}^{-1}{b}_{1/2}^{(j)}(\alpha ))$ (MD 6.63) and ${b}_{s}^{(j)}$ is a Laplace coefficient (MD 6.67). We define ${\tilde{A}}_{{jmn}}={a}_{2}{A}_{0,j,m,n}$ to simplify the expressions below.

We rewrite ${{ \mathcal R }}_{D}$ in complex notation (in the plane-parallel limit, expanded to second order in ${\epsilon }_{i}$ ), giving

where we have used the fact that ${A}_{0,j,k,l}={A}_{0,-j,k,l}$ since ${b}_{s}^{(j)}={b}_{s}^{(-j)}$ . Likewise,

Taking the derivative of ${{ \mathcal R }}_{D}$ and ${{ \mathcal R }}_{E}$ with respect to ${\theta }_{1}$ gives to first order in ${\epsilon }_{1}$

where we have used ${n}_{i}^{2}={{Gm}}_{\star }/{a}_{i}^{3}\quad =$ and ${\mu }_{i}={m}_{i}/{m}_{\star }$ . The derivative of ${{ \mathcal R }}_{1}$ with respect to ${\epsilon }_{1}$ is

For the outer planet,

The derivative of ${{ \mathcal R }}_{2}$ with respect to ${\epsilon }_{2}$ is

The angles and radii in the derivatives of the disturbing function can be expanded to first order in eccentricity, yielding

where $\psi ={\lambda }_{1}-{\lambda }_{2}$ . We have also combined thej = 0 eccentricity terms since $\mathrm{Re}({z}_{i}^{*}{e}^{{\mathbb{i}}{\lambda }_{i}})=\mathrm{Re}({z}_{i}{e}^{-{\mathbb{i}}{\lambda }_{i}})=\frac{1}{2}({z}_{i}^{*}{e}^{{\mathbb{i}}{\lambda }_{i}}+{z}_{i}{e}^{-{\mathbb{i}}{\lambda }_{i}})$ .

For the outer planet,

A.5. Trial Solution

The derivatives of the complex disturbing function contain terms that are proportional to ${e}^{{\mathbb{i}}j\psi }$ , ${e}^{{\mathbb{i}}j\psi }{e}_{i}{e}^{\pm {\mathbb{i}}({\lambda }_{i}-{\varpi }_{i})}$ for $j\geqslant 0$ . We treat these terms as harmonic driving terms and solve the inhomogeneous partial differential equations term by term. We expand the complex solutions for $\delta {\theta }_{1}$ and $\delta {\epsilon }_{1}$ as trial solutions:

and

We also define the solutions to zeroth order in eccentricity as

Then, the (real) TTVs are equal to

and

A.6. Inner Planet Coefficients

Substituting the $(j,0,\pm k)$ trial solutions into the above differential equations, to zeroth order in eccentricity we find for the inner planet

where ${\beta }_{j}=j({n}_{1}-{n}_{2})/{n}_{1}\quad =\quad j(1-{\alpha }^{3/2})$ and we have divided the equations by ${n}_{1}^{2}$ to make them dimensionless.

Similarly, we can write down the equations for the coefficients to first order ine₁. Note that to solve for $\delta {t}_{1,j}^{(\pm 1)}$ , we need to compute $\delta {\theta }_{1,j}^{(\pm 1)}-\delta {\theta }_{1,j}^{(0)}$ . Hence, we can subtract the matrix on the left times the vector $\{\delta {\theta }_{1,j}^{(0)},0\}$ from both sides of the equation. This may be rewritten in dimensionless form as

The solutions of Equation (31) for $(\delta {\theta }_{1,j}^{(0)},\delta {\epsilon }_{1,j}^{(0)})$ may be plugged into the first term on the right-hand side of this equation to solve for the first-order eccentricity terms.

To first order ine₂ in dimensionless form (for $j\geqslant 1$ ),

where ${\eta }_{\pm }={\beta }_{j}\pm {\alpha }^{3/2}$ $=\;j({n}_{1}-{n}_{2})/{n}_{1}\pm {n}_{2}/{n}_{1}$ $=\;j(1-{\alpha }^{3/2})\pm {\alpha }^{3/2}$ and ${\tilde{A}}_{j11}=-(2{\tilde{A}}_{j10}+{\tilde{A}}_{j20})\;=$ $-2\alpha \partial {b}_{1/2}^{(j)}$ $/\partial \alpha -{\alpha }^{2}{\partial }^{2}{b}_{1/2}^{(j)}/\partial {\alpha }^{2}$ . Note that forj = 1 the ${\eta }_{+}$ term becomes ${\eta }_{+}=1$ . The frequency dependence of this term is at the Keplerian frequency of the inner planet,n₁, and the determinant of the left-hand matrix becomes zero as the second row of the matrix is equal to $2{\mathbb{i}}$ times the first row. This singularity occurs as a result of the fact that the equations become those of a resonantly driven oscillator, which means that the amplitude grows linearly with time (or, equivalently on short timescales, the Keplerian frequency is shifted). This term is not relevant for transit-timing analyses as it occurs at the frequency of the transiting planet and grows on the secular timescale; consequently, we will drop this term for now and discuss below in Appendix B.

Forj = 0,

A.7. Outer Planet Coefficients

Defining ${\kappa }_{j}=j({n}_{1}-{n}_{2})/{n}_{2}={\alpha }^{-3/2}{\beta }_{j}$ , dividing this equation by ${n}_{2}^{2}$ gives the dimensionless form of

and for the equations to first order ine₁ for $j\geqslant 1$ ,

where ${\xi }_{\pm }={\kappa }_{j}\pm {\alpha }^{-3/2}$ , while forj = 0,

and for first order ine₂,

The functionu originates from the inversion of the matrices on the left-hand side of Equations (31)–(38), which each have the form of

whereγ is the dimensionless frequency in units of the orbital frequency of the transiting planet. Since TTVs only depend on $\delta \theta$ , only the first term in the inverse of this matrix times the right-hand side of the equation gives the functionu. These terms are driven by the disturbing function only.

The function ${v}_{\pm }$ results from the driving terms caused by appearance of the zeroth-order eccentricity equation on the right-hand side of the linearized Equations (16). These can be expressed as the inverse of the zeroth-order matrix times the coefficients of the disturbing function driving the TTVs at zeroth order; the first term on the right-hand sides of Equations (32) and (36) times the inverse of the matrix on the left-hand side yields the functions ${v}_{\pm }$ .

The expressions foru and ${v}_{\pm }$ result from the coefficients in Equations (31)–(38), which can be found by inverting the matrices. Note that the coefficients of $\delta {\theta }_{i}$ are imaginary, while $\delta {\epsilon }_{i}$ are real; thus, $\delta {\theta }_{i}$ , and hence $\delta {t}_{i}$ , will always have a sine dependence, while $\delta {\epsilon }_{i}$ will always have a cosine dependence.

As with the inner planet, forj = 1 the ${\xi }_{-}$ term becomes ${\xi }_{-}=1$ . The frequency dependence of this term is at the Keplerian frequency of the outer planet,n₂, and hence the determinant of the left-hand matrix becomes zero as the second row of the matrix is equal to $2{\mathbb{i}}$ times the first row (as occurs for the inner planet). We will drop this term for now and discuss below in Appendix B.

APPENDIX B: SECULAR TERMS

In the foregoing analysis we neglected the presence of secular terms, which in the disturbing function appear at zero frequency, as well as at the Keplerian frequency of the planet that is being perturbed. These terms cause corrections of ${ \mathcal O }({\mu }^{1})$ to the ephemeris of the planet and thus cause an error of ${ \mathcal O }({\mu }^{1})$ to the computation ofα from the best-fit mean period. The correction toα affects the coefficients of the TTVs at order ${ \mathcal O }({\mu }^{2})$ , and so it can be neglected for the purposes of the first order in eccentricity transit timing solution. However, the solution we present here may have other applications, such as for radial-velocity planets or astrometric motion, which are not aliased at the orbital frequency of the planets, and so these secular terms enter at the ${ \mathcal O }({\mu }^{1})$ level. In this appendix we compute these secular terms to first order inμ ande.

In the equations of motion for the inner planet, to include the secular and Keplerian frequency terms, we will use angular momentum,l₁ in lieu of angle ${\theta }_{1}$ . The equations of motion become

where we have made the substitution ${r}_{1}={a}_{1}(1+{\epsilon }_{1})$ into Equations (11) and we have divided bya₁. In Equation (24) we keep only the secular terms and terms at frequencyn₁, giving

for the inner planet, and similarly for the outer planet

where we have taken the complex conjugate since this does not change the real component.

The solution for the inner planet’s angular momentum to first order in eccentricity is

This can be substituted into the equation for ${\epsilon }_{1}$ , keeping terms of order eccentricity, to obtain

Note that we have chosen the constant of integration inl₁ to cancel the constant term in the disturbing function derivative that appears in the equation for ${\epsilon }_{1}$ so that the ${\epsilon }_{1}$ does not have an offset. This is because we prefer to specify the value of the semimajor axis in the initial conditions.

The solution to this equation is

Note that this solution grows in amplitude linearly with time; however, the growth is slow, occuring on the secular timescale times the inverse of the eccentricity.

With these solutions in hand, we can solve for ${\dot{\theta }}_{1}={l}_{1}/{r}_{1}^{2}\approx {l}_{1}{a}_{1}^{-2}(1-2{\epsilon }_{1})+{ \mathcal O }({e}_{1}^{2})$ . We then integrate this with respect to time, giving

A similar solution can be derived for the outer planet, with

As before, this can be substituted into the equation for ${\epsilon }_{2}$ :

This equation has solution

Substituting this into the relation ${\dot{\theta }}_{2}={l}_{2}/{r}_{2}^{2}\approx {l}_{2}{a}_{2}^{-2}(1-2{\epsilon }_{2})+{ \mathcal O }({e}_{2}^{2})$ and integrating ${\dot{\theta }}_{2}$ with respect to time yields

We have verified these solutions by plugging them back into the differential equations, both analytically and numerically. The solutions for ${\theta }_{i,\mathrm{sec}}$ can be transformed to timing variations with $\delta {t}_{i,\mathrm{sec}}={-\dot{\theta }}_{i,\mathrm{sec}}^{-1}({\theta }_{i,\mathrm{sec}}-{\lambda }_{i})$ .

Please wait… references are loading.

10.3847/0004-637X/818/2/177

Movatterモバイル変換

TRANSIT TIMING TO FIRST ORDER IN ECCENTRICITY

Article metrics

Permissions

Share this article

Dates

ABSTRACT

1. INTRODUCTION

2. FIRST-ORDER SOLUTION:

3. COMPARISON WITH OTHER FORMULAE

3.1. Comparison with First-order Resonant Equations

3.2. Comparison with N-body Integrations

3.2.1. Eccentricity and Period Ratio Dependence

3.2.2. Mass Dependence

3.3. Comparison with Two-planet Systems

3.4. Comparison with Multi-planet Systems

4. NUMERICAL IMPLEMENTATION AND SPEED

5. DISCUSSION AND CONCLUSIONS

APPENDIX A: A NEW APPROACH TO TTVs

A.1. TTVs from Angular and Radial Variations

A.2. Perturbed Equations of Motion

A.3. First Order in Eccentricity Equations

A.4. Expansion of the Disturbing Functions

A.5. Trial Solution

A.6. Inner Planet Coefficients

A.7. Outer Planet Coefficients

APPENDIX B: SECULAR TERMS