Planet masses with TTV

# Measurement of planet masses with transit timing variations due to synodic “chopping” effects

## Abstract

Gravitational interactions between planets in transiting exoplanetary systems lead to variations in the times of transit that are diagnostic of the planetary masses and the dynamical state of the system. Here we show that synodic “chopping” contributions to these transit timing variations (TTVs) can be used to uniquely measure the masses of planets without full dynamical analyses involving direct integration of the equations of motion. We present simple analytic formulae for the chopping signal, which are valid (generally error) for modest eccentricities . Importantly, these formulae primarily depend on the mass of the perturbing planet, and therefore the chopping signal can be used to break the mass/free-eccentricity degeneracy which can appear for systems near first order mean motion resonances. Using a harmonic analysis, we apply these TTV formulae to a number of Kepler systems which had been previously analyzed with full dynamical analyses. We show that when chopping is measured, the masses of both planets can be determined uniquely, in agreement with previous results, but without the need for numerical orbit integrations. This demonstrates how mass measurements from TTVs may primarily arise from an observable chopping signal. The formula for chopping can also be used to predict the number of transits and timing precision required for future observations, such as those made by TESS or PLATO, in order to infer planetary masses through analysis of TTVs.

planetary systems
1

## 1. Introduction

In a multi-planet system, mutual gravitational interactions between planets lead to deviations from Keplerian orbits. In particular, the instantaneous orbital periods are no longer constant, which in turn implies that transiting planets in multi-planet systems will not transit at a fixed, constant rate. The detection of these changes in the transit rate, or ‘transit-timing variations’ (TTVs), was initially recognized as a way to infer the presence of non-transiting planets in systems with at least one other transiting planet (Schneider, 2003; Agol et al., 2005; Holman & Murray, 2005; Miralda-Escudé, 2002). TTVs have since been used to confirm that a transit light curve signal is due to a planetary transit (e.g. Holman et al., 2010), to constrain planetary orbital elements and measure planetary masses using photometry alone (e.g. Carter et al., 2012), and to detect and characterize non-transiting planets (e.g. Ballard et al., 2011; Nesvorný et al., 2012).

TTV data are most commonly analyzed through inversion, a process through which observed transit times are fit using a model of gravitationally interacting planets in order to determine the system parameters, including planetary masses relative to the mass of the star, as well as orbital elements.2 Transiting exoplanets generally have well constrained radii, so measurements of their masses yields information regarding densities, bulk compositions and gravities. This in turn can be used to identify promising targets for atmospheric characterization and to constrain planetary formation and dynamical evolution (e.g. Hansen & Murray, 2013). Of the orbital elements, constraints on the planetary eccentricities in particular are necessary to understand the importance of interaction with the protoplanetary disk, interactions with remnant planetesimals, and tidal dissipation (e.g. Lithwick & Wu, 2012; Hansen & Murray, 2013; Batygin & Morbidelli, 2013; Hansen & Murray, 2014; Mahajan & Wu, 2014).

However, the TTV inversion problem is often complicated by strong nonlinear correlations between parameters in a large dimensional space, and as a result precise planetary mass and orbit measurements can be difficult to make. Many of the Kepler multi-planet systems with partially characterized planetary orbits and masses are those near first order mean motion resonances, a configuration in which the period ratio of two planets is close to , where is an integer greater than unity, is the period of the inner planet, and is the period of the outer planet. Indeed, for nearly circular orbits and given planet-to-star mass ratios, TTVs are largest in amplitude near first-order mean motion resonances (e.g. Agol et al., 2005; Holman & Murray, 2005). If the planets are near to, but not in resonance, then the planets show sinusoidal variations with a period equal to the ‘super-period’,

 Pj=1|jR/P2−(jR−1)/P1| (1)

(Agol et al., 2005; Lithwick et al., 2012). However, the amplitude of this TTV signal depends on both the mass of the perturbing planet and the eccentricity vectors of both planets (Lithwick et al., 2012). This degeneracy can be broken statistically with analyses of a large number of planetary systems (Hadden & Lithwick, 2013) or for systems with very precisely measured transit times; however in practice it inhibits the measurement of the masses and limits our knowledge of the eccentricities of individual planetary systems.

In spite of this mass-eccentricity degeneracy (and others3), it has been possible in some cases to precisely measure the masses of planets using TTVs (e.g. Carter et al., 2012; Nesvorný et al., 2013; Masuda, 2014; Dreizler & Ofir, 2014a; Nesvorny et al., 2014). The successful mass measurements in these systems is due to the fact that an additional, independent periodic component of the TTVs, with a timescale other than the super-period, was resolved. Other components of TTVs have amplitudes that depend on the orbital parameters and masses in different ways, and so the measurement of secondary components leads to additional, independent constraints on orbital parameters. In particular, the so-called short-timescale4 “chopping” TTV associated with the planetary synodic timescale has been identified as an important feature for unique characterization of systems (Holman et al., 2010; Nesvorný et al., 2013). More recently, Nesvorný & Vokrouhlický (2014) studied this chopping TTV to clarify how, despite degeneracies between parameters, the TTV method can be used to measure planetary masses in the case of low-eccentricity orbits.

In Section 2, we begin by describing a harmonic approach to analyzing TTVs. We review the work of Lithwick et al. (2012) in Section 3, and discuss the TTV signal for a system near a first order mean motion resonance (referred to hereafter as the “Lithwick et al. formula”). We then introduce the conjunction effect and give analytic formulae for the chopping signal in Section 4; these were derived first in Agol et al. (2005) and more recently, using a similar approach, in Nesvorný & Vokrouhlický (2014). We show how the chopping formula encompasses near-resonant effects, and more generally consider the range of validity of the synodic TTV formulae. We then place the synodic TTV expression in context with that of Lithwick et al. and discuss the regimes in which each should be used. In Section 4.4, we address the more general problem of using these formulae to predict, given the timing precision on the transits, how many observations are required to infer the masses of a particular system. This will be important in the planning of follow-up observations of partially characterized systems and for estimating the timing precision required for future surveys to obtain a measurement of planetary masses through chopping.

In Section 5, we apply the synodic TTV formula to Kepler data, and use it to infer planetary masses for systems both near and far from mean motion resonance. The synodic formulae can be used alone to determine planetary masses, or in combination with the Lithwick et al. formula for systems near first order mean motion resonances, in which case the mass-free eccentricity degeneracy can be broken and a constraint on the free eccentricities can be determined as well. We present our conclusions in Section 6. In the Appendix we give an alternate derivation of the synodic TTV formulae based on Hamiltonian perturbation theory, and we discuss the convergence of the series for the synodic chopping signal.

## 2. Transit aliasing and harmonic analysis of TTVs

The TTVs of a planet can be written as a combination of periodic components with frequencies that are integer combinations of the two interacting planets’ orbital frequencies, , where are integers, and are the mean motions of the two planets (Nesvorný & Morbidelli, 2008; Nesvorný, 2009; Nesvorný & Beaugé, 2010). More explicitly, transit-timing variations for a two planet system can be expanded as:

 δt1,k = P1m2M⋆∑p,q[a1,p,qcos[(pn1−qn2)t1,k]+b1,p,qsin[(pn1−qn2)t1,k]], (2) δt2,k = P2m1M⋆∑p,q[a2,p,qcos[(pn1−qn2)t2,k]+b2,p,qsin[(pn1−qn2)t2,k]], (3)

where denotes the transit number, () are coefficients which are functions of the orbital elements of the planets (except the mean longitudes ), and therefore vary on timescales long compared to the orbital period; are the masses of the two planets; is the mass of the star; and is the th transit time of the th planet. We assume that the observation baseline is short compared to the secular timescales (which are typically “long” since they are proportional to ) so that treating the coefficients as constant is justified. We will discuss this further in Section 4.2. Note that the transit timing variations scale in proportion to the orbital period of each transiting planet and in proportion to the mass ratio of the perturbing planet. (These equations and scaling relations do not apply in mean-motion resonance (Agol et al., 2005).)

The transit times, are converted to transiting timing variations after removing a mean ephemeris, , where is an integer (Agol et al., 2005). Since transit timing variations are typically much smaller than the planetary orbital periods, , the planets’ transits are (nearly) sampled on their orbital frequencies, so

 δt1,k ≈ P1m2M⋆∑p,q[a1,p,qcos[(pn1−qn2)(t1,0+P1k)]+b1,p,qsin[(pn1−qn2)(t1,0+P1k)]], (4) δt2,k ≈ P2m1M⋆∑p,q[a2,p,qcos[(pn1−qn2)(t2,0+P2k)]+b2,p,qsin[(pn1−qn2)(t2,0+P2k)]], (5)

where we have dropped terms of order . However, since , these equations can be rewritten as:

 δt1,k ≈ P1m2M⋆∑q[a′1,qcos2πq(P1/P2)k+b′1,qsin2πq(P1/P2)k)], (6) δt2,k ≈ P2m1M⋆∑p[a′2,pcos2πp(P2/P1)k+b′2,psin2πp(P2/P1)k)], (7)

where, e.g. for the inner planet, the sum over is now implicit: for a specific value of , the coefficient for each possible value of has been absorbed into the new coefficients (for the outer planet, the sum over is now implicit). Therefore, each coefficient now contains variations due to multiple linearly independent frequencies of the form . But because of the sampling of the TTV on the transiting planet’s orbital period, frequencies that differ by integer multiples of the orbital frequency of the transiting planet are indistinguishable in transit timing variations. Therefore short-period TTVs (arising from short-timescale periodic terms in the variation of the orbital elements) can contribute to the same harmonic of the perturbing planet as resonant TTVs, those arising from any near-resonant configuration. For example, in a system near the 2:1 resonance, the fast frequency will, due to the discrete sampling at every transit of the inner planet, be aliased such that it appears as a sinusoid with the super-period, i.e. the same timescale as the resonant frequency (akin to the stroboscopic effect in which a spinning car wheel appears to be rotating at a slower rate due to the sampling rate).

The inner planet will have apparent TTV frequencies that are integer multiples of the perturbing planet’s orbital frequency, , while the outer planet TTVs will show frequencies that are integer multiples of the inner planet’s frequency, . If the perturbing planet does not transit the host star, this aliasing can mask the true period of the perturbing planet since we can add integer multiples of the transiting planet’s orbital frequency without affecting the TTV. If both planets transit the star, though, then both periods and orbital phases are known. Consequently, TTVs can be fit using a harmonic analysis of the perturbing planet’s orbital frequency if that planet has been identified. We utilize this harmonic analysis below to identify components of the TTV that are caused by the conjunctions of the planets (Nesvorný & Vokrouhlický, 2014), and thus depends on the difference between their mean longitudes, the synodic angle . But first we turn to a review of TTVs due to first-order resonant terms.

## 3. First-Order Resonant TTV signal

Here we summarize the formula for the largest amplitude component of transit timing variations of a pair of planets near a first order mean motion resonance, as derived by Lithwick et al. (2012), and remind the reader of the origin of the mass-eccentricity degeneracy.

Lithwick et al. (2012) considered a system of two planets near the : first order resonance, in which case the TTVs take the following approximate form:

 δt1 = R[−iV1expiλjR] (8) δt2 = R[−iV2expiλjR]; (9)

here denotes and

 V1 = P1π1j2/3R(jR−1)1/3Δm2M⋆[−f−3Z∗free2Δ] (10) V2 = P2π1jRΔm1M⋆[−g+3Z∗free2Δ] (11) λjR = jRλ2−(jR−1)λ1, (12)

with

 Zfree = feiϖ11+geiϖ22 (13) Δ = P2P1jR−1jR−1 (14)

and and are combinations of Laplace coefficients (note that and ), the free eccentricities, and the longitudes of pericenter of each of the planets. The reference direction from which longitudes are measured is the observer’s line of sight, so that at transit. Note that , , and are complex quantities; however, all observables are found by taking the real components as in Equation 8. ( is the complex conjugate of , and and denote the real and imaginary components of a complex number ). In the derivation of these formulae, only the two resonant terms associated with the : resonance at first order in eccentricity are considered. All other terms in the gravitational potential between the planets, even those of independent of eccentricity, are neglected since near resonance the TTVs they produce are very small compared to those caused by the resonant terms (due to the small denominator appearing in the amplitudes given in Equations 10).

The total eccentricity of a planet near a first-order mean motion resonance is made up of “free” and “forced” components; the forced eccentricity is driven by resonant interactions between the planets while the free eccentricity is determined by matching the initial eccentricity vector (the initial conditions). The free eccentricity vector can be approximated as constant on observational timescales (it varies on the long secular timescale) while the forced eccentricity vector precesses on the resonant timescale (the super-period). Both components affect the TTVs: the forced eccentricity is responsible for the first term in the brackets in Equations 10 (depending on the coefficients and ), while the free eccentricities are responsible for the term. The amplitude of the free eccentricity is not necessarily smaller or larger than that of the forced eccentricity.

In the case of two transiting planets, the unknown quantities are , , , and , and, in principle, there are four measurable quantities: both amplitudes and both phases. However, as pointed out in Lithwick et al. (2012), in both the case of and the TTVs are approximately anti-correlated, and so in either of these regimes the relative phase is no longer a constraining quantity (if the signal-to-noise is insufficient to measure a phase offset). This leads to the degeneracy between mass and free eccentricities, and hence from this first-order TTV alone one in practice cannot usually determine the masses or the combination of free eccentricities .

Finally, since the TTVs are sampled for the inner [outer] planet at [], these expressions become:

 δt1 = R[−iV1exp(ijRλ2)] (15) δt2 = R[−iV2exp(i(jR−1)λ1)], (16)

Actually, Lithwick et al. have taken advantage of the aliasing effect to create these expressions. The component of the TTVs due to the forced eccentricity has a dependence for the inner planet and a dependence of the outer planet (see their Equations (A15) and (A24)). Because these appear at the same aliased frequency as the resonant term, they are included with the resonant TTV. Since the reference direction was chosen such that at transit, these two terms have the same phase and can be grouped together in this way.

## 4. The Synodic Chopping Signal

A pair of planets interacts most strongly when the distance between them is smallest. For low eccentricity and nearly coplanar orbits, this occurs at conjunction, when the synodic angle (). Conjunctions occur periodically, with a timescale of , and it is intuitive to expect that this timescale would appear in the TTVs. In fact, at zeroth order in the eccentricities, the transit timing variations only depend on the synodic angle and its harmonics,

 ψj=j(λ1−λ2); (17)

see Appendix for more detail.

### 4.1. Synodic chopping signal formulae to zeroth order in the free eccentricities

The synodic chopping signal is included in the computations in Agol et al. (2005) and Nesvorný & Vokrouhlický (2014). We also give a distinct derivation in the Appendix here which casts the expressions in a form useful for the harmonic analysis described in Section 2. All three formulae agree in the limit that the reflex motion of the star can be ignored (after correcting a typo in Nesvorný & Vokrouhlický (2014); see Appendix for more detail).

For an inner transiting planet, the synodic component of the transit timing variations takes the form:

 δt1=∞∑j=1P12πm2M∗f(j)1(α)sinψj, (18)

where

 f(j)1(α)=−αj(β2+3)bj1/2(α)+2βDαbj1/2(α)−αδj,1(β2+2β+3)β2(β2−1), (19)

where , is the Laplace coefficient

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

and is the derivative operator .

The Laplace coefficients can be evaluated in terms of complete elliptic integrals; for example:

 b(1)1/2(α) = 4(K(α)−E(α))/(απ), (21) b(1)1/2(α)+α∂b(1)1/2(α)∂α = 4αE(α)π(1−α2), (22)

where and are the complete elliptic integrals of the first and second kinds, respectively.

For an outer transiting planet, the synodic component of the transit timing variations takes the form:

 δt2=∞∑j=1P22πm1M∗f(j)2(α)sinψj, (23)

where

 f(j)2(α)=j(κ2+3)bj1/2(α)+2κ(Dαbj1/2(α)+bj1/2(α))−α−2δj,1(κ2−2κ+3)κ2(κ2−1), (24)

where .

In the limit of large , the function has a leading coefficient which scales like and the function has a leading coefficient which scales like (see Appendix), so that the largest contributions to the sum in Equation 18 and in Equation 23 in general come from smaller values of . Because of the behavior of the leading coefficients, more terms are necessary to faithfully approximate the TTV signal as because the convergence of the sums is slow.

Closely spaced planets will, at a given , have a larger synodic TTV than more widely spaced planets because the parameters and appearing in the denominators approach zero as . Additionally, if the pair is near the mean motion resonance, the denominators and will be close to zero for the term with (for the inner planet) and (for the outer planet). Near the resonant configuration, then, these terms dominate the synodic TTV because of these small denominators. This reflects the fact that near the mean motion resonance, the time between conjunctions approximately corresponds to orbits of the outer planet and orbits of the inner one, such that most of the TTV amplitude is incorporated in these harmonics of the synodic angle. Additionally, these particular harmonics of the synodic angle contribute to the TTV as a long period effect. Due to aliasing, they have a timescale given by the super-period of the first order resonance.

The magnitude of the functions and are plotted in Figure 1 for . It is clear that the synodic amplitude generally grows as , for all , and therefore more terms must be included in the sum for close pairs of planets. Furthermore, when the period ratio is close to the mean motion resonance, the functions with (inner planet) and (outer planet) peak for the reasons discussed above. This happens for all at some period ratio except for the synodic TTV of the inner planet since . As noted in Agol et al. (2005), the dip near a period ratio of 2.5 in the synodic signal of the outer planet is due to the fact that the TTV caused by the motion of the star about the barycenter of the inner planet-star binary subsystem is opposite in sign and comparable to that caused by direct interaction with the inner planet at this period ratio. For larger period ratios, the TTV of the outer planet is dominated by the component due to the motion of the star about the barycenter of the inner planet-star binary subsystem.

Note that these relations have no eccentricity dependence; they only depend on the semimajor axis ratio of the two planets, , and the planet-star mass ratios. The phase and period of the chopping signal are straightforward to measure if the perturbing planet also transits the star; thus, if this synodic chopping signal can be measured, the mass ratio of the perturbing planet can be immediately inferred if the parameters of the system satisfy the major assumptions made in this paper (see Section 4.2).

### 4.2. Range of validity of the synodic formula

The intrinsic assumptions made in deriving these TTV expressions were 1) that the system has low eccentricities and nearly coplanar orbits; 2) that the system is not in or too near resonance; 3) that is not too close to unity; 4) that the masses of the planets are small compared to that of the star; and 5) that the coefficients themselves can be treated as constant in time. Additionally, if two perturbing planets contribute to the TTVs of a third planet, the true TTV of this third planet cannot be written as only a sum of one of the perturbing planet’s harmonics. In most configurations, we hypothesize that the TTV could be approximated as a simple sum of two “single-perturber” contributions.

We will focus on testing the validity of the formula in light of the intrinsic sources of error in the formula. First, we consider the qualification of “small” eccentricities and inclinations. A simple first guess of the order of magnitude error of these neglected terms would be that the coefficient would change from a value to , in which case if the eccentricities and inclinations are only a few percent than the error will also be a few percent. However, whether or not the eccentricities and inclinations can be considered “small” depends on how close the system is to resonance.

In general, if one derived a formula for the TTV good to first order in eccentricity, it would have terms at first order in eccentricity only depending on the angles of (not including the longitude of pericenter piece). In fact, the TTV formulae themselves should have the d’Alembert characteristics since they only depend on angles referenced to a fixed direction (longitudes) (Hamilton, 1994). Due to the transit aliasing effect discussed above, these first-order frequencies will be indistinguishable from variations of frequency for the inner planet since these appear as the same harmonic of the outer planet, (see Equation 15). For the outer planet, these first order (in eccentricity) terms will indistinguishable from those at the harmonic of .

These first-order terms will have amplitude proportional to , which, away from resonance, should be negligible for all if the eccentricity is low. But near the : resonance, the (for the inner planet) and the term (for the outer planet) will have amplitudes proportional to (appearing as in the formulae of Lithwick et al, Equation 10). The resonant combination of frequencies in the denominator mitigates the effects of the eccentricity coefficient, so that these terms are large corrections to the TTV formulae at these frequencies, and hence they make the synodic formula for these values of less accurate. Note that since , the synodic frequency, , has no resonant aliases for the inner planet. Thus the harmonic component which depends on can be uniquely identified with the synodic frequency variation.

If the system is in a mean motion resonance, the dominant, resonant TTV period will be related to the libration time. In this case, the resonant contribution to the TTVs will not appear at the frequency for the inner planet and for the outer planet, since the libration time cannot be written simply as an integer times the orbital frequency. In the resonant case, then, the harmonic analysis approach cannot be used.

It is important to point out that in each of these cases - both near the : resonance, where eccentricity errors are larger, and in the mean motion resonance, where the derivation in the Appendix must be modified, the synodic TTV formula for values of (inner planet) and (outer planet) will still apply.

In Figure 2, we show how the predicted amplitude of the synodic chopping term compares with that calculated numerically, for both the inner and outer planet, assuming the pair is near the 2:1 resonance. We varied only the period of the outer planet and the eccentricity of both planets (assumed to be equal; the vertical scale is the square root of the sum of the squares of the eccentricities, ). The mass of each planet was set to be . For each planet, we first determined the TTVs by fitting a line to the transit times, and we then used the computed average orbital periods to determine the value of used in the synodic formulae with for each planet. We then determined the numerical amplitude of the chopping signal by fitting 10 harmonics of the perturber’s period along with a linear term to 1,000 simulated TTVs of the inner planet and between 460 and 540 simulated TTVs of the outer planet and selecting the harmonic. This experiment therefore represents an ideal case where the main source of error comes from the assumptions made in deriving the formulae, and not from issues with not having enough data and/or precision to resolve the amplitudes of the various harmonics.

The error in the formula for the inner planet is below across the entire range studied, except very near the 2:1 resonance itself. This is as expected - near the 2:1 resonance, only the contribution of the TTV of the inner planet is expected to have large corrections due to eccentricity effects. The configurations with errors of are likely those in the mean motion resonance, where the harmonic approach does not apply (the libration frequency is not simply aliased with the frequency ), or those where the super period is significantly longer than the simulation time (to be discussed at the end of this section). For the outer planet, the synodic term is aliased with the resonant term, and hence we expect large errors as a function of eccentricity. Indeed, even relatively far from resonance the synodic formula with fails (errors larger than 10%) for eccentricities larger than 0.04.

We performed the same numerical experiment varying the number of harmonics fit to the TTVs, and we also decreased the simulation time. We found the same results even with fewer harmonics (in all cases, the number of data points was much larger than the number of free parameters). With shorter simulation times, we also recovered the same results except when the simulation time became significantly shorter than the TTV period (in which case the error is not due to the formula, but due to insufficient coverage).

Although we have not calculated the effects of inclination on the TTVs, if the TTVs in fact follow the d’Alembert characteristics then the correction due to inclination will only appear at second order. Therefore, one would expect that the inclination terms neglected will be a smaller source of error than the eccentricity terms neglected. Indeed, Nesvorný & Vokrouhlický (2014) did not find errors larger than 20% in the synodic chopping formula (away from resonance) until the mutual inclination was larger than .

In deriving the TTV formula, we began from the disturbing function, which assumes that the interaction between the two planets can be written as a converging series in . This means that the formulae will not work for co-orbital planets (see Vokrouhlický & Nesvorný (2014) for an analysis of that case). Additionally, as , the Laplace coefficents converge less quickly, and so higher order eccentricity and inclination terms, ignored in the derivation, are potentially more important. A different issue is that as , mean motion resonances are densely spaced (for arbitrary , the system is more likely to be close to a resonance if is closer to unity). In principle then the (neglected) effects of eccentricity could be more important for closer pairs of planets due to the effect of the small denominators discussed above. Note, however, again the special case of the synodic TTV of the inner planet. At there are no possible small denominators in the TTV of the inner planet, since all resonances (except the 1:1) require .

In Figure 3, we show the error in the formula for the synodic TTV, for the inner planet, across a wide range of orbital separation and eccentricities. In making this plot, 500 transits of the inner planet were simulated and 10 harmonics were fit to the resulting TTVs. Except at resonance, the error is less than , regardless of the eccentricity. This indicates that the neglected eccentricity terms are in fact small across a wide range of . Note that the error is unbiased, in that it is not typically positive or negative.

In deriving the synodic formulae, we assumed that the masses of the planets, compared to that of the host star, are small, so that neglected corrections of order are small. In Figure 4, we show how the fractional error in the synodic TTV formula with for the inner planet grows as we increase the masses of the planets to and then to . Some of these configurations tested were unstable, and some were perturbed enough to change the number of transits by more than 1, and those were not fit (though shown in magenta). Note that the widths of resonances grow as the mass of the planets grow, and so more systems are in resonance, where the harmonic analysis will not work, than in the fiducial case shown in Figure 3. In the case of Jupiter mass planets, there are configurations where the chopping formula may be too approximate, as even outside of resonance, between the 3:2 and the 2:1, the error is on the order of . As Nesvorný & Vokrouhlický (2014) point out, in these cases the chopping formula primarily provides motivation and understanding as to how mass measurements from TTVs arise.

Lastly, we assume that the coefficients of the synodic TTV, which depend on , can be treated as fixed in time at their average values. It is possible, especially in the near resonant case where the TTV timescale is long, that the observations will cover only a small fraction of the TTV cycle. In this case, the observed semimajor axis ratio may be different than the average one. The fractional error resulting from using the “incorrect” value of in the formula for will be of order . Although is small, of the same order of the TTV compared to the orbital period, the derivative of the coefficient can be large near resonance. Therefore, if one observes only a small fraction of the super-period of a near resonant system, there will be errors relating to an incorrect estimate of the average value of . These errors are larger for the terms in the synodic sum aliased with the resonant frequency.

The long-period oscillations in the eccentricity and inclinations due to secular effects will not be resolved, since the secular timescale is in general very long compared to observational timescales (for example, on secular timescales will change). The synodic chopping formula, which is independent of eccentricity and inclination, will therefore have slowly varying error terms, but as long as the eccentricities and inclination remain small over the secular timescale the error terms will remain small as well.

### 4.3. Comparison to the Lithwick et al. formula with Zfree=0

How do the synodic chopping formulae compare with the Lithwick et al. formulae? First, the Lithwick et al. formulae applies for pairs of planets near a first order mean motion resonance, and it includes the contribution to the TTV at the frequencies for the inner planet and for the outer planet. However, since only the resonant terms were used in deriving the Lithwick et al. formulae, some synodic effects which get aliased to the resonant frequencies were neglected (these neglected terms are small amplitude near resonance since they do not have the small denominator of ). On the other hand, the synodic formula applies for pairs both far and near mean motion resonances, and encompasses the effects of conjunctions at every harmonic in the TTVs of the inner planet and in the TTVs of the outer planet. However, the Lithwick et al. formulae includes the approximate first order eccentricity correction for pairs near a mean motion resonance, while the synodic formulae only hold at zeroth order in eccentricity.

We now compare the two sets of formulae in the regime where they should agree: near resonance (where the synodic chopping terms without the small denominator ignored by Lithwick et al. (2012) are negligible), with , and for the correct value of chosen in the synodic sum (that of for the inner planet and for the outer planet). Our expectation is borne out by a numerical comparison between the two, the results of which are shown in Figure 5. In short, the Lithwick et al. expression is an excellent approximation to the TTV of systems near first order mean motion resonances, while further away from resonance it becomes a worse approximation to the chopping signal with for the inner planet and for the outer planet. Correspondingly, the synodic formulae will be a worse approximation at these specific values of for systems closer to a first order resonance because of eccentricity effects, as we neglect the correction - see Section 4.2.

### 4.4. Measurement precision

One can use the synodic TTV formulae to estimate the precision on the mass measurement of a perturbing planet. For example, the expected mass precision (of the outer planet) due to synodic chopping with in the inner planet is given by

 σM2=σt1M∗2π√(Ntrans−Nparam)/2P1f11(α), (25)

where is the timing precision of the inner planet, is the number of model parameters, and is the total number of transits observed. This formula assumes that the phase of the sine function is adequately sampled so that the RMS of can be assumed, that there is no uncertainty on the mass of the star, , or , and that there is no covariance with other harmonics being fit. This formula also assumes that transits are observed continuously over the full super-period.

## 5. Applications

In this section we apply the chopping and resonant formulae to several planetary systems which have been analyzed using a full dynamical model involving numerical integration of the gravitational equations of motion. We instead analyze these systems using the synodic formulae for TTVs, in combination with the component of the Lithwick et al. formula for near resonant systems. As a result, we can measure planetary masses for these systems without full numerical analysis and without the complication of the mass-free eccentricity degeneracy inherent in the Lithwick et al. formula alone. This allows us to demonstrate empirically the validity of the synodic TTV signal and to strengthen the understanding of what information in a TTV signal leads to mass measurements.

For a system far from a first order resonance, we would recommend using the synodic formulae for the terms (for which ), and using the second component of the Lithwick et al. formula for the term for the resonant for the inner planet and for the outer planet. If the system is close to resonance, the synodic chopping formulae can be used for values of not aliased with the resonant frequencies. Note again that the resonant term and the synodic terms aliased with the resonant term will in general have different phases, unless the reference direction is chosen so that at transit. Additionally, while only the sinusoidal component of the harmonics should have nonzero amplitude, the resonant component has in general both sine and cosine components due to the complex quantity .

### 5.1. PH3/Kepler-289

The Kepler-289 (PH3/KOI-1353/KIC 7303287) planetary system was identified by the Kepler pipeline (Borucki et al., 2011; Batalha et al., 2013; Tenenbaum et al., 2013) and by the Planet Hunters crowd-sourced project (Fischer et al., 2012). This system consists of three planets with orbital periods near 35 days (Kepler-289b/PH3b), 66 days (PH3c), and 126 days (PH3d); each adjacent pair of planets is close to a period ratio of 1:1.9 (Schmitt et al., 2014). The outer two planets both display large amplitude transit-timing variations with a timescale of the super-period of the nearby 2:1 resonance, and the middle planet shows a strong chopping signal caused by the outermost planet. The masses of the outer two planets were measured by the transit timing variations through a full numerical analysis of the (assumed coplanar) system, performed by EA, as part of Schmitt et al. (2014). The inner planet does not have significantly detected TTVs, and does not significantly affect the outer two planets’ transit times, resulting in an upper limit on its mass only.

In this work, we returned to the published transit times and uncertainties of Schmitt et al. (2014) and analyzed them using the harmonic-fitting approach described in section 2 in order to measure the masses of PH3c and PH3d. Figure 6 shows our initial harmonic fit to the data; two harmonics are required for PH3c, while only one is required for 3d. We ignore the innermost planet (PH3b) in the analysis. We used the component of the Lithwick et al. formula with for modeling the component of the outer planet’s TTVs and for the component of the middle planet’s TTVs. We included the synodic chopping signal with for the middle planet PH3c, and the component for the outer planet (which encompasses the contribution from the Lithwick et al. formula).

As these three terms together only constrain a linear combination of the free eccentricities of the planets (), we also enforced as a prior the Hill stability criterion (Gladman, 1993) to prevent the eccentricities from growing too large. We added a systematic error parameter, , in quadrature to the measured timing errors, such that larger values of are penalized in the likelihood function while smaller values of require a closer fit to the transit times in order to have a high likelihood. We carried out an affine-invariant Markov chain analysis (Foreman-Mackey et al., 2013) with eleven free parameters: the ephemerides (, ), eccentricity vectors (, ), and masses of each planet, plus .

Figure 7 shows the confidence limits in the planet masses from the harmonic analysis with and without the constraint from the synodic chopping signal, as well as the confidence limits from the full dynamical analysis of Schmitt et al. (2014). As explained in Lithwick et al. (2012), the 2:1 resonant signal constrains a combination of the mass ratios of the planets and the free eccentricity, . Without the synodic chopping signal, our analysis shows a banana-like degeneracy between the two planet masses which is due to the trade-off between their masses and the free eccentricity (Figure 7, light blue), giving and . When the chopping signal in the TTVs of the middle planet is included, the mass ratio of the outer planet becomes constrained; this then breaks the mass/free eccentricity degeneracy, and allows the mass of the inner planet to be determined as well. The derived error ellipse is similar to that from the full dynamical analysis: Schmitt et al. (2014) report masses of and , while the harmonic analysis yields and .

This analysis demonstrates the power of the chopping signal in constraining planetary masses near a first-order mean motion resonance. In principle, although the chopping components (for the outer planet) and chopping components also provide independent constraints on the planetary masses, they are smaller in amplitude and not detected in this case.

### 5.2. Kepler 11d/e

The Kepler-11 system (Lissauer et al., 2011) is a system with six transiting planets. A full dynamical analysis has been carried out for this system, giving constraints on the masses of all planets (Migaszewski et al., 2013; Lissauer et al., 2013). Several of the planets, despite being only a few Earth masses, have low densities which require H/He atmospheres; this result is puzzling in light of core-accretion theory, which would not predict planets so low in mass to accumulate substantial gaseous envelopes. Here we validate the existing mass measurements of these two planets and we show that the mass constraints of Kepler-11d and Kepler 11-e largely result from the chopping TTV signal.

Kepler 11d and 11 e are the two of the three most massive planets in the Kepler-11 system, with periods near 23 and 32 days, respectively, in close proximity to 3:2 commensurability. Each of these planets have transit timing variations that are dominated by the other; thus they can be dynamically ‘decoupled’ from the rest of the planets, and treated as a two-planet system. Note, however, that the decoupling is “one-way”: Kepler 11e affects the TTVs of planet Kepler 11-f, and hence there is more information with regards to the masses of d and e to be gained by fitting the entire system instead of treating the (d,e) pair in isolation.

Figure 8 shows the harmonic fitting results for Kepler 11d/e using the transit times due to Jason Rowe presented in Lissauer et al. (2013), to be compared to the dynamical constraints in Table 7 of that paper. In black circles, we show the actual TTV measurements for each planet, with corresponding uncertainties. In this case we simply fit for the harmonics of the TTV with the frequency of the companion planet up to (dotted lines), which resulted in excellent fits for both; the synodic chopping signal from the harmonic fit is also plotted (in red). Note that near the 3:2 mean motion resonance, the synodic signal is not aliased with the resonant frequencies, and so we expect the chopping signal to be well approximated by our formula, as discussed in 4.2. We over-plot the predicted synodic chopping signal based on the the best-fit mass ratios from Table 7 in Lissauer et al. (2013), shown as the blue curves. This shows that the chopping signal is detected for both planets, and that it is consistent with the chopping signal predicted by the full dynamical analysis.

Next, we carried out a Markov chain analysis for these planets including the Lithwick et al. resonance formulae with , relevant to the 3:2 commensurability (Section 3), as well as both the inner and outer chopping formulae, summed to . Figure 9 shows the constraints on the masses of the two planets. The black curve shows the 1 confidence limit from dynamical analysis in Lissauer et al. (2013). In dark[light] blue is the 1[2] confidence limit for our analysis with the resonant and full chopping signals included for both planets; this is consistent with the dynamical analysis at the level, albeit with a larger uncertainty (recall again that the TTVs of planet Kepler 11-f are affected by Kepler 11-e, and so there is more information as to the (d,e) subsystem, available when fitting the whole system).

Another interesting byproduct of the chopping signal, when using in conjunction with the Lithwick et al. formula, is that it allows for a measurement of the quantity . In the case of Kepler 11 d/e, we find a value of . Since dissipation of eccentricities first damps the free eccentricities, it is interesting that the here, though modest, is distinctly nonzero.

In general, because the chopping amplitude function is smaller in magnitude for larger values of , the components with larger may not be measurable. However, if they are, they can be used to provide additional constraints on the masses as consistency checks. For example, in the case of Kepler 11 d, the chopping signals with all independently constrain the mass of planet e, while the chopping signals of present in the TTVs of e all independently constrain the mass of planet d. For planet d, the inferred mass (assuming in this case a one solar mass star) is or , compared to in Lissauer et al. (2013). For planet e the inferred mass is or , compared to in Lissauer et al. (2013). Given that these estimates agree at the level, this indicates that the dynamical analysis yields masses that are consistent with the chopping amplitudes.

### 5.3. Kepler 9

In some cases, the Lithwick et al. formula alone can be used to determine the masses of the planets. In practice this requires that the amplitudes and phases of the TTVs must be measured with high accuracy.

The first system with detected transit timing variations is the Kepler-9 system (Holman et al., 2010), which consists of three planets, the outer two of which are close to 2:1 period commensurability with periods near 19 and 39 days, respectively. The outer pair is dynamically decoupled from the inner planet in that the inner planet does not measurably affect their TTVs. A recent dynamical analysis of Kepler-9 shows that the transit timing variations yield masses of the outer two planets of and (Dreizler & Ofir, 2014b). We carried out a harmonic fit to the transit times for each planet, and find an excellent fit to planet Kepler-9c[d]’s transit times with four[six] harmonics of the period of Kepler-9d[c]. Figure 10 shows the transit timing variations for both planets along with the harmonic fit. The synodic chopping () amplitude of the inner planet matches the phase and amplitude predicted based on the mass inferred by Dreizler & Ofir (2014b).

We carried out a Markov chain analysis on the set of transit times published by Dreizler & Ofir (2014b) for this system with the 2:1 resonant term and chopping, as well as additional harmonics with amplitudes that were not constrained by the physical parameters of the model: for the inner planet we added harmonics at , and , while for the outer planet we added harmonics at , and to our model. We find similar masses as Dreizler & Ofir (2014b) and , albeit with larger uncertainties. When we remove the constraint on the amplitude of the chopping signal, we find comparable masses and uncertainties; we suspect that the reason for this is that in this case the amplitudes of the TTV are measured with sufficient precision that the and terms that occur in the first order resonant TTV formula can be distinguished in amplitude. The imaginary component of the TTV is significant for both planets, so in this case one can break the degeneracy between the planet masses with just the resonant term. The amplitudes for the resonant term for both planets have four (well-measured) constraints, the real and imaginary amplitude for each planet, while there are four unknowns: the mass ratios for each planet, and the real and imaginary component of . This gives a unique solution, so the masses are well determined without the need for the chopping constraint.

Although chopping is not required to determine the mass of the planets, we can show that in this case it is consistent with the masses inferred from the resonant terms alone. The measured mass of the outer planet predicts the amplitude of the synodic chopping signal of the inner planet, which only has dependence. In Figure 11 we show that the measured sine amplitude of the inner chopping signal is consistent to with the predicted amplitude. The amplitude of the synodic chopping term gives a mass of the outer planet of , while the mass estimated from the resonant term is . In addition, the amplitude of the cosine term is consistent with zero at , as it should be for the synodic chopping term.

When radial velocities are included in the analysis, a larger mass is derived for the planets, yielding about 55 for the outer planet (Dreizler & Ofir, 2014b). This is inconsistent with the chopping signal, at about the 3 level (see green point in Figure 11) and inconsistent with the masses derived from resonant TTV alone. This discrepancy indicates that there is still some tension between the TTV data and the RV data, possibly due to RV jitter, additional planets (causing perturbations of the RV velocities), systematic errors in the transit times (perhaps due to star spot crossings), or, perhaps, simply statistical fluctuations. The fact that the resonant and chopping terms give similar estimates of the outer planet mass increases our confidence that the transit timing analysis is not strongly affected by additional planets in this system.

For this system we also tried to use the amplitude of the chopping signal for the outer planet to constrain the mass of the inner planet. However, the amplitude is much too large to be due to chopping, and instead we believe is due to the 2:4 resonant term (which is of order , but these planets are so close to 1:2, that the term compensates for this).

### 5.4. Mass precision of KOI-872c

For planets that are not near a first-order mean-motion resonance (), the synodic chopping amplitude can provide the strongest constraint upon the planet masses. If there are a large number of transits, then the signal-to-noise of the planet mass ratio can be estimated with equation (25).

The first non-transiting planet found with transit timing variations (with a unique identification of the perturbing planet’s period) occurred in the system KOI-872 (Nesvorný et al., 2012), in which the period ratio of the two planets is close to 5:3. We carried out a harmonic analysis of the transit times of the inner planet with harmonics up to , starting with the period of the perturbing planet, KOI-872c, from the published dynamical analysis. The synodic chopping signal, the coefficient of the term, was measured to have an amplitude of days for the sinusoidal component (detected at 32), and for the cosine component (consistent with zero at ). Applying the synodic chopping amplitude formula (Equation 18, with ), we find a mass ratio of: . This compares favorably with the mass ratio measured from the full dynamical analysis (Nesvorný et al., 2012) of: , with a similar magnitude uncertainty. Figure 12 shows the harmonic fit to the transit timing variations of KOI-872, with the synodic chopping signal shown with red points. The predicted synodic chopping based upon the dynamical solution is shown in blue, demonstrating agreement with the derived signal (albeit discrepant by ).

The expected precision in the planet mass is given by equation 25; for KOI-872 there are 37 transit times with a typical precision of 0.0015 days, 13 model parameters, giving an expected mass precision of , which matches well that found with the harmonic fits, and is close to the uncertainty found with the full dynamical analysis. This indicates that the timing precision along with the total number of transits observed and the number of free parameters fit can be used to forecast the mass measurement precision in the non-resonant case. A caveat is that this formula applies when chopping dominates the mass uncertainty, which may not be true for large eccentricities of the planets which can cause higher order resonant terms to play a more important role.

### 5.5. Predicting planetary mass precision inferred from the synodic TTV

The chopping effect potentially allows for mass measurements of planets with TTVs observed by a TESS-like mission because the synodic TTV, when unaliased with a resonant frequency, is a short-period effect. However, the amplitude of the effect is also considerably smaller than that due to the long-timescale resonant variations, and so more transits are needed to build up signal to noise. Here we consider a system with 2 planets with a period ratio of 1.5, with an inner orbital period of 20 days. We assume that the system has been observed for 1 year (the baseline TESS will have for stars near the celestial poles), and assume that the timing uncertainty on the transit times of the inner planet are 1 minute.

In this case, and the formula in Equation 25 yields, for 10 free parameters and assuming a solar mass star, a mass precision of the outer planet based on the chopping of the inner planet of . The mass precision on the inner planet due to the chopping in the outer planet is . Wider pairs will have larger mass uncertainty since the function is smaller. For example, this same pair moved to the 2:1 resonance will allow an mass uncertainty for the outer planet.

Note that the mass uncertainty scales like the inverse of the square of the orbital period, and, given an orbital period, like the inverse of the square of the observation time. Therefore, for a longer mission like PLATO, with an observational baseline of 2 or 3 years, the uncertainty on the mass of the perturbing planet will be smaller by a factor of or , respectively, compared with that estimated for 1 year of data above.

## 6. Conclusions

In this paper we have written down expressions for transit-timing variations in the plane-parallel limit, in the limit of zero free eccentricity of both planets, outside of resonance, and for timescales shorter than the secular timescale. Despite these assumptions, these terms have important consequences for analysis of transit-timing variations of multi-planet systems: 1) the TTVs have a dependence on , for to ; 2) the amplitude of these terms only depends on the mass-ratio of the perturbing planet to the star, and the semi-major axis ratio of these planets. Although other papers have presented formulae in the limit of zero-eccentricity (Agol et al., 2005; Nesvorný & Vokrouhlický, 2014), this is the first time that the coefficients for each have been written down explicitly. This allows for harmonic analysis of transit times in terms of the period and phase of the perturbing planet; the coefficients of each harmonic can then be related to the planet mass ratios and eccentricities via these formulae (except for the harmonics affected by resonant terms). When the period and phase of the perturbing planet are known, then fitting for the harmonic coefficients (and the mean ephemeris) is a linear regression problem; thus a global solution can be found by simple matrix inversion, yielding a unique solution for the coefficients. This means that there cannot be multi-modal degeneracies in the derived masses of the planets. The amplitudes and uncertainties of these coefficients can then be translated into planet masses using the formulae given above. Alternatively the coefficients can computed from the physical properties of the planets (masses, eccentricity vector, and ephemerides), and then the model can be fit to the data with a non-linear optimization or Markov chain, as we have carried out for several of the examples presented above.

In particular, we have highlighted the importance of the ‘synodic chopping’ signal of the inner planet, which has a weak dependence on the eccentricities of the planets, and is not aliased with any first-order resonant terms. The formula for this term in particular has very little error even for very compact orbits. For the outer planet, the chopping signal can be used if the period ratios are distant from the 2:1 resonance, though if the transit times are precisely measured the synodic TTV resulting from values of may be used to constrain the mass of the inner planet instead. In general, if the system is near or in a first order resonance, with , only the synodic TTV component with (inner planet) and (outer planet) will be altered; the formulae with not equal to these values will apply. Other limitations of this formula is that it cannot be directly used when more than one planet strongly perturbs the transiting planet, it breaks down for very massive planets () with period ratios less than 2, and that as the period ratio approaches unity the error due to neglected eccentricity terms can become more important in the formula (though this does not apply to the term in the TTVs of the inner planet). Following Nesvorný & Vokrouhlický (2014), we conclude that in these regimes the chopping formula provides insight into mass measurements with TTVs even though the analytic formulae may be too approximate to be applied directly.

We have applied these formulae to existing transiting planet systems that have been analyzed in prior publications, recovering the mass measurements, but using harmonic fits and the analytic chopping/resonant expressions rather than full N-body integrations. In some cases (KOI-872) this shows that the primary constraint on the mass of the planets comes from the synodic chopping signal. In other cases (KOI-1353c/d, Kepler 11d,e) the primary mass constraint comes from the combination of the first-order resonant signal and the synodic chopping of the inner planet. In the case of Kepler-9, the primary constraint comes from the resonant signal, although the chopping component gives a consistent constraint on the mass of the outer planet.

In future applications, we expect that these formulae can be used for rapid fitting and estimation of transit timing variations, for rapid estimation of planet masses, for initialization of N-body integration, for determining the requisite timing precision to measure planet masses with future follow-up observations of multi-planet transiting systems, and for forecast of transit timing variations. It should be possible to apply these formulae in systems of more than two planets using linear-combinations of the TTVs induced by more than one perturbing planet.

We would like to thank Dan Fabrycky, Eric Ford, Matt Holman, Daniel Jontof-Hunter, Jack Lissauer, Jason Steffen, and the Kepler TTV group for helpful conservations. E.A. acknowledges funding by NSF Career Grant AST 0645416, NASA Astrobiology Institute’s Virtual Planetary Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C, and NASA Origins of Solar Systems Grant 12-OSS12-0011. K.M.D. acknowledges supports from JCPA fellowship at Caltech. Work by K.M.D was supported by NASA under grant NNX09AB28G from the Kepler Participating Scientist Program and grants NNX09AB33G and NNX13A124G under the Origins program.

## Appendix A Derivation of synodic TTV

Here we provide an alternate derivation of the synodic TTV to zeroth order in planetary eccentricities. We follow the method first described in Nesvorný & Morbidelli (2008), and work at zeroth order in eccentricities and inclinations and first order in the parameter .

As mentioned in the main text, the synodic TTV has been derived before, by Agol et al. (2005) and Nesvorný & Vokrouhlický (2014). These, and the derivation below, give agreement in the limit that the reflex motion of the star can be ignored (after correcting a typo in Nesvorný & Vokrouhlický (2014)). For the outer planet, the reflex motion of the star dominates at large period ratios. This effect is accounted for in the equations in the appendix of Agol et al. (2005) and in the alternate derivation in the appendix here, but is missing from the terms in Nesvorný & Vokrouhlický (2014) (the Agol et al. (2005) calculation utilized heliocentric coordinates, while that of Nesvorný & Vokrouhlický (2014) used Jacobi coordinates and thus computed the TTVs of the outer planet relative to the center of mass, not relative to the star. Here we also employ Jacobi coordinates but correct for the reflex contribution afterwards.

TTVs are determined by taking an observed set of transit times and performing a linear fit to them, such that the transit timing variation can be written as

 δtn =tn−(t0+n¯P) (A1)

The slope is the average time between successive transits during the observational baseline, or the average period over this timespan. In this sense, transit timing variations are the deviations from the transit times predicted by the average of the true (non-Keplerian) orbit.

Transits occur when the true longitude of the transiting planet is equal to a particular value. The expression for is, to first order in eccentricity,

 θ =f+ϖ =M+2esinM+ϖ+O(e2) =λ+2esin(λ−ϖ)+O(e2) (A2)

where is the eccentricity, is the true anomaly, the longitude of periastron, the mean anomaly, and the mean longitude of the planet. Because we will perform the calculation using a Hamiltonian formalism, we now switch to canonical coordinates. Our variables will be

 Λi=mi√GM⋆ai λi xi=√2Picospi yi=√2Pisinpi

where:

 Pi=Λie2i2+O(e4i) pi=−ϖi

where for the two planets, is the mass of the planet, the mass of the star, and all orbital elements are Jacobi elements. Variables in the left columns are the canonical momenta, while those in the right columns are the conjugate coordinate. In terms of this canonical set, Equation A for becomes

 θ =λ+2√Λ(xsinλ+ycosλ) (A3)

We will perturb Equation A3 about the averaged orbit, and keep terms only at zeroth order in and :

 δθ =δλ+2√Λ(δxsinλ+δycosλ) (A4)

where as we will see the parameters are going to be of order and independent of eccentricities to lowest order. It is important to note that we define for any arbitrary function , where over-bar denotes a time-average. In that case,

 δθ =−nδt (A5)

where is the average mean motion and is the timing variation. Note that 1) we are consistent with the style that a negative timing variation corresponds to an “early” transit with and 2) we can use the average mean motion because Equation A4 holds only at first order in the perturbations; in the same sense, we can treat quantities without a on the right hand side as averaged variables as well.

The question now arises - how do we determine the perturbations to the average orbit and ?

The Hamiltonian for a system of two planets of mass orbiting a much more massive star of mass , written in Jacobi coordinates and momenta , takes the form, to first order in combinations of , of

 H =H0+H1 H0 =HKepler,1+HKepler,2 HKepler,i =p2i2~mi−G~Mi,⋆~mi|ri| H1 (A6)

where , and denote Jacobi masses. Note that the perturbation takes the functional form of the disturbing function with an exterior perturber (Murray & Dermott, 1999). We set and also ignore the difference between Jacobi and physical masses. In the Keplerian piece , this approximation corresponds to a (constant) change in the mean motions of the planets by order , but since we are interested in TTVs this constant change does not matter. The correction between Jacobi and physical masses in the perturbation generates only terms, and we ignore these.

When expressed in terms of the canonical set given in Equation A and Equation A, the Hamiltonian takes the form:

 H =H0(Λi)+ϵ1H1(Λi,Pi,λi,pi)+O(ϵ2) (A7)

where e.g. , , and

 H0 =−μ12Λ21−μ22Λ22 H1 =−μ2Λ22[j=∞∑j=−∞gj,0(α)cosj(λ1−λ2)+ j=∞∑j=−∞gj,27(α)√2P1Λ1cos(jλ2−(j−1)λ1+p1)+ j=∞∑j=−∞gj,31(α)√2P2Λ2cos(jλ2−(j−1)λ1+p2)] (A8)

and , , , and , and are functions of Laplace coefficients with the indirect terms included (Murray & Dermott, 1999)5. The relevant functions of Laplace coefficients can be written as

 gj,0(α) =12bj1/2(α)−δj,1α gj,27(α) =12(−2j−αddα)bj1/2(α)+δj,132α−δj,−112α gj,31(α) =12(−1+2j+αddα)bj−11/2(α)−δj,22α bj1/2(α) =1π∫2π0cos(jθ)√1−2αcosθ+α2dθ. (A9)

In writing the gravitational potential in this form, we have assumed that is not too close to unity, in which case Laplace coefficients converge slowly. For , the analogs of and appearing with higher powers of eccentricity are larger and so the neglected terms in Equation A become more important.

We now convert the variables (Equation A) to the set (Equation A), and also set and . Then the perturbation Hamiltonian takes the form:

 H1=−μ2Λ