QUANTIFYING THE CHALLENGES OF DETECTING UNSEEN PLANETARY COMPANIONS WITH TRANSIT TIMING VARIATIONS

1.1. Motivation

1.2. Observational Studies

1.3. Theoretical Studies

1.4. Our Methodology

As a first step toward understanding the large parameter space for two-planet systems, we establish a fiducial set of simulations which broadly characterize the TTV signals produced from a transiting hot Jupiter and an external terrestrial-mass planet. We consider an M_i = 1 M_J inner planet which transits its parent M_⋆ = 1 M_☉ star on an initially circular orbit (e_i = 0) at a_i = 0.05 AU, and an external M_o = 1 M_⊕ planet with the semimajor axis, a_o, and eccentricity, e_o. The subscripts "i" and "o" are abbreviations for "inner" and "outer." The mean longitude, mean anomaly, and longitude of pericenter will be denoted by λ, Π, and ϖ, respectively. For any system, only one of λ or Π needs to be specified; either determines a planet's location along a Keplerian orbit, and the former is often used in planar resonance studies (Murray & Dermott 2000). Henceforth, all stated orbital parameters are assumed to be initial values unless an explicit time dependence is included.

After establishing values for a_o/a_i, e_o, ϖ_o, Π_i or λ_i, and Π_o or λ_o, we integrate the system over 10 yr and tabulate the times of transits, and subtract the best-fit Keplerian model to generate a TTV curve. The deviations of the transit times from strict periodicity tend to form periodic patterns which can have timescales extending to the duration of the simulations. We display a sample of these patterns and periodicities in Figure 1 over the course of the nominal Kepler mission lifetime (3.5 yr) in order to motivate the other figures in this paper. The dots represent transit times, and the amplitude of their variations ranges over three orders of magnitude across the panels. The variety of the patterns seen hints at the difficulties in deriving masses and orbital parameters of the external planet from TTV curves alone. The top panel curve (a_o/a_i = 3.6593, e_o = 0.596, Π_o = 70°, Π_i = 0°, ϖ_o = 180°) may be fitted well with a single sinusoid. The curve in the next panel below (a_o/a_i = 2.3313, e_o = 0.395, Π_i = Π_o = 0°, ϖ_o = 0°) exhibits several modulations in addition to a simple sinusoid. Up until ≈1 yr, and between ≈2.5and3.5 yr, this curve appears to exhibit a long period (>3.5 yr) trend with a modulated sinusoid of amplitude less than 10 s. In between, at ≈1and2.5 yr, the curve appears to oscillate with an amplitude of 20 s without modulation. The curve in the second panel from bottom (a_o/a_i = 1.5812, e_o = 0.204, Π_i = Π_o = 0°, ϖ_o = 0°) appears to exhibit a very long periodic trend (⩾3.5 yr) with a high amplitude (hundreds of seconds) and a sawtooth-like modulation of amplitude <100 s. Note the distribution of spacings between the transit times on the modulated sinusoids. These spacings may provide hints as to the architecture of the system (see Section 7). The bottom panel (a_o/a_i = 1.5812, e_o = 0.219, Π_i = Π_o = 0°, ϖ_o = 0°) displays TTVs which vary by thousands of seconds over the first ≈2 yr before abruptly taking on the form of a simple sinusoid with amplitude <10³ s. This panel especially highlights the importance of considering a sufficient number of observed transits with any type of analysis in order to avoid spurious conclusions. We emphasize that the four curves in Figure 1 are not representative of all curves henceforth studied, but rather illustrate different types of trends and patterns that one might uncover.

Zoom In Zoom Out Reset image size

In principle, we may select any combination of transits to mimic observed data. For our fiducial case, we select the TTV curves based on the first 10 yr of transits (N ≈ 874, where N is defined as the number of transits). Figure 2 displays the results of integrations of 120,000 systems. We sampled 400 logarithmically spaced values of the semimajor axis ratio and 60 uniformly spaced values of e_o for each semimajor axis ratio. For every given pair of a_o/a_i and e_o values, we simulated five systems initialized with random mean anomalies of both planets and a random longitude of pericenter of the outer planet. The reported values represent the median RMS TTV deviation amplitude (in seconds) of each of these sets of five simulations. We only sampled systems guaranteed to be stable according to the Hill stability limit (Gladman 1993), i.e., where e_o < e_H. The value of e_H is a function of both planets' masses, semimajor axes, and eccentricities, and for e_o < e_H, the orbits are guaranteed to never cross. As we vary these parameters, the bounding curve (where e_o = e_H) in the figure will change slightly. We caution that (1) some (unsampled) systems which lie above the curve in Figure 1 may be stable, and that (2) some systems close to the Hill stability boundary which we do sample might be Lagrange unstable (where the outer planet drifts outward, generally causing an increase in the RMS TTV signal). We explore the extent of the systems that are likely to be Lagrange unstable in Section 8.1, but note briefly here that such systems are very unlikely to be observed. The inner semimajor axis ratio bound of 1.3 was chosen to roughly correspond to the point where the Hill stability curve intersects the x-axis. The outer semimajor axis ratio bound of 5.0 is arbitrarily chosen to allow one to consider highly hierarchical (widely separated) systems. Additionally, 5.0 is the ratio beyond which analytic formulas achieve success at reproducing the RMS TTV amplitude given the orbital parameters of the exterior planet (Figure 2 of Agol et al. 2005).

Zoom In Zoom Out Reset image size

Figure 2 is contoured on a logarithmic scale. If one takes 10 s as the current detectability threshold for TTV signals, then any region on the plot that is not white, pink, or red should contain a detectable signal. Hence, widely separated (a_o/a_i > 3) planets containing an external planet on a low-eccentricity orbit (e₁ ≲ 0.2) produce currently undetectable signals. Additionally, for almost any a_o/a_i ⩽ 5 with an external planet whose orbital eccentricity is close to the Hill stability limit, the resulting TTV signal will be easily (⩾10³ s) detectable. The figure demonstrates that the signal amplitude is highly sensitive (on the 10⁻³–10⁻² AU scale) to the semimajor axis ratio.

The flame-like features on the plot, similar to those found in Figure 5 of Agol et al. (2005), indicate regions where the two planets are near a period commensurability (PC) that can be expressed as a ratio between two small integers (p:q, where p, q ⩽ 20). The longest of these "flames of resonance" corresponds to PCs of the form p:1, where p ⩽ 11. These PC locations well approximate the regions where both planets might be locked in or reside just outside of an MMR. Whether or not the planets with a given p:q ratio are in an MMR is subject to individual system study, and is often highly dependent on its orbital angle architecture. In the exoplanet dynamics community, a widely used criterion for determining whether a system is "in" MMR is to determine if at least one "resonant angle" is librating. Consideration of multiple resonant angles for a given p:q can be found in Figure 3 of Laughlin & Chambers (2001), Section 4.2 of Kley et al. (2005), Section 4.1.1 of Raymond et al. (2008), Figure 9 of Crida et al. (2008), and Section 5 of Fabrycky & Murray-Clay (2010). However, Michtchenko et al. (2008a, 2008b) claim that for a given p:q, an MMR is characterized just by one resonant angle, whereas in the planar case, the other independent angular variable should be a secular angle. Mardling (2008, 2010) asserts that for a given p:q, the resonant angle can be expressed in terms of a single index that traces the contribution from a particular semimajor axis ratio order in the gravitational potential.

For our purposes, we can neglect the contribution of the inner planet's longitude of pericenter (ϖ_i) because e_i = 0 and the mass ratio of the outer planet to inner planet (M_o/M_i) is negligible (∼0.001). Hence, e_i(t) will remain low (⩽0.01) as t increases, and the only resonant angle we consider is

$$\begin{equation} \phi _{p,q} \equiv p \lambda _o - q \lambda _i - (p-q) \varpi _o. \end{equation} \tag{ 1 }$$

The libration width of this resonance, defined here to be a range of a_o at a given e_o, bounds a region of phase space where this angle might librate. We use the formalism in Murray & Dermott (2000) to compute libration widths for a sampling of first- to fourth-order resonances, where order ≡p − q. S. F. Dermott (2010, private communication) also provided us with the coefficients needed to compute libration widths for fifth- to eighth-order resonant libration widths, with q = 1. The horizontal spacing between the thin vertical lines extending from the x-axis on the plots in Figure 2 represents these libration widths. Note that the spacing is negligible until e_o ≳ e_H/2. Furthermore, the lines do not extend to the Hill stability boundary because of the Sundman convergence criterion (Ferraz-Mello 1994; Sidlichovsky & Nesvorny 1994). See the Appendix for additional details and a derivation of this libration width.

The location and extent of the libration widths help confirm that pronounced RMS TTV amplitudes (≡S(Q) ≡ S(t, a_o/a_i, e_o, Π_o, Π_i, ϖ_o)) are likely caused by PCs and that MMR configurations are achieved only for a high enough eccentricity of the outer planet. Additionally, the libration widths for the 9:1 MMR demonstrate that two planets may be locked in that resonance despite the weakness of an eighth-order resonance (resonant strength scales as ∼e^p−q). Such planets would likely produce a detectable TTV signal. Also, note that the length of each "flame" corresponds to the value of q. Planets locked in a high-order (⩾10) MMR with q = 1 or 2 can produce a detectable signal (S(Q)>10s) for sufficiently large e_o. Although the highest order MMR which would produce a detectable signal is highly dependent on the characteristics of the systems studied, planets in MMR of order 20–30 can in principle produce detectable signals if e_o is high enough. Thus, a transiting hot Jupiter at 0.05 AU could have detectable TTVs due to an Earth-mass planet in the habitable zone of a K- or M-type star. Because S(Q) is highly sensitive (on the 10⁻³–10⁻² AU scale) to the semimajor axis ratio, detectable MMRs may help identify planets with particular orbital parameters. However, when considering high q MMRs, one should also take into account the time dependence of resonant locations; as a system evolves, libration widths would likely oscillate in a_o(t)–e_o(t) space.

Figure 2 displays the median RMS TTV deviation amplitude (in seconds) of five different orbital angle configurations, whereas Figure 3 reports the minimum and maximum deviations. We emphasize that the following conclusions are for our fiducial case, with a given perturbing mass of 1M_⊕. The differences in the two plots in Figure 3 are striking. For a given a_o and e_o, S(Q) may vary by orders of magnitude based on the initial values of the mean longitudes/mean anomalies and the outer planet longitude of pericenter (ϖ_i is meaningless when e_i = 0). The range of TTV signal amplitudes achieved from the different combinations of initial angles is comparable to the lower panel of Figure 3, the maximum signal achieved. Furthermore, the median signals are on average greater than the mean signals (not shown). The secular regime, away from PC, demonstrates little variation (at most a few seconds) in signal with initial orbital configuration. However, detections in these regimes are unlikely without the advent of better precision observational techniques. The upper panel in Figure 3, which displays the minimum signal achieved, best shows that for many strong (order ⩽7) MMRs with q = 1 or 2, the RMS signal is lower than that from the high-eccentricity (at e_H) secular regime. Therefore, for observations of a hot Jupiter, (1) moderate (tens of seconds) values of S(Q) could indicate that e_o may be as high as e_H, (2) high amplitude (≳10³) signals are not necessarily indicative of close proximity (within a few libration widths) to a strong PC, and (3) high amplitude (≳10³) S(Q) values are always indicative of a high value (≳0.8e_H) of e_o.

Zoom In Zoom Out Reset image size

The choice of contour levels is arbitrary, but meant to show phase space structure and establish the detectability threshold of 10s. Increasing each contour level by one order of magnitude indicates that little structure can be discerned in the highest signal (purple) regions of Figures 2 and 3 and that S(Q) can reach well over 10⁴ s. In contrast, lowering each contour level by one order of magnitude demonstrates a clear structure in the lowest signal regime (pink) of Figures 2 and 3. TTV curves in this regime have been successfully correlated with the mass and orbital elements of a hypothetical planet (Agol et al. 2005; Nesvorný & Morbidelli 2008). However, characteristic exoplanet detection thresholds (≈10 s) exceed typical signal amplitudes in this regime.

TTV signal amplitudes, S(Q), crucially depend on the number of observed transits, N, and the sampling rate. Although one can define TTVs solely as the difference between consecutive transit times, Agol et al. (2005) define TTVs as the deviation from an overall linear fit to the transit times of a system. This definition does not rely on having successive transits. Before missions such as CoRoT and Kepler, ground-based observations have struggled to yield N ⩾ 100. These already operating space-based missions, however, hope to observe consecutive transits for a half year (CoRoT) or multiple years (Kepler). Therefore, we consider our fiducial system (Figure 2; we henceforth use the median value when averaging over multiple random angle initial configurations) for various durations corresponding to N = 874, 313, 100, 50, 30, and 10 (Figures 4–6) and by considering consecutive transits. The first two values of N listed correspond to observing campaigns of approximately 10 yr and 3.5 yr (roughly the maximum and nominal lifetime of Kepler).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 illustrates the TTV amplitude variation as a function of N. For the N = 10 plot, almost no systems are detectable. At N = 30, only a few systems featuring high eccentricity may be detectable. For N ⩾ 100, systems close to a q = 1 PC generate signals which change little with N for e_o < 0.75e_H. However, S(Q) continues to increase for increasing N for e_o values closer to e_H.

Because S(Q) is sensitive to the semimajor axis ratio on scales potentially smaller than 10⁻³ AU, higher resolution sampling of phase space will reveal additional features. Similarly, as N increases, larger sections of a_o − e_o phase space will exhibit TTV curves with clearly discernible periods. Hence, we have performed high-resolution (200 points on the x-axis and 120 points on the y-axis) explorations of the phase space around the 3:1 PC in Figure 5. The figure demonstrates that as N increases, the contour outlines become sharper and other weaker commensurabilities which neighbor the 3:1 PC become more apparent. Additionally, the double-lobed feature on the 3:1 PC location at e_o = 0.1 disappears with increasing N.

Importantly, the figure demonstrates that planets in a PC, and possibly in an MMR, do produce a distinct TTV signature, but not necessarily a high amplitude signature compared to its near-PC surroundings. At N = 30, S(Q) at PC (yellow) is lower than the signal near PC (light green). At N = 313, S(Q) at PC (light green) is lower than the signal near PC (olive), and this result is insensitive to further increases of N. The near-PC regime is sharply divided from the secular regime (note the thin light green contour between the regimes). Additionally, after 10 yr, there is a sharp boundary (marked by sparse blue dots) between the near-PC regime and the high-eccentricity (at ≈e_H) secular regime. Ten years might not be long enough for the some TTVs to produce periodic signals in the blue or purple areas.

The variation of TTV signal with N is not necessarily monotonic, and is a function of the masses and orbital parameters of the planets. We can look more closely at this S(Q) versus N dependence by fixing all other orbital parameters in specific cases. Although identifying "representative" systems is difficult because of the large size of the phase space, we sample eight different particular systems in each of three regimes: "in-PC," "near-PC," and "secular." Table 1 contains the a_o/a_i and e_o values for all 24 systems, and Figure 6 plots S(Q) versus N for these systems.

We first consider the "in-PC" systems plotted in the upper plot of Figure 6. In order to model the signal amplitude of any in-PC system with continuous observational coverage, one must observe at least 50 transits. Only for N > 50 do any of the curves level off. Generally, the higher the amplitude of S(Q), the greater the number of observations needed to level the curves (for a given perturber mass). Only the red (4:1), salmon (11:1), and purple (11:3) curves are level at N = 874, and only for the blue (2:1) curve is S(Q) decreasing. The signal curves are generally non-monotonic with N; additional observations may cause the signal amplitude to increase or decrease. Note additionally that although eccentricities are high (0.47e_H ⩽ e_o ⩽ 0.80e_H), the signal amplitudes are relatively low (S(Q) < 200 s) for all resonant curves for all t ⩽ 10 yr. Systems in a PC do not necessarily achieve signal amplitudes as high as those in other regimes, and preferentially have lower amplitudes, even after 10 yr of observational sampling.

The prospects for detecting systems near but not on PC is more promising. The middle panel of Figure 6 indicates that most of the curves level off after N ≈ 300 (a little over 3 yr). However, at N = 800, one can see the aqua (5:2) and olive (3:1) curves beginning to trend upwards. The system near the 2:1 PC (blue) exhibits similar oscillatory behavior from the 2:1 in-PC system in the upper panel. The two curves corresponding to the highest e_o values (salmon with e_o = 0.85e_H and magenta with e_o = 0.80e_H) maintain S(Q) < 20 s for all N. The salmon curves in the upper two panels of Figure 6 contain the same value of e_o and both level at 20 s. However, the in-PC system levels off N ≈ 230, whereas the near-PC system levels off at N ≈ 110, reflecting the general trend of near-PC systems having robust values of S(Q) with fewer observations.

The two secular systems with the greatest signals in the lower panel of Figure 6 (olive and magenta) undergo variations on the order of their initial values for at least 5 yr. Both of these systems feature e_o > 0.80e_H and a_o/a_i ≲ 2.2. The strongly hierarchical (a_o/a_i = 4.405) system indicated by the gray curve raises its TTV signal by nearly one order of magnitude between N = 250and300, and then doubles this value over the next 6 yr. The variation with N of all three of these systems dwarfs the variation by the other five systems on the plot, including two with high eccentricities (aqua and purple, with e_o = 0.73e_H and 0.71e_H, respectively). These two high eccentricity systems show little (a few percent) variation for N > 200. Additionally, the signal for two of the lower eccentricity systems indicated by the blue and red curves (e_o = 0.14e_H and 0.25e_H) appears to show negligible variation with N for 50 ≲ N ≲ 400. However, both signals begin to increase for N ≳ 400 such that S(N = 874)/S(N = 400)>1.5. The qualitative difference in all these secular curves helps illustrate that S(Q) is non-trivially dependent on N, even for secular systems far from a PC.

Although we have already shown in Figure 2 and Section 2 that how randomly chosen sets of initial orbital angles (mean anomalies and longitudes of pericenter) for a given a_o and e_o can produce order of magnitude variations in TTV signals, here we study this dependency systematically. In the following analyses, we keep all orbital parameters but one fixed in order to gain insight into the physical origin of the TTV signals.

Although such insight is difficult to discern when considering how the entire (a_o, e_o) phase space varies with individual orbital angles, line plots for fixed (a_o, e_o) values show more structure. First, however, we consider "flame" plots as a function of Π_o. In each plot of Figure 7, the initial orbital angles are Π_i = ϖ_o = 0°, with Π_o = 0°, 10°, ⋅ ⋅ ⋅, 80°. Systems with planets initially at conjunction have high signal amplitude at high (near e_H) eccentricity, and, unlike in all other plots, saturate PCs at these eccentricities with S(Q)>10³ s at N = 874. For Π_o = 10°–80°, in each plot, several PCs produce signals which are over two orders of magnitude lower than 10³ s at high e_o. The "flames" which produce these effects appear at different sets of PCs in each plot, and have no immediately recognizable pattern. The lack of an apparent pattern in these plots, despite fixing all but one parameter and despite the high value of N = 874, expresses well the difficulty in trying to use transit timing variations to deduce the properties of unseen planets.

Zoom In Zoom Out Reset image size

An alternative and perhaps more demonstrative method of exhibiting this signal decrease at high eccentricities is through sequences of individual transit curves. Figure 8 displays nine sequences (panels) of scaled curves, where each curve traces 3.5 yr, as e_o is increased from 0 to e_H at the 7:1 PC (a_o/a_i = 3.66). In each panel, e_o is increased from the top left curve to the bottom left, to the upper right to the bottom right. The numbers accompanying each curve represents the RMS TTV amplitude of that curve in seconds. The curves are scaled such that their extrema are the plot boundaries, so that one could better see detail and modulation.

Zoom In Zoom Out Reset image size

The figure demonstrates that at this PC: (1) although low-eccentricity curves (e_o ∼ 0) are highly structured and modulated, their signal amplitude is negligible, (2) at mid-to-high (e_o ∼ 0.7e_H) eccentricity, the frequency of the curves suddenly increases, (3) in only some cases, eccentricities closest to e_H will exhibit long-period (>3.5) yr signals of ⩾10³ s, (4) the transitions to these high signals are often sudden (caused by a difference of e ≲ 0.018), and (5) in no case does the signal monotonically increase with increasing e_o. In particular, for Π_o = 10°, 40°, and 60°, TTV curves exhibit low (<10 s) amplitudes at e_H. For the other values of Π_o, TTV signals at these high eccentricities exhibit long (>3.5 yr) period variations and very high (>10³ s) amplitudes. This qualitative difference might vanish at high enough eccentricities, beyond the Hill stability limit. Note additionally that the eccentricity range which allows for high-frequency (featuring several crests over 3.5 yr) TTV signals varies significantly (over a factor of two) for different Π_o values. For Π_o = 50°, the long-trend signals near e_H are punctured by an orderly, periodic 89 s signal amplitude, showcasing the unpredictability of TTV signals near MMR.

The above analysis keeps ϖ_o fixed. Although this angle refers to the outer planet's orbit and not the planet's location on that orbit, the value of ϖ_o may play a crucial role in the dynamics. In fact, systems with ϖ_o = 90° and ϖ_o = 180° (not plotted) do show completely different patterns of low amplitude (<100s) "flames" near the Hill stability limit from those shown in Figure 7. The starkest difference appears in the Π_o = 0° case, where these flames are absent for ϖ_o = 0°, fully saturate the Hill stability region for ϖ_o = 90°, and appear only for a_o/a_i ⩽ 2.5 when ϖ_o = 180°.

The above analysis does not necessarily extend to other PCs. As the coefficients of the terms in disturbing function, or the different shapes of the libration widths in Figure 2 would indicate, the resonant structure of each PC is different. Analysis of planets thought to be near or in a particular MMR hence would benefit from a high resolution and complete exploration of the phase space at that location.

In order to better determine if the variations with Π_o in Figure 7 have structure, we now consider the signal variations while fixing (a_o, e_o). We plot this variation for selected systems (I, II, VII, XVII, XX, XXII from Table 1) in Figure 9, where the solid, dotted, and dashed curves in each plot correspond to ϖ_o = 0°, 90°, and 180°. The y-axes of all plots are logarithmic. At the coarse resolution of 10° per data point, the curves in the upper panels are vaguely oscillatory. For both the 2:1 (I) and 3:1 (II) systems, the S(Q) sharply (by several hundred seconds over 10°) peaks at Π_o = 90°, but for a different value of ϖ_o at each PC. At the 7:2 PC, this peak does not occur; curves for the three values of ϖ_o sampled have S(Q) with N = 874 values all within 10 s of one another at Π_o = 90°. All three upper panels indicate an anticorrelation for the blue and olive curves, with the most pronounced difference (over 400 s) for the 3:1 PC at Π_o = 90°.

Zoom In Zoom Out Reset image size

The variation of S(Q) with Π_o and ϖ_o for the secular systems sampled (lower panels of Figure 9) qualitatively differs from those for the in-PC systems. For the low-eccentricity system XVII (e_o = 0.14e_H), the curves have fewer extrema and have S(Q) < 10 s except for the dotted and dashed curves peaking at 30° and 40°. The other two secular systems sampled are at high (e_o ⩾ 0.73e_H) eccentricity, and both feature curves with multiple signal extrema, despite the 3–4 orders of magnitude difference in S(Q) between those two plots. The greatest extrema occurs at 60° for the dotted curve in system XX and at 70° for the solid curve in system XXII. The apparent anticorrelation between the solid and dotted curves from the in-PC systems appears to vanish for these secular systems.

In order to observe the variation of S(Q) with Π_o at a higher resolution of Π_o, we performed additional simulations. For the secular systems XVII, XVIII, and XIX, we varied Π_o over the entire [0°, 360°] range by sampling the angle at intervals of 01, and report the results in Figure 10. The three panels show qualitatively different behavior. The bottom panel (with S(Q) < 1 s) features smooth curves, the top panel curves (with 2.8 s <S(Q) < 4.4 s) show spikes and dips which could be hidden in broader resolution studies, and the middle panel (with S(Q) up to 7 × 10³ s) shows little structure. The figure demonstrates that TTV amplitude can vary by seconds due to changes of a few degrees in Π_o in some regimes, and vary by thousands of seconds due to changes of a fraction of a degree in Π_o in others.

Zoom In Zoom Out Reset image size

We plot these curves primarily to demonstrate the difficulties inherent in solving for the mass and orbital parameters of the unseen planet. Because of the broad parameter space, one is hard-pressed to determine "representative" systems in each regime. Additionally, because of the semimajor axis resolution of our simulations, the typical difference between nominal MMR a_o values and those from our simulations is ∼10⁻³–10⁻⁴. This difference may be significant due to the sensitivity of TTV signal profiles on the planetary semimajor axis ratio. Thus, we conclude that any TTV signal will need individual analysis, as opposed to being characterized by a few summary statistics. Furthermore, limits on the mass of the planets possible for a given TTV data set must be mindful of the possibility that a putative planet could have orbital elements that result in a signal much smaller than is typical for a planet of a given mass, semimajor axis, and eccentricity.

Thus far, we have fixed the masses of the star and both planets in all systems studied. Because these three masses are hierarchical (M_☉ ≫ M_J ≫ M_⊕), varying any one of them by a factor of a few would not significantly alter the contour phase space structure of Figure 2. However, mass variation might qualitatively affect individual systems at the edge of a secular, near-PC, or in-PC regime, or with an outer planet at a moderate-to-high (>0.5e_H) eccentricity.

We consider four different external perturber masses (1 M_⊕, 5 M_⊕, 10 M_⊕, and50 M_⊕) and five different transiting planet masses (0.1 M_J, 0.5 M_J, 1 M_J, 5 M_J, and10 M_J). Because Hill stability is a function of these masses, transiting inner planets more massive than 1 M_J will most significantly restrict the semimajor axis–eccentricity space in which the system is guaranteed to be stable. Therefore, in Figure 11, the phase space plotted is smaller than that of Figure 2 for all 20 combinations of planetary masses.

Zoom In Zoom Out Reset image size

Increasing the mass of the outer planet increases the TTV signal in almost all areas of phase space. More subtly, increasing mass causes the signal amplitude at resonance to be higher than that from near-resonance, in contrast to Figure 5. Fewer PC signatures are discernible at higher-mass-transiting planets; for M_i = 10 M_J, only PCs with q = 1 or 2 can be identified by inspection in the figure, independent of the value of M_o. Correspondingly, one can distinguish the greatest number of PCs for the lowest M_i values. Additionally, for the most massive interior and exterior planets, the highest signal stable (purple) region is greatest in extent.

Unlike the correlation between TTV signal and orbital angles, the correlation between TTV signal and planetary mass is robustly exponential in most stable regimes. Plotted in Figure 12 are six systems (I, IV, IX, XIV, XVII, and XVIII) from Table 1, with the solid, dotted, dashed, dot-dashed, and triple-dot–dashed curves corresponding to M_i = 0.1 M_J, 0.5 M_J, 1 M_J, 5 M_J, and 10 M_J. All six systems illustrate explicitly a roughly exponential increase in S(Q) as a function of M_o. The upper panel systems, which correspond to systems in the 2:1 and 8:1 PCs, demonstrate that S(Q) decreases with increasing M_i, unlike in the panels below, and that for 8:1, all curves monotonically increase, whereas the curves for 2:1 do not. The middle panel of the figure displays 2:1 and 5:2 near-PC systems with e_o = 0.44e_H and e_o = 0.49e_H, respectively. The difference in S(Q) upon varying the planetary masses is striking: six orders of magnitude for the 2:1 near-PC system, and just two orders of magnitude for the 5:2 near-PC system. The secular systems in the lowest panel show two anomalous features: (1) for system XVII, S(Q) for the M_i = 10 M_J curve is two orders of magnitude higher than any of the other M_i curves, and (2) system XXIII joins the in-PC 8:1 system as the only ones which contain a curve (solid; M_i = 0.1 M_J), where S(Q) decreases as M_o is increased from 1 M_⊕ to 50 M_⊕.

Zoom In Zoom Out Reset image size

Over the course of a multi-year observing campaign, the semimajor axis and eccentricity of the transiting planet might vary noticeably. N-body simulations show that an exterior terrestrial planet could induce the semimajor axis and eccentricity variations of ∼0.003 AU and ∼0.005, respectively, of the massive transiting Jupiter. If this variation is observed over time, then in principle these profiles could supplement TTVs as a method for identifying or confining the parameters of unseen planets.

Away from PC, although various formulations (Veras & Armitage 2007) of Laplace–Lagrange secular theory (see Murray & Dermott 2000) may be used to provide approximations of e_i(t) profiles which do not exceed low-to-moderate eccentricities ∼0.2–0.4e_H, these secular timescales (often 10^4–8 yr) exceed the duration of optimistic observing campaigns (10 yr) by many orders of magnitude. One should instead consider evolution on orbital timescales. One could apply perturbation theory without averaging over mean longitude in order to model these systems. Alternatively, suites of N-body simulations are prudent since the CPU time required to evolve two-planet systems on the timescale of 10 yr is on the order of μs.

At PC, and especially in MMR, theory may yield an additional useful constraint. Because of the mass hierarchy in our system, and the initially circular orbit of M_i, the resonant state of the systems studied here is likely to be dominated by a single resonant angle of the form in Equation (1). In this situation, constants of motion exist which relate just two orbital parameters to each other through conservation of angular momentum and energy. For any eccentricity type of resonance, one of these constants, C₁, may be expressed (Veras 2007) in terms of a_i(t) and e_i(t) only (recall that e_i ≡ e_i(t = 0) and a_i ≡ a_i(t = 0)):

$$\begin{eqnarray} C_1 & \equiv & \sqrt{a_i} \big[ p \big(\sqrt{1 - e_{i}^2} - 1 \big) - \big(p-q \big) \big] \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} & = & \sqrt{a_{i}(t)}[ p (\sqrt{1 - e_{i}(t)^2} - 1) - (p-q) ]. \end{eqnarray} \tag{ 3 }$$

Another constant, C₂, relates a_i(t) to a_o(t):

$$\begin{eqnarray} C_2 & \equiv & \sqrt{a_o} \left(\frac{q M_o}{p M_i} \right) \frac{M_{\star } \sqrt{M_{\star } + M_o + M_i} }{ \left(M_{\star } + M_1 \right)^{3/2} } - \sqrt{a_i} \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} & = & \sqrt{a_o(t)} \left(\frac{q M_o}{p M_i} \right) \frac{M_{\star } \sqrt{M_{\star } + M_o + M_i } }{ \left(M_{\star } + M_o \right)^{3/2} }\nonumber\\ && - \sqrt{a_i(t)} \equiv \sqrt{a_{o}(t)} Y - \sqrt{a_{i}(t)}. \end{eqnarray} \tag{ 5 }$$

We rewrite Equations (2)–(5) as

$$\begin{eqnarray} && \sqrt{1 - e_{i}(t)^2} = 1 + \frac{1}{q} \bigg[ \left(Y \left[ \frac{\sqrt{a_o(t)} - \sqrt{a_{o}}}{\sqrt{a_{i}}} \right] + 1 \right)^{-1}\nonumber\\ && \times \big(q \big(\sqrt{1 - e_{i}^2} - 1 \big) - \big(p-q \big) \big) + \big(p-q \big) \bigg], \end{eqnarray} \tag{ 6 }$$

which illustrates how a transiting planet's eccentricity is predicted to change as a function of the variation in the external perturber's semimajor axis for a given MMR. However, this expression does require knowledge of the initial semimajor axes and eccentricities of both planets.

Assuming that the outer planet is much more massive than the inner planet, the unseen outer planet will experience greater variations of its orbital elements over time than the transiting inner planet. Here, we consider only the semimajor axis ratio range of 1.3–3.2. For e_o = 0, the variation of e_o(t) typically does not exceed 0.01, but may reach 0.1 for systems near the 2:1 PC. The e_o(t) profiles in this regime include sinusoidal, lightly modulated, and heavily modulated curves with characteristic periods ranging from a few days to a few years. For e_o = 0.5e_H, typically 0.01 ⩽ e_o(t) ⩽ 0.1, with a greater proportion of sinusoidal-like curves. For e_o = 0.75e_H, typically e_o(t) ⩾ 0.1, with variations up to ≈0.6. Many e_o(t) profiles have little apparent structure and are periodic only on long (>10 yr) timescales. The outer planets with the smallest eccentricity variation, at e_o(t) ≈ 0.01, have corresponding TTV RMS amplitudes of just 10 s. RMS signal amplitudes of ∼10³ s may correspond to e_o variations of <0.10 or >0.40. For e_o = e_H, e_o(t) profiles vary chaotically and have typical amplitudes >10⁴ s. For a_o/a_i > 3.2, e_o(t) profiles are generally flat except for those corresponding to e_o ≈ e_H, in which they vary chaotically.

Transiting systems at or near PCs are perhaps of the greatest observational interest. Therefore, we trace the time evolution of the resonant angle in each system to determine if and over what timescale, and with what amplitude, this angle librates. We consider the strongest (lowest order) PCs for q = 1, 2, and 3, and fix Π_i = ϖ_o = 0° and Π_o = 60°. For each PC, we sample the resonant argument evolution for all simulations with the three values of a_o/a_i simulated which are closest to the value of (p/q)^2/3. Here, we describe the results qualitatively for this narrow region of resonant phase space.

All these systems are near a PC, and at eccentricities near e_H are all within the libration width of the corresponding commensurability. However, only a fraction of these systems actually exhibit a librating resonant angle, and hence are "in" MMR. The resonant angle may librate for different timespans: throughout an observational campaign or periodically during that campaign. The importance of a system actually in MMR versus just having a period ratio close to that commensurability is arguable, but a topic we now address.

First, we consider a commensurability that is saturated by resonant systems. We plot the systems sampled for the 2:1 PC in Figure 13. The upper, middle, and lower panels represent systems at the three different initial semimajor axis ratios (1.5812, 1.5865, and 1.5919), and the left and right panels represent the ranges 0>e_o > 0.5e_H and 0.5e_H > e_o > e_H, with e_o decreasing downward. Every row of each panel represents a different system; the left column represents the TTV signal, with corresponding amplitude in seconds to the left; the middle column represents the "apsidal angle" (≡ϖ_o − ϖ_i); and the right column represents the resonant angle (2λ_o − λ_i − ϖ_o) in red and another angle (2λ_o − λ_i − ϖ_i) in green. The middle and right columns have 0°–180° plot ranges, whereas the left column's range is scaled according to S(Q). The curves are shown for a time evolution of about 1.1 yr (corresponding to N ≈ 100).

Zoom In Zoom Out Reset image size

The red curve librates for all 180 systems in the figure, demonstrating that all these planets are locked in MMR, at least for about 1 yr. At the lowest values of e_o, this libration occurs about 0°, and is hence said to be "symmetric" and "aligned." For other values of e_o, the libration is about 45° and is hence said to be "antisymmetric." The green and blue curves librate for just a few systems at each a_o value. The libration of these angles means that the system is in apsidal corotation resonance (ACR; Beaugé & Michtchenko 2003; Ferraz-Mello et al. 2003; Lee 2004; Kley et al. 2005; Beaugé et al. 2006; Michtchenko et al. 2006a; Voyatzis & Hadjidemetriou 2006; Sándor et al. 2007; Michtchenko et al. 2008a, 2008b). The location of these systems occurs at ≈0.2e_H, but does vary sensitively on the value of a_o. The 2:1 MMR is the only resonance studied which features a librating apsidal angle.

The primary difference between the upper, middle, and lower panels is the amplitudes of the TTV signals, which can vary by several factors despite the small (<0.006 AU) difference between a_o values. The highest signals (with >100 s) appear in two eccentricity groupings in the upper and middle panels, and just one eccentricity grouping in the lower panels. The first group of high signals in the upper and middle panels ends when the libration of the resonant angle becomes antisymmetric and the TTV signal qualitatively changes character. The gap between groups of high signals for the upper panels is five times smaller than that for the middle panels. For 0.5 < e_o < e_H, the libration curves, TTV signal profiles, and TTV amplitudes are relatively insensitive (varying by a few percent) to changes in a_o, unlike for 0 < e_o < 0.5e_H. In this low-eccentricity regime, the highest amplitudes seen in the middle panels occur at ACR; the same is not true for the lower or upper panels. In none of those panels does the TTV signal profile noticeably change at ACR.

Unlike for the 2:1 commensurability, other p:q commensurabilities feature systems that are not in resonance, a situation corresponding to a circulation of the red angle. For all q = 1, and p = 2, ...11, only the 6 : 1 PC contains no MMRs. Systems also appear in MMRs for the 3 : 2, 5 : 2, 7 : 2, 9 : 2, 7 : 3, 8 : 3, 10 : 3, and 11 : 3 PCs, but none for the 5:3 PC. However, recall importantly that these results hold for a narrow region of phase space where Π_i = ϖ_o = 0° and Π_o = 60° and for just one set of masses.

As e_o increases, the transition between non-resonant and resonant systems produces a sudden change in TTV signal amplitude (by a factor of at least several) for all commensurabilities where this transition can be seen. Therefore, TTV signals could clearly distinguish resonant from non-resonant systems, even if both are within the libration width of a particular commensurability. Furthermore, the transition between symmetric and antisymmetric libration corresponds to a sudden (by a factor of at least two) change in the TTV signal in several cases around the 2:1, 3:1, 4:1, 5:1, 7:1, 8:1, and 10:1 PCs. The signal jump caused by this shift in the libration center is typically not as great as the jump caused by the transition to resonance from circulation. For 0.75e_H ⩽ e_o ⩽ e_H, only the 2:1, 5:1, 7:1, 3:2, and 9:2 PC systems are in MMR. At the other commensurabilities, the resonant angle varies in a chaotic fashion, a qualitatively different behavior from the smooth circulation (featuring a continuous curve) found at non-resonant lower eccentricity regimes.

In previous sections, we have focused our investigations on TTV signal amplitudes. Their importance stems from the ability of observers to distinguish a TTV curve from measurement uncertainty. As we have demonstrated, the TTV signal amplitude is highly sensitive to the orbital parameters of both planets, and is of limited utility when attempting to identify an unseen planet. In order to help break the degeneracy, one can employ the shape of the TTV curve. As seen in Figure 1, TTV curves may take on a variety of forms. By incorporating this form into a TTV analysis, one may be able to better pinpoint the mass and orbital parameters of the unseen planet.

A quantifiable way of obtaining TTV shape data is to consider the autocorrelation function $\mathcal {A}$ at a given (time) lag L:

$$\begin{equation} \mathcal {A}_L = \frac{ \sum _{k=1}^{N-L} x_k x_{k+L} }{ \sum _{k=1}^{N} x_{k}^2 }, \end{equation} \tag{ 7 }$$

where x_k represents the deviation of the kth transit time from the constant period model (to be distinguished from S(Q) at N = k, which represents the RMS amplitude from the first kth observations). The function $\mathcal {A}_L$ is bounded by [ − 1, 1] and can be computed solely from the transit observation data alone. In principle, $\mathcal {A}_L$ provides N constraints, in addition to S(Q), on the mass and orbital parameters of the unseen planet.

We find that the autocorrelation function may provide a useful summary of the TTV signal shape (see Figures 14 and 15). Figure 14 plots the autocorrelation as contours for nine different values of L (2, 3, 4, 5, 6, 20, 50, 100, and 200) for N = 313 (≈3.5 yr) and Π_o = Π_i = ϖ_o = 0°. One notes immediately that although some "flame-"like features from Figure 2 are apparent for low (L ⩽ 6) and high (L = 200) lags, the contours generally sculpt the phase space differently than does S(Q). This difference is most apparent in the low-eccentricity (e_o ⩽ 0.2) regimes, which demonstrate rich structure in autocorrelation space. For L = 100, sharp "flames" puncture this regime for initial semimajor axis ratios of up to ≈4. This low-eccentricity regime is transformed drastically as one increases L from L = 2 to L = 6. Note importantly that a purple flume (with $\mathcal {A}_L$ close to unity) is almost precisely centered on the L : 1 PC for L = 2, 3, 4, 5, and 6 (corresponding to a_o/a_i = 1.59, 2.08, 2.52, 2.92, 3.30, and 3.66). Such a correlation is a promising sign of the utility of $\mathcal {A}_L$ in helping to diagnose plausible integrations for future study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above contours fixed the initial orbital angles, which we have shown to sensitively affect the TTV signal amplitude. When we instead fix the lag, e.g., at L = 3, and vary Π_o as in Figure 7, then at the contour resolution of Figure 14, the change in initial angles appears to have no effect on the value of $\mathcal {A}_L$. However, at finer resolutions, there is a variation. For some systems, this variation in $\mathcal {A}_L$ strongly mimics the variation in S(Q) as Π_o is changed. Figure 15 displays three secular systems: the system (XVII from Table 1) in the upper panel with L = 3 shows remarkable agreement in the shapes and relative extrema of the solid, dotted, and dashed curves (corresponding to ϖ_o = 0°, ϖ_o = 90° and ϖ_o = 180°) to the signal from the bottom left panel of Figure 9. In contrast, the middle panels, which feature the secular system XVIII, do not share the same agreement with the corresponding signal (not shown). In between these two extremes, the bottom panels, featuring the secular system XIX, illustrate different levels of correlation to the signal variation (not shown): strong for the dashed curve, moderate for the dotted curve, and none for the solid curve. Therefore, the RMS amplitude and autocorrelation value might be strongly or weakly correlated depending on the particular orbital configuration. The right panels of Figure 15 demonstrate how the shape and magnitude of the $\mathcal {A}_L$ curves can qualitatively change for two different (arbitrarily chosen) values of L.

One may consider how $\mathcal {A}_L$ varies with planetary mass as well as initial orbital angles. In many cases, we do find a correlation between the variation of $\mathcal {A}_L$ and S(Q) with M_o, even if the agreement fails to come close to that from some orbital angle regimes. The variations induced by changing M_o for a given L tend to be greater by an order of magnitude than those induced from orbital angles alone. In a sense, the variation caused by the angles is a modulation to those caused by the masses.

These preliminary results suggest a promising avenue in which to pursue the inverse TTV problem, particularly if additional effects discussed in Section 8 are coupled with amplitude and shape data. Future work will entail identifying the critical lags which produce the autocorrelation extrema, quantifying how low (L < 10) lags track L:1 PCs, and incorporating S(Q) values and $\mathcal {A}_L$ for N − 1 values of L all as independent parameters to help remove the degeneracy of the inverse TTV problem.

8.1. Lagrange Unstable Systems

The highest signature TTV amplitudes (S(Q) ≳ 1000 s), like those generated from the curve in the fourth panel of Figure 1, might be the result of erratic changes in the orbital parameters of the planets. First, we confirmed that the highest eccentricity systems simulated for this work are provably Hill stable. In principle, such systems might be manifestly unstable even though they satisfy Hill stability. Although the planetary orbits will never cross in a Hill stable system, the outer planet might get ejected because the outer planet is not bounded, or is Lagrange unstable. If this is the case, then the timescale for this ejection is expected to be orders of magnitude less than the lifetime of the system, so observing planets in such configurations would be rare.

We attempt to place an upper limit on the fraction of systems which exhibit this manifestly unstable behavior by investigating configurations which produce sudden time variations of orbital elements such as the semimajor axis, periastron and apastron, mean longitude, and longitude of pericenter. Such estimates may be used as a Lagrange stability filter in additional statistical studies of TTVs. We find that comparing the maximum semimajor axis difference or the ratio of both planets over the length of the simulation with their initial values provides a useful diagnostic. In particular, we define

$$\begin{equation} \chi \equiv \frac{{\rm max}\left(a_o(t) - a_i(t)\right)}{a_o - a_i}. \end{equation} \tag{ 8 }$$

Therefore we compute χ for every one of the 120,000 simulations studied in our fiducial configuration of Figure 2, plus in selected configurations from Figure 11, representing extreme values of the transiting and exterior planet masses. We report the percent of all systems with χ less than given values from 1.1 to 10.0 in the left panel of Figure 16. The figure demonstrates that the large majority of Hill stable systems have planets which do not vary their semimajor axis difference by more than 10%; in these relatively quiescent systems, the outer planet is unlikely to be ejected. All but a few systems simulated with a transiting planet of 0.1 M_J exhibit this quiescent behavior. Alternatively, systems with a more massive transiting planet are the most likely to exhibit instability. Nearly 20% of all systems sampled with a 10 M_J transiting planet exhibit changes in their semimajor axis ratio approaching 1000%. Such systems could become unstable (in the sense of the outer planet being ejected) on timescales much shorter than the system lifetime.

Zoom In Zoom Out Reset image size

In order to help determine what Hill stable initial conditions most likely generate this suggested Lagrange instability, one can consider the initial eccentricity of the outer planet. We discover that highly eccentric outer planets close to the Hill stability limit produce the sudden changes in the semimajor axis characterized by χ ≳ 2.0. In order to quantify this finding, we plot in the right panel of Figure 16 the median (dots) and mean (squares) values of e_o as a percent of e_H which cause χ to be less than fixed values for the same mass configurations as in the left panel. Note that these eccentricities correspond to the highest TTV signals (S(Q) ≳ 1000 s; purple and blue regions in "flames" plots) in the fiducial case. We caution that these results are statistical in nature. The evolution of a particular system is highly dependent on the initial orbital element configuration and will likely require long term (≫10 yr) study in order to determine whether they will individually tend toward instability. We find this sensitive dependency by noting that for e_o = e_H and for given (initial) sets of a_o/a_i, altering initial sets of orbital angles can cause χ to vary from ≈1.1 to ≈10.0.

8.2. Radial Velocity Follow-up

Ideally, the existence of unseen planetary companions suggested from TTVs can be confirmed with radial velocity observations. However, the low mass and orbital separation of such unseen perturbers may preclude them from radial velocity detection (as well as observational limitations due to stellar properties). In this section, we estimate the maximum magnitude of the radial velocity semiamplitude produced by an external perturber in a guaranteed-to-be-stable multi-planet exosystem with a hot Jupiter. Because radial velocity semiamplitude $\equiv K_o \propto 1/\sqrt{1 - e_{o}^2}$, a planet could most easily be detected by radial velocity surveys when its eccentricity is highest. At the Hill stability limit (Gladman 1993),

$$\begin{eqnarray*} & & \left(\mu _o + \mu _o \mu _i + \mu _i \right)^{-3} \left(\mu _o + \frac{\mu _i}{\alpha } \right) \big(\mu _i \sqrt{\alpha } + \mu _o \sqrt{1 - e_{o}^2} \big)^2 \nonumber \end{eqnarray*}$$

$$\begin{eqnarray} & & =1 + \frac{3^{4/3} \mu _o \mu _i}{\left(\mu _o + \mu _i \right)^{4/3}} - \frac{\mu _o \mu _i \left(11 \mu _i + 7 \mu _o \right)}{ 3 \left(\mu _o + \mu _i \right)^2 }, \end{eqnarray} \tag{ 9 }$$

where α ≡ a_i/a_o, μ_i ≡ M_i/M_⋆, and μ_o ≡ M_o/M_⋆. Furthermore, the radial velocity semiamplitude of the outer planet in the coplanar edge-on case is (Lee & Peale 2003)

$$\begin{equation} K_{o,\rm max} = \left(\frac{2 \pi G}{P_o} \right)^{1/3} \frac{M_o}{\left(M_{\star } + M_o + M_i \right)^{2/3} } \frac{1}{\sqrt{1 - e_{o}^2}}, \end{equation} \tag{ 10 }$$

$$\begin{equation} P_o = \frac{2 \pi a_{o}^{3/2}}{G \left(M_{\star } + M_i + M_o \right)}. \end{equation} \tag{ 11 }$$

Therefore, the maximum semiamplitude of the external perturber, assuming M_⋆ ≫ M_i, M_o, is

$$\begin{equation} K_o \approx \frac{\sqrt{\frac{G M_{\star }}{a_o}} \mu _{o}^2}{ \left(\mu _i + \mu _o \right)^{3/2} \left(\frac{\mu _i}{\alpha } + \mu _o \right)^{-1/2} - \mu _i \sqrt{\alpha } }. \end{equation} \tag{ 12 }$$

This limiting equation, plotted in Figure 17, assumes that (1) the system is seen edge-on and (2) e_o = e_H. The figure, which displays curves for given masses of 1 M_⊕, 5 M_⊕, 10 M_⊕, 20 M_⊕, 50 M_⊕, and 1 M_Saturn (in ascending order of the curves on the plot) demonstrates that given a hot Jupiter at 0.05 AU around a solar-like star, K_o is a stronger function of M_o than of α. A single Earth-mass external perturber is not yet detectable by radial velocity measurements, but a super-Earth may be detectable by this means, independent of planetary separation. Once detected by radial velocity measurements, those observations can be combined with TTVs in order to help characterize the architecture of the system (Meschiari & Laughlin 2010).

Zoom In Zoom Out Reset image size

8.3. Light Travel Time

The variation in LTT during each transit may exceed the previously cited detectability threshold of 10 s depending on the planetary masses and orbital configuration. Therefore, TTV models might need to incorporate leading-order LTT effects in order to correctly describe the motion. All our simulations incorporate this effect. Here, we show where the effect might be important. In order to help determine when LTT effects may be neglected, we plot this contribution in Figure 18. Note the domain of the plot is three times as wide as those from Figures 2 and 3, but still sampled with 400 points along the X-axis. Figure 18 quantifies the contribution of LTT effects for N = 874 (left panels) and N = 50 (right panel). The contours plot the quantity log(S(Q)/T(Q)), where T(Q) is the TTV RMS amplitude produced by LTT effects alone. The figure demonstrates that the LTT's contribution for M_o = M_⊕ for an idealized 10 yr of coverage can be at the few percent to tens of percent level, but does not dominate the TTV signal. However, for a more massive external perturber, LTT can dominate the signal, especially for hierarchical (widely separated) systems with an outer planet eccentricity up to 0.5e_H.

Zoom In Zoom Out Reset image size

8.4. Extensions

Our goal was to illustrate and quantify the challenges involved in characterizing the mass and orbital parameters of an external perturber in a system with an observed transiting hot Jupiter through the use of TTVs. TTV signals may vary by orders of magnitude due to ∼10⁻³ AU shifts in the semimajor axes, ∼0.005 shifts in eccentricity, or ∼1° shifts in orbital angles. However, TTVs are sensitive probes of hidden planets, and can suggest the existence of external perturbers due to signatures produced at high-order period commensurabilities (PCs; p:q, where p ≲ 20 and q ≲ 3). More specific conclusions from this work are: (1) moderate (tens of seconds) values of signal amplitudes could indicate that the outer planet eccentricity could be as high as the Hill eccentricity, (2) high amplitude (>10³ s) signals do not necessarily imply close proximity to a PC, (3) high amplitude (>10³ s) signals are always indicative of a high value (at least 4/5th of the Hill eccentricity) of the outer planet eccentricity, (4) near-PC amplitudes are often higher than in-PC amplitudes, a result largely independent of the number of observed transits, N, (5) signal amplitude is a non-monotonic function of N, (6) increasing the external mass generally increases the TTV signal, (7) increasing the internal (hot Jupiter) mass generally decreases the number of easily discernible PCs, and (8) the distinctive signal amplitudes for systems at PCs do not necessarily imply that those systems are actually captured in an MMR. We propose using the shape data of a TTV curve through the autocorrelation function as a method to help characterize the external planets, and believe that the method could provide a promising avenue of future study.

We thank the referee for insightful observations and Eric Agol, Dan Fabrycky, Nader Haghighipour, Matt Holman, David Nesvorný, and Jason Steffen for valuable and extensive discussions. This material is based upon work supported by the National Science Foundation under grant No. 0707203.

Here we derive the estimates for the libration widths drawn in Figure 2 for selected PCs of up to eighth order and discuss the goodness of the approximation. We take "libration width" to mean roughly the a_o range for a given e_o in which MMR locking between two planets is likely. Because of the hierarchical nature of the three masses for the systems we consider, and the initially circular orbit of the more massive planet, we assume that the one argument given by Equation (1) of the disturbing function will dominate the gravitational potential. This approximation is vital to obtaining a tractable analytic expression.

We can assess the goodness of the approximation of neglecting terms that include e_i(t) by using Veras's (2007) term-based integrator. We tested systems locked in the 3:1 and 5:1 MMRs by including different numbers of terms up to fourth order. For all systems, we set M_⋆ = M_☉, M_i = M_J, M_o = M_⊕, a_i = 0.05, and a_o = a_i(p/q)^2/3. For the 3:1 MMR, we set e_o = 0.30, λ_i = 100°, λ_o = 0°, ϖ_i = 340°, and ϖ_o = 180° in order to produce resonant behavior. For the 5:1 MMR, we set e_o = 0.45, λ_i = 5°, λ_o = 170°, ϖ_i = 3°, and ϖ_o = 19° in order to produce resonant behavior. Importantly, for the 3:1 MMR, we set e_i = 0.00, but for the 5:1 MMR, we set e_i = 0.01. Initializing the hot Jupiter with a nonzero eccentricity might cause additional disturbing function terms to become significant, and is realistic given that an external terrestrial-mass perturber can typically force the hot Jupiter's eccentricity to values of ≈0.005.

We test these systems by including all arguments that, in their coefficients, include e_i(t) (Case 1), and by neglecting all such arguments (Case 2). For both cases and for both MMRs, the resonant angle librates, and does so about 180°. In Case 1, for the 3:1 and 5:1 MMRs, the libration amplitudes are ≈55° and ≈65°; Case 2 changes these values by at most a few degrees. The range of e_o(t) over 3 yr in Cases 1 and 2 is 0.3040–0.2639 = 0.0401 and 0.3035–0.2586 = 0.0449 for the 3:1 MMR, and 0.4669–0.4346 = 0.0323 and 0.4669–0.4346 = 0.0323 for the 5:1 MMR. These results demonstrate that neglecting these additional arguments affects the motion in too small a manner to necessitate inclusion in this general, qualitative study.

With a single disturbing function argument, ϕ, we can approximate a libration width according to the same prescription found in Murray & Dermott (2000). The angle behaves like a pendulum so that

$$\begin{equation} \ddot{\phi } = \omega ^2 \sin {\phi }. \end{equation} \tag{ A1 }$$

One can find the value of ω through a combination of (1) the definition of ϕ, (2) use of the disturbing functions, R = R(ϕ), and (3) use of Lagrange's planetary equations (Murray & Dermott 2000). We can define two disturbing functions for the systems we study here as

$$\begin{eqnarray} & & R_1 = \frac{G \!M_i}{a_i} e_{o}^{p-q} \left(\alpha f_d + f_i \right) \cos {\phi } \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} & & R_2 = \frac{G\! M_o}{a_o} e_{o}^{p-q} \left(f_d + f_e \right) \cos {\phi }, \end{eqnarray} \tag{ A3 }$$

where G is the gravitational constant, α = a_i/a_o, and f_d, f_i, and f_e are the functions of α and e_o. By defining σ_k = G(M_⋆ + M_k), for k = i, o, we find

$$\begin{eqnarray} & & \ddot{\phi } \approx p \dot{n_{o}} - q \dot{n_{i}} \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} & & = -\frac{3}{2} p \sigma _{o}^{1/2} a_{o}^{-5/2} \dot{a}_o - \frac{3}{2} q \sigma _{i}^{1/2} a_{i}^{-5/2} \dot{a}_i \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} & & = \frac{-3 p}{a_{o}^2} \frac{\partial R_1}{\partial \lambda _o} - \frac{3 q}{a_{i}^2} \frac{\partial R_2}{\partial \lambda _i} \end{eqnarray} \tag{ A6 }$$

so that

$$\begin{equation} \omega ^2 = \frac{3 G e_{o}^{p-q}}{a_o a_i} \left(\frac{p^2 M_i}{a_o} \left[ \alpha f_d + f_{id} \right] + \frac{q^2 M_o}{a_i} \left[ f_d + f_e \right] \right) \end{equation} \tag{ A7 }$$

The energy of a pendulum is $E = (1/2)\dot{\phi }^2 + 2 \omega ^2 \sin ^2(\phi /2)$. Hence, the maximum energy is 2ω², and equating this value with E yields a relation between $\dot{\phi }$ and ϕ. One can combine this relation with the Lagrange planetary equation for $\dot{a}$ in order to obtain

$$\begin{equation} da_o = \pm \frac{p G a_{o}^{1/2} M_i}{\omega a_i \sigma _{o}^{1/2}} e_{o}^{p-q} \left[ \frac{a_o}{a_i} f_d + f_{id} \right] \frac{\sin {\phi }}{\cos {(\phi /2)}} d\phi \end{equation} \tag{ A8 }$$

which can be integrated to finally obtain the (maximum) libration width δ_a:

$$\begin{equation} \delta _a = \frac{ \frac{2}{\sqrt{3}} \frac{p a_o M_i}{a_{i}^{1/2} \left(M_{\star } + M_o \right)^{1/2} } e_{o}^{(p-q)/2} \left(\alpha f_d + f_{id} \right) }{\sqrt{ \frac{p^2 M_i}{a_{o}} \left(\alpha f_d + f_{id} \right) + \frac{q^2 M_o}{a_{i}} \left(f_d + f_e \right) } }. \end{equation} \tag{ A9 }$$

In the approximation M_o = 0,

$$\begin{equation} \delta _a = a_o e_{o}^{(p-q)/2} \sqrt{\left(\frac{4}{3 \alpha } \right) \left(\frac{M_i}{M_{\star }} \right) \left(\alpha f_d + f_{id} \right) }. \end{equation} \tag{ A10 }$$

For most MMRs, f_id = f_e = 0. However, when q = 1, one of the most important classes of MMRs for TTVs, both f_id and f_e are nonzero. Formulas for f_d, f_id, and f_e are provided in Murray & Dermott (2000) up to fourth order in e_o. S. F. Dermott (2009, private communication) has provided us with terms up to eighth order. Explicit formulas for f_d in terms of α may be found in Veras & Armitage (2007).

The eccentricity at which Equation (A10) holds is restricted by the Sundman criterion (Ferraz-Mello 1994; Sidlichovsky & Nesvorny 1994), a fundamental convergence criterion on the planar expansion of the disturbing function from Ellis & Murray (2000). This criterion can be expressed as

$$\begin{equation} a_i D(e_i) < a_o d(e_o), \end{equation} \tag{ A11 }$$

where

$$\begin{eqnarray} D(y) &=& \sqrt{1 + y^2} \cosh {z} + y + \sinh {z}, \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} d(y) &=& \sqrt{1 + y^2} \cosh {z} - y - \sinh {z} \end{eqnarray} \tag{ A13 }$$

such that z is implicitly defined as z = qcosh z. This restriction prevents computation of the libration width for all eccentricities up to the Hill stability limit. Figure 5 of Nesvorný & Morbidelli (2008) shows a detailed view of the difference between the Hill stability limit and the Sundman convergence limit for several different mass ratios.