Prospects for TTV Detection and Dynamical Constraints with TESS

The generation of the synthetic population of planetary systems follows the same basic procedure employed in Barclay et al. (2018). We provide a broad outline of the procedure here and detail some modifications to the original procedure of Barclay et al. (2018) that we make in order to more accurately model the architectures of multiplanet systems. We refer the reader to the original text for more details.

Barclay et al. (2018) simulated planets around stars in version 6 of the TESS Input Catalog (TIC) Candidate Target List (CTL). There are a total of 3.8 million stars in this catalog and all have computed stellar properties such as mass, radius, and brightness. Of these 3.8 million, 3.2 million were computed as observable to TESS, and were assumed to be observed by the mission using the 30 minute integration time FFI mode. From these, 214,000 were assumed to be observed using the spacecraft two-minute cadence mode. The stars observed at a two-minute cadence were selected using the CTL prioritization metric (Stassun et al. 2018), but we forced the stars to be evenly distributed over the observing sectors.

We modify the method of Barclay et al. (2018) for assigning planets to host stars in order to more accurately capture the architectures of multiplanet systems since planets' spacings strongly influence their TTVs. In particular, we aim to reproduce the period–ratio distribution of multiplanet systems inferred from Kepler using the debiased period–ratio distribution derived by Steffen & Hwang (2015). As in Barclay et al. (2018), we begin by assigning each star in the sample a number of planets drawn from a Poisson distribution. AFGK dwarfs were drawn assuming a mean of λ = 0.689 for periods ≤85 days (Fressin et al. 2013), while M dwarfs are drawn assuming λ = 2.5 for periods ≤200 days (Dressing & Charbonneau 2015). Next, we assign periods and radii to the innermost planet from the measured period–radius distribution of planets from Fressin et al. (2013) for AFGK dwarfs and from Dressing & Charbonneau (2015) for M dwarfs. We differ from Barclay et al. (2018) by drawing the periods of subsequent planets in the system from the period–ratio distribution in Steffen & Hwang (2015). Planets in a single system are drawn iteratively, the period ratio being with respect to the previous outermost planet at each step. Radii of these subsequent planets are drawn from the appropriate period bins of the inferred period–radius distributions of Fressin et al. (2013) or Dressing & Charbonneau (2015). Because giant planets almost never occur in multitransiting systems, planets in multitransiting systems initially assigned radii larger than 6 R_⊕ are reassigned new random radii from the period–radius distribution until a radius smaller than 6 R_⊕ is drawn.

The eccentricities of the planets were assumed to follow a Rayleigh distribution with a scale parameter of σ_e = 0.02, consistent with the eccentricity distribution of Kepler  multiplanet systems (Hadden & Lithwick 2014, 2017; Van Eylen & Albrecht 2015; Xie et al. 2016). The inclination of the innermost planet was drawn from a distribution uniform in $\cos i$, where i is the inclination with respect to the plane of the sky. All additional planets in any multiplanet system have sky-projected inclinations drawn from a Gaussian distribution with σ_i = 2° centered on the innermost planet's inclination. This inclination selection strategy mimics the mutual inclination distribution seen in Kepler's multis (Fang & Margot 2012; Fabrycky et al. 2014). Because constraints on the intrinsic mutual inclination distribution of Kepler multis are uncertain we also simulate a second, low-inclination population with σ_i = 15 to assess the influence of the assumed inclination distribution on our results. Results of our analysis for this low-i population are described in Appendix A. All results in the main body of the paper pertain to our fiducial σ_i = 2° population unless explicitly stated otherwise. Finally, each planet is assigned a mass randomly drawn using the probabilistic mass–radius relationship code, forcaster, of Chen & Kipping (2017).

We perform a basic check for stability of the multiplanet systems by determining whether each system's adjacent planet pairs satisfy the two-planet stability criterion of Hadden & Lithwick (2018). We remove 10 systems from our sample that host adjacent planet pairs predicted to be chaotic and unstable.

Having generated our sample of planetary systems, we use the detection model of Barclay et al. (2018) to determine which stars host at least one planet that is expected to be detectable by TESS. In brief, the detection model of Barclay et al. (2018) determines which planets will have transit detections with signal-to-noise ratio S/N ≥ 7.3 based on their transit properties, the predicted TESS photometric noise level (dependent on their host star's brightness; Stassun et al. 2018), and the flux contamination reported in the CTL. Applying the detection model, we find a total of 3756 stars host one or more planets with detectable transits. Table 1 summarizes the multiplicity statistics of all systems in which at least one transiting planet is detected. Of these 3756 stars, 1814 host more than one planet. These multiplanet systems are the focus of our TTV study. We note that the stellar hosts of the multiplanet systems are split such that 437 are M dwarfs, and 1377 of earlier type.

Of the 1814 multiplanet systems, 147 are multi-"tranet" systems (Tremaine & Dong 2012) containing two or more planets with detected transits while the other 1667 are single-"tranet" systems in which only one of the planets exhibit detectable transits while the other planets either do not transit or have transits below the S/N ≥ 7.3 threshold. These systems are detected in both two-minute cadence data and thirty-minute cadence FFI data: of the 1814 multiplanet systems, 775 are in the two-minute cadence in the nominal mission, while 1039 are in the FFI data.

In order to study what dynamical information may be gleaned from TTVs, we need an estimate of the precision with which TESS will be able to measure transit mid-times. Building on the work of Carter et al. (2008, 2009), Price & Rogers (2014) developed analytic expressions to approximate the uncertainties in various quantities derived from transit fitting as functions of planet properties and photometric precision. Their formula for the uncertainty in mid-transit times can be expressed as

$$\begin{eqnarray}{\sigma }_{{t}_{c}}=\displaystyle \frac{\sigma }{\delta }\sqrt{\displaystyle \frac{\tau }{2}}\times \left\{\begin{array}{ll}{\left(1-\tfrac{I}{3\tau }\right)}^{-1/2}, & \mathrm{if}\ \tau \geqslant I\\ {\left(1-\tfrac{\tau }{3I}\right)}^{-1/2}, & \mathrm{if}\ \tau \lt I\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

where σ is the uncertainty of individual TESS photometric measurements, δ is the planet's transit depth, I is the cadence (i.e., 2 minutes for pre-selected TESS targets and 30 minutes for planets discovered in the FFIs), and τ is the ingress time. We use Equation (1) to estimate the expected transit mid-time uncertainty, ${\sigma }_{{t}_{c}}$, for each planet in our sample, by combining the photometric uncertainty, σ, of the host star with the transit depth, duration, and ingress times computed from the planet's radius and orbital properties. Goldberg et al. (2019) demonstrate that Equation (1) accurately predicts the precision of transit mid-time measurements when transit S/Ns are sufficiently high (S/N ≳ 3; see their Appendix A).

A central aim of this work is to predict the expected improvements to the detection of TTVs in TESS data (and the mass measurements that can be extracted from these data) under various extended mission scenarios. We consider four distinct 3 yr extensions of the TESS mission, which we summarize in Table 2. Our extended mission scenarios cover two possible camera configurations: the first, C₄, has camera 4 centered on the ecliptic pole as in the nominal mission. The second, C₃, has camera 3 centered on the ecliptic pole and provides a larger area of sky with multiple pointings, at the expense of coverage near the ecliptic equator. (Huang et al. (2018b) refer to this configuration as C3PO.) We also consider two possible extended mission pointing sequences: one in which TESS remains pointed in the northern ecliptic hemisphere for the entire extended mission (NNN) and one in which TESS alternates hemispheres each year after starting in the north (NSN).

We approximate the TESS sky footprint as a rectangular region spanning $24^\circ \ \times 96^\circ$ in ecliptic longitude and latitude in order to determine which planets are observed during the various extended missions considered.

We ignore any potential TTVs of planets that are not detected during the nominal mission but later discovered during the extended mission. We do not expect such planets to contribute significantly to TESS's TTV yield for the same reasons they are missed during the nominal mission: either their transits are of low S/N and therefore their transit mid-times will be measured with poor precision or their orbital periods are longer than their nominal mission observational time baseline and they will have only a few transits during an extended mission.

To synthesize TTVs, we use the analytic TTV model described below in Section 3. We account for interactions between all planets with period ratios of P'/P < 2.2 in a given system, both transiting and nontransiting. Planets' initial mean anomalies are selected randomly from a uniform distribution. For simplicity, we assume all planets lie in the same plane in our dynamical integrations since the small mutual inclinations typical of multiplanet systems (e.g., Fang & Margot 2012; Figueira et al. 2012; Fabrycky et al. 2014) have negligible influence on planets' TTVs (e.g., Hadden & Lithwick 2016). We integrate each simulation for 5 yr and produce the times of mid-transit for every transiting planet. In order to generate statistics we repeat 50 simulations for each system, drawing new random mean anomalies each time.

We also compute a set of synthetic TTVs using the TTVFast N-body code of Deck et al. (2014) which we analyze in Appendix A. These N-body-synthesized TTVs exhibit signs of strong resonant interactions much more frequently than is observed in the existing TTV planet population: mean motion resonances (MMRs) are rare among Kepler's population of multiplanet systems (Fabrycky et al. 2014) and very few Kepler TTV systems show the effects due to second- or higher-order resonances (Hadden & Lithwick 2017). The prevalence of these resonant effects in our synthesized TTVs likely owes to imperfect modeling of the joint distribution of planetary systems' masses, period ratios, and eccentricities, upon which these effects sensitively depend. While inferring a joint distribution of planet masses, period ratios, and eccentricities that successfully reproduces the prevalence of resonant effects in the observed population is a potentially fruitful avenue of future research, it is beyond the scope of this work. Instead, we opt to use the analytic model as a means of producing a TTV data set that is similar to the observed TTV sample: our analytic model captures the effects of proximity to first-order MMRs and such TTVs comprise the majority of significant TTVs in the observed planet population (e.g., Lithwick et al. 2012).

We expect our analysis with the analytic model to provide an reliable estimate of how many planets are amenable to TTV mass measurements even though it may not faithfully reproduce every planet's TTV signal. This is because measurements of the chopping signal (Nesvorný & Vokrouhlický 2014; Deck & Agol 2015; Hadden & Lithwick 2016), which our model accounts for, are usually necessary for deriving mass constraints from TTVs whether planets are near first-order MMRs, higher-order MMRs (see, e.g., Kepler-26 b and c in Hadden & Lithwick 2016), or even librating in MMR (Nesvorný & Vokrouhlický 2016). Our analytic TTV model will faithfully reproduce the chopping component of the TTV as part of the δt⁽⁰⁾ basis-function component (see Linial et al. 2018) even when the full TTV signal is not accurately modeled.

Planets' TTVs can, to a good approximation, be treated as linear combinations of basis functions (Hadden & Lithwick 2016; Linial et al. 2018). These basis functions depend only on planets' periods and are thus known when all dynamically relevant planets in a system transit. Specifically, the nth transit of the ith planet is given by

$$\begin{eqnarray}\begin{array}{rcl}{t}_{i}(n) & = & {{nP}}_{i}+{T}_{i}\\ & & +\,\displaystyle \sum _{j\ne i}{\mu }_{j}\left(\delta {t}_{i,j}^{(0)}(n)+({{ \mathcal Z }}_{i,j}\delta {t}_{i,j}^{(1)}(n)+c.c.)\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where P_i is the ith planet's period, T_i sets the time of the ith planet's first transit, μ_j is the planet–star mass ratio of the jth planet,

$$\begin{eqnarray}&&{{ \mathcal Z }}_{i,j}\approx \displaystyle \frac{{e}_{j}{e}^{i{\varpi }_{j}}-{e}_{i}{e}^{i{\varpi }_{i}}}{\sqrt{2}}\end{eqnarray} \tag{ 3 }$$

is the combined complex eccentricity of planets i and j as defined in Hadden & Lithwick (2016) where e_i and ϖ_i are the eccentricity and argument of pericenter, c.c. denotes the complex conjugate of the preceding term, and $\delta {t}_{i,j}^{(0)}$ and $\delta {t}_{i,j}^{(1)}$ are TTV basis functions that capture the TTV induced by the jth planet. It will frequently be convenient to refer to the real and imaginary components of $\delta {t}_{i,j}^{(1)}$ separately as $\delta {t}_{i,j}^{(1,x)}$ and $\delta {t}_{i,j}^{(1,y)}$, respectively. A least-squares fit of a planet's TTV yields one real and one complex amplitude, μ_j and ${\mu }_{j}{{ \mathcal Z }}_{i,j}$, for each perturbing planet.⁵

The TTV basis functions, $\delta {t}_{i,j}^{(0)}$, ${t}_{i,j}^{(1,x)}$ and ${t}_{i,j}^{(1,y)}$, can be computed using a perturbative solution of the planets motion (Hadden & Lithwick 2016) or extracted from numerical integrations (Linial et al. 2018). Higher-order terms in ${{ \mathcal Z }}_{i,j}$ can also be included as needed using the perturbative method described in Hadden & Lithwick (2016). Details on computing the TTV basis functions are given in Appendix B.

A least-squares fit of a set of transit times provides the joint distribution of the basis-function amplitudes $\mu ,\mu \mathrm{Re}[{ \mathcal Z }],$ and $\mu \mathrm{Im}[{ \mathcal Z }]$ as a multivariate normal distribution. The quantities of interest in most applications will instead be μ and $Z\equiv | { \mathcal Z }|$, who's joint distribution can be obtained from the integral

$$\begin{eqnarray}&&P(\mu ,Z)={\mu }^{2}Z{\int }_{0}^{2\pi }N(\mu ,\mu Z\cos \theta ,\mu Z\sin \theta )d\theta ,\end{eqnarray} \tag{ 4 }$$

where $N(\cdot )$ is the Gaussian joint probability density of the TTV amplitudes, $\mu ,\mu \mathrm{Re}[{ \mathcal Z }],$ and $\mu \mathrm{Im}[{ \mathcal Z }]$, derived from the least-squares fit and the prefactor of ${\mu }^{2}Z$ comes from the Jacobian of the variable transformation $(\mu ,\mu \mathrm{Re}[{ \mathcal Z }],\mu \mathrm{Im}[{ \mathcal Z }])\to (\mu ,Z,\cos \theta )$.

TTV parameter inference is frequently hampered by a degeneracy between planet masses and eccentricities (Lithwick et al. 2012). This is because, often, the basis functions δt⁽⁰⁾ and $\delta {t}^{(1,x)}$ are both approximately sinusoids with identical phases. Consequently, there is a strong anticorrelation between the inferred amplitudes of these two basis functions, μ' and $\mu ^{\prime} \mathrm{Re}[{ \mathcal Z }]$.

The TTVs of two or more planets are determined by fewer free parameters than assumed by a linear model approach. Consider the TTVs of an interacting pair of planets: the TTVs are determined by four free parameters (two planet–star mass ratios and two components of the planets' ${ \mathcal Z }$). However, the linear model includes six independent TTV amplitudes: two planet–star mass ratios and two complex numbers of the form $\mu { \mathcal Z }$. Simply fitting two linear models to a pair of planets' TTV ignores the fact that the same combined eccentricity components appear in the basis-function amplitudes of both the inner and outer planet. Therefore, the linear model approach will likely overestimate the uncertainties in planet parameter constraints in some situations.

The linear TTV model described above can be used to predict the expected precision of planet parameter constraints that can be derived from a set of transit timing observations. Let us rewrite Equation (2) for the transit times of the ith planet in matrix form as

$$\begin{eqnarray}&&\left(\begin{array}{c}{t}_{i}(1)\\ {t}_{i}(2)\\ \vdots \\ {t}_{i}({N}_{\mathrm{trans}.})\end{array}\right)={{\boldsymbol{M}}}^{(i)}\cdot \left(\begin{array}{c}{T}_{i}\\ {P}_{i}\\ {\mu }_{1}\\ {\mu }_{1}\mathrm{Re}{[{ \mathcal Z }]}_{i,1}\\ {\mu }_{1}\mathrm{Im}{[{ \mathcal Z }]}_{i,1}\\ \vdots \end{array}\right),\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{M}}}^{(i)}$ is a $[2+3{N}_{\mathrm{pert}.}]\times {N}_{\mathrm{trans}.}$ matrix for a planet with ${N}_{\mathrm{trans}.}$ transit observations perturbed by ${N}_{\mathrm{pert}.}$ companion planets. Let ${\sigma }_{{t}_{c},l}$ be the uncertainty in the planet's lth transit time observation. We will assume that errors in transit time observations are normally distributed and independent. If we define the design matrix as

$$\begin{eqnarray}&&{A}_{{lm}}^{(i)}={M}_{{lm}}^{(i)}/{\sigma }_{{t}_{c},l}\,,\end{eqnarray} \tag{ 6 }$$

then the $[2+3{N}_{\mathrm{pl}.}]\times [2+3{N}_{\mathrm{pl}.}]$ covariance matrix between all basis-function amplitudes of the planet's TTV is given by

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={[{({{\boldsymbol{A}}}^{(i)})}^{T}\cdot {{\boldsymbol{A}}}^{(i)}]}^{-1}\,.\end{eqnarray} \tag{ 7 }$$

Importantly, ${\boldsymbol{\Sigma }}$ does not depend on the values of the basis-function amplitudes. Only knowledge of the planets' periods and initial orbital phases, which come directly from observations, is required to determine the TTV basis functions, and hence ${\boldsymbol{\Sigma }}$. Therefore, it is possible to predict the expected precision of TTV amplitude measurements for a set of transit observations without knowledge of the masses and eccentricities of the planets in the system.

Figure 1 illustrates the application of the linear TTV model with an example TTV data set taken from our synthesized TESS sample. The top panel shows the TTV of a planet near a 2:1 MMR with an exterior $37{M}_{\oplus }$ perturber, along with the decomposition of the TTV into its constituent basis functions. This particular system is located in the continuous viewing zone of the northern ecliptic hemisphere and is observed at a thirty-minute cadence during the nominal mission and then receives two-minute cadence observations during the first and third year of an extended mission. In the bottom panel we plot the predicted fractional uncertainty in the basis-function amplitudes versus the observing baseline. The amplitude uncertainties are determined by computing appropriate diagonal elements of Σ in Equation (7). In this example, the additional observations provided by the extended mission lead to a measurement of the perturbing planets' mass with a fractional uncertainty of $\sim 40 \%$ by constraining the amplitude of the δt⁽⁰⁾ basis function plotted in blue in the top panel of Figure 1.

Zoom In Zoom Out Reset image size

After our initial selection of possible TTVs (described above), we attempt to identify multitranet systems with significant TTV signals, in which the TTVs can be attributed to interactions among the observed planets. Such systems are generally the richest targets for TTV analysis since mass measurements can yield density constraints for the observed planets. To identify these systems, we fit the transit times of every multitranet TTV system in our catalog via least-squares with the linear model described in Section 3. Interactions between planets separated by a period ratio greater than $P^{\prime} /P\gt 2.2$ are ignored in the least-squares fits. These linear models account only for those planets that are observed to transit and ignore any potential nontransiting or undetected perturbers. We do not fit the TTVs of planets for which our analytic TTV model is underdetermined due to insufficient transits.⁷ After fitting a multitranet system's TTVs with our linear model, it falls into one of three categories.

• 1.  
  Category 1 systems yield direct mass measurements. The system hosts at least one planet for which a δt⁽⁰⁾ TTV component amplitude is nonzero at >1σ confidence. This means that one or more planet's TTVs constrain the mass of its perturbing companion(s).
• 2.  
  Category 2 systems exhibit significant TTVs, but mass inference is hampered by a mass-eccentricity degeneracy. For one or more planets in the system, the error ellipse defined by best-fit TTV amplitudes and their covariance matrix excludes zero at >2σ confidence. However, all δt⁽⁰⁾ amplitudes are consistent with zero, so the planets' masses are poorly constrained and strongly degenerate with the eccentricities.
• 3.  
  Category 3 systems show significant timing residuals but are poorly fit by the analytic model. The planets timing residuals, after subtracting off the best-fit analytic model, indicate a poor fit which we define as a χ² value inconsistent with random Gaussian noise at >3σ significance. For TTVs synthesized with the analytic model, poor fits are caused either by the presence of unseen additional perturbers or where the analytic model is underdetermined due to an insufficient number of transits. For the TTVs synthesized via N-body simulation analyzed in Appendix A, this category also includes TTVs exhibiting effects not captured by the analytic TTV model (such as proximity to second- or higher-order resonances). In some instances of this latter case, it may be possible to extract dynamical mass constraints from the TTVs.

Our categories are mutually exclusive, meaning that systems satisfying the criteria of multiple categories simultaneously are included in only one of them. Any system satisfying the Category 3 criterion is automatically classified as Category 3 irrespective of whether it also satisfies the criteria of either of the other two categories. In principle, a TTV system may or may not simultaneously meet the criteria for classification in Categories 1 and 2: a TTV can have δt⁽⁰⁾ amplitude detected with a >1σ significance (Category 1) independent of whether the total TTV signal is nonzero at >2σ confidence (Category 2). We find that, in practice, nearly every Category 1 TTV also satisfies the criterion for inclusion in Category 2 and we classify such systems as Category 1. Additionally, all Category 1 and 2 systems exhibit >3σ timing residuals relative to a linear ephemeris, and would therefore be included in Category 3 if their interactions were not modeled by the analytic model (though this is not an explicit requirement of our categorization scheme).

Figure 2 shows representative TTVs from the first two categories. For planets in Category 1, the TTV will provide a measurement of the perturbing companion's mass as well as the planet pair's combined eccentricity. Planets in Category 2 exhibit a degeneracy between the perturber mass and eccentricity. An independent measurement of the perturber's mass, e.g., via radial velocity, can yield eccentricity constraints in such systems.

Zoom In Zoom Out Reset image size

Figure 3 plots the masses and radii of the TTV detections from one of the 50 realizations of our synthetic TTV catalog. The left-hand panel shows the results of the nominal mission: eight mass measurements derived from Category 1 TTVs as well as five Category 2 TTVs that provide a degenerate joint mass-eccentricity constraint similar to the example in the right-hand panel of Figure 2. The masses and radii of all additional detected planets in multitransiting systems are plotted in gray.

Zoom In Zoom Out Reset image size

The right-hand panel of Figure 3 shows the mass–radius results for the same realization of the TTV catalog after three additional years of an E_4,NSN extended mission. The extended mission yields six additional Category 1 mass measurements and improves the precision of the mass measurements obtained during the nominal mission. The extended mission also adds 12 new Category 2 planets along with four Category 3 systems. All but two of the detected TTVs are induced by planets smaller than 4R_⊕.

Table 3 summarizes the predicted yields of detected TTV systems, separated into the above classifications, in both the nominal and extended mission scenarios.⁸ We list the median, minimum, and maximum number of planets falling within each category from the 50 simulation iterations randomized in planets' initial orbital phases. Overall, the various extended mission scenarios yield similar numbers of TTV detections.

In addition to the multitransiting TTV systems described above, TESS is expected to discover a number of singly transiting planets that exhibit TTVs. While such TTVs generally offer little information beyond an indication of additional planets in the system, it may be possible to extract dynamical constraints in some instances (e.g., Nesvorný et al. 2013). Table 3 summarized the number of such systems whose transit timing residuals are inconsistent with Gaussian random noise at >3σ for various mission scenarios. In all extended mission scenarios, the expected number of singly transiting TTV planets outnumbers the total combined multitranet systems in Categories 1–3.

Many of the multiplanet systems hosting single-tranet TTV planets contain additional transiting planets whose transits fall below TESS's 7.3 S/N detection criterion (e.g., Sullivan et al. 2015). Figure 4 shows the transit S/Ns of undetected transiting companions to the single-tranet TTV planets discovered in various mission scenarios. Many such companions could be recoverable upon relaxing the 7.3 transit S/N requirement. If these planets are responsible for the TTV of the initially detected planet then the TTVs could yield dynamical constraints on the mass and eccentricities of these planets.

Zoom In Zoom Out Reset image size

Table 4 summarizes the planet and tranet multiplicities of systems in the low-i population. As expected, there are more multitransiting systems in the low-i population, with 23 (=166–143) more two-tranet systems and 4 (=7–3) more three-tranet systems. Table 5 summarizes TTV detections and dynamical constraints for the low-i population. The median number of multitransiting TTV detections of each type increases slightly for the low-i (with a corresponding decrease in the number of singly transiting TTV) though the ranges produced by randomizing over orbital phases largely overlap between our fiducial and low-i populations.

In Section 2.4 we synthesized planets' TTV using analytic formulas. We also have run simulations that instead use the TTVFast code of Deck et al. (2014) to generate transit times via N-body integration. Table 6 summarizes TTV detections and dynamical constraints derived from TTVs synthesized via N-body simulation. As shown by Table 6, the N-body sample contains a significant number of Category 3 systems poorly fit by the analytic model. All Category 3 TTVs generated with the analytic formulas are poorly fit due to the presence undetected perturbers inducing TTVs that are unaccounted for when fitting them with the analytic model. By contrast, N-body TTVs with significant residuals can be caused by undetected perturbers as well as failures of the analytic formulas to accurately model the planets' interactions. The analytic model fails to accurately capture TTVs of planets that are strongly affected by a second- or higher-order MMRs or librating in an MMR. Inspecting the Category 3 TTVs of the N-body-generated sample, we find that significant residuals are most frequently caused by these resonant effects.

We write the time dependence of the planet's complex eccentricity, z, and mean longitude, λ, of the planet as

$$\begin{eqnarray}&&z(t)={z}_{0}+\delta z(t)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \frac{2\pi }{P}(t-T)+\delta \lambda (t),\end{eqnarray} \tag{ 15 }$$

where $\delta z(t)$ and $\delta \lambda (t)$ represent the deviations from a purely Keplerian orbit induced by the perturber. The planet's TTV is then given as a function of time by

$$\begin{eqnarray}&&\delta t(t)=-\displaystyle \frac{P}{2\pi }\left[\delta \lambda (t)+i\delta z(t)-i\delta {z}^{* }(t)\right]+{ \mathcal O }({e}^{2})\,.\end{eqnarray} \tag{ 16 }$$

The equations of motion are

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}a}{{dt}}=2n^{\prime} \displaystyle \frac{\mu ^{\prime} }{\sqrt{\alpha }}\displaystyle \frac{\partial R}{\partial \lambda }\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\delta \lambda }{{dt}}=-\displaystyle \frac{3n^{\prime} }{2{\alpha }^{3/2}}\displaystyle \frac{\delta a}{a}-2n^{\prime} \mu ^{\prime} \sqrt{\alpha }\displaystyle \frac{\partial R}{\partial \alpha }\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dz}}{{dt}}=2{in}^{\prime} \displaystyle \frac{\mu ^{\prime} }{\sqrt{\alpha }}\displaystyle \frac{\partial R}{\partial {z}^{* }},\end{eqnarray} \tag{ 19 }$$

where $\alpha =a/a^{\prime}$ and R is the disturbing function,

$$\begin{eqnarray}&&R=a^{\prime} \left(\displaystyle \frac{1}{| r^{\prime} -r| }-\displaystyle \frac{r\cdot r^{\prime} }{| r^{\prime} {| }^{3}}\right).\end{eqnarray} \tag{ 20 }$$

The disturbing function expanded to the first order in eccentricity can be written as

$$\begin{eqnarray}\begin{array}{rcl}R & = & \displaystyle \frac{1}{\sqrt{1+{\alpha }^{2}-2\alpha \cos \psi }}-\alpha \cos \psi \\ & & +\,\left({\left.{z}^{* }\displaystyle \frac{\partial R}{\partial {z}^{* }}\right|}_{z,{z}^{* }=0}+z{\left.\displaystyle \frac{\partial R}{\partial z}\right|}_{z,{z}^{* }=0}\right)+{ \mathcal O }({e}^{2})\end{array}\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{\partial R}{\partial {z}^{* }}\right|}_{z,{z}^{* }=0} & = & \left(\displaystyle \frac{{\alpha }^{2}+\alpha \cos \psi +2i\alpha \sin \psi }{{\left(1+{\alpha }^{2}-2\alpha \cos \psi \right)}^{3/2}}\right.\\ & & +\ \left.\displaystyle \frac{\alpha }{2}\cos \psi +i\alpha \sin \psi \right)\exp \left[i\lambda \right]\end{array}\end{eqnarray} \tag{ 22 }$$

and ψ = λ'–λ.

We begin by solving for δλ(t) to zeroth order in eccentricity. Let us define ${\rm{\Delta }}n=n^{\prime} -n$ so that, ignoring the effect of planet–planet interactions, $\psi ={\rm{\Delta }}{nt}+{\psi }_{0}$, where ψ₀ is determined by initial conditions. Inserting Equation (21) into Equation (17), to zeroth order in eccentricity and first order in μ' we have

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta a(t)}{a} & = & \displaystyle \frac{2n^{\prime} \mu ^{\prime} }{\sqrt{\alpha }{\rm{\Delta }}n}\left\{\displaystyle \frac{1}{\sqrt{1+{\alpha }^{2}-2\alpha \cos ({\rm{\Delta }}{nt}+{\psi }_{0})}}\right.\\ & & -\ {\left.\alpha \cos ({\rm{\Delta }}{nt}+{\psi }_{0})\right\}}_{\mathrm{osc}.}\end{array}\end{eqnarray} \tag{ 23 }$$

where the braces indicate the oscillating part of the enclosed expression, i.e.,

$$\begin{eqnarray}&&{\left\{f(t)\right\}}_{\mathrm{osc}.}=f(t)-\displaystyle \frac{{\rm{\Delta }}n}{2\pi }{\int }_{0}^{2\pi /{\rm{\Delta }}n}f(t^{\prime} ){dt}^{\prime} .\end{eqnarray} \tag{ 24 }$$

Inserting Equation (23) into Equation (18) we derive

$$\begin{eqnarray}&&\delta \lambda (t)=\mu ^{\prime} {\left\{{\left(\displaystyle \frac{n^{\prime} }{{\rm{\Delta }}n}\right)}^{2}A({\rm{\Delta }}{nt}+{\psi }_{0})+\displaystyle \frac{n^{\prime} }{{\rm{\Delta }}n}B({\rm{\Delta }}{nt}+{\psi }_{0})\right\}}_{\mathrm{osc}.}\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&A(x)=\displaystyle \frac{3}{{\alpha }^{2}}\left[\displaystyle \frac{2}{(1-\alpha )}K\left(\left.\displaystyle \frac{x}{2}\right|m\right)-\alpha \sin \psi \right]\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{rcl}B(x) & = & -\displaystyle \frac{2\sqrt{\alpha }}{1-{\alpha }^{2}}\left[\left(1-{\alpha }^{-1}\right)E\left(\left.\displaystyle \frac{x}{2}\right|m\right)\right.\\ & & +\,(1+{\alpha }^{-1})K\left(\left.\displaystyle \frac{x}{2}\right|m\right)\\ & & +\,\left.\left(1-{\alpha }^{2}-\displaystyle \frac{2}{\sqrt{1+{\alpha }^{2}-2\alpha \cos \psi }}\right)\sin \psi \right],\end{array}\end{eqnarray} \tag{ 27 }$$

K and E are incomplete elliptic integrals of the first and second kind, respectively, and $m=-\tfrac{4\alpha }{{(1-\alpha )}^{2}}$. Note that the oscillating parts of elliptic functions are given by

$$\begin{eqnarray}&&{\left\{K(x| m)\right\}}_{\mathrm{osc}.}=K(x| m)-\displaystyle \frac{2x}{\pi }K(m)\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\left\{E(x| m)\right\}}_{\mathrm{osc}.}=E(x| m)-\displaystyle \frac{2x}{\pi }E(m),\end{eqnarray} \tag{ 29 }$$

where $K(m)$ and E(m) are complete elliptic integrals of the first and second kind, respectively. Agol et al. (2005) derive a similar closed-form expression for δλ in their TTV formulas (see also Malhotra 1993).

Inserting Equation (21) into Equation (19), the equation of motion for z(t) to zeroth order in eccentricity is given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dz}}{{dt}}=2{in}^{\prime} \displaystyle \frac{\mu ^{\prime} }{\sqrt{\alpha }}\left(\displaystyle \frac{{\alpha }^{2}+\alpha \cos \psi +2i\alpha \sin \psi }{{\left(1+{\alpha }^{2}-2\alpha \cos \psi \right)}^{3/2}}\right.\\ \quad +\left.\displaystyle \frac{\alpha }{2}\cos \psi +i\alpha \sin \psi \right)\exp \left[i\lambda \right]\,.\end{array}\end{eqnarray} \tag{ 30 }$$

Unlike δλ(t), Equation (30) cannot be integrated to produce a simple closed-form solution for δz(t). Rather than relying on a Fourier expansion of Equation (30), we simply integrate Equation (30) numerically after inserting ψ = Δnt+ψ₀ and λ = nt+λ₀ to obtain δz(t) to zeroth order in eccentricity and first order in μ'. Combining our numerical solution to Equation (30) with Equation (18), we obtain the ${t}_{i,j}^{(0)}$ via Equation (16).

Here we provide the TTV of a planet perturbed by an interior companion. The derivation is essentially the same as above so we do not reproduce all the steps but instead mention some minor modifications and quote the final results. First, for an outer planet the disturbing function is given by

$$\begin{eqnarray}&&R^{\prime} =a^{\prime} \left(\displaystyle \frac{1}{| r^{\prime} -r| }-\displaystyle \frac{r\cdot r^{\prime} }{| r{| }^{3}}\right)\,\end{eqnarray} \tag{ 31 }$$

and the equations of motion are

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}a^{\prime} }{{dt}}=2n^{\prime} \mu \displaystyle \frac{\partial R^{\prime} }{\partial \lambda ^{\prime} }\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\delta \lambda }{{dt}}=-\displaystyle \frac{3n^{\prime} }{2}\displaystyle \frac{\delta a^{\prime} }{a^{\prime} }-2n^{\prime} \mu \left(1+\displaystyle \frac{\partial }{\partial \alpha }\right)R^{\prime} .\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dz}^{\prime} }{{dt}}=2{in}^{\prime} \mu \displaystyle \frac{\partial R^{\prime} }{\partial {z}^{{\prime} * }}\end{eqnarray} \tag{ 34 }$$

The outer planet's disturbing function, Equation (31), expanded to first order in z' becomes

$$\begin{eqnarray*}\begin{array}{rcl}R^{\prime} & = & \left({\left(1+{\alpha }^{2}-2\alpha \cos \psi \right)}^{-1/2}-{\alpha }^{-2}\cos \psi \right)\\ & & +\,\left(\displaystyle \frac{\alpha \cos \psi -1+2i\alpha \sin \psi }{2{\left(1+{\alpha }^{2}-2\alpha \cos \psi \right)}^{3/2}}+\displaystyle \frac{1}{2{\alpha }^{2}}\left(\cos \psi -2i\sin \psi \right)\right)\\ & & \times \,{z}^{{\prime} * }{e}^{i\lambda ^{\prime} }+c.c.\,.\end{array}\end{eqnarray*}$$

The solution for δλ' is given by

$$\begin{eqnarray}\begin{array}{rcl}\delta \lambda ^{\prime} (t) & = & -\mu \left\{{\left(\displaystyle \frac{n^{\prime} }{{\rm{\Delta }}n}\right)}^{2}A^{\prime} ({\rm{\Delta }}{nt}+{\psi }_{0})\right.\\ & & {\left.+\displaystyle \frac{n^{\prime} }{{\rm{\Delta }}n}B^{\prime} ({\rm{\Delta }}{nt}+{\psi }_{0})\right\}}_{\mathrm{osc}.}\end{array}\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&A^{\prime} (\psi )=-{\alpha }^{2}A(\psi )+3({\alpha }^{-2}-\alpha )\sin \psi \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&B^{\prime} (\psi )=-\sqrt{\alpha }B(\psi )-\displaystyle \frac{2{\alpha }^{2}}{3}A(\psi )-(4\alpha +2{\alpha }^{-2})\sin \psi \end{eqnarray} \tag{ 37 }$$

and the solution for $\delta z^{\prime} (t)$ is obtained by numerically integrating

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dz}^{\prime} }{{dt}}=2{in}^{\prime} \mu \left(\displaystyle \frac{\alpha \cos \psi -1+2i\alpha \sin \psi }{2{\left(1+{\alpha }^{2}-2\alpha \cos \psi \right)}^{3/2}}\right.\\ \quad \left.+\displaystyle \frac{1}{2{\alpha }^{2}}\left(\cos \psi -2i\sin \psi \right)\right)\exp \left[i\lambda ^{\prime} \right]\,.\end{array}\end{eqnarray} \tag{ 38 }$$

with respect to time after substituting ψ = Δnt + ψ₀ and $\lambda ^{\prime} =n^{\prime} t+{\lambda }_{0}^{{\prime} }$.