{'fraction_non_alphanumeric': 0.05020021654321412, 'fraction_numerical': 0.026517950745011773, 'mean_word_length': 3.767944936086529, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics and Transit Variations of Resonant Exoplanets

The Transit Timing Variations (TTVs) are deviations of the measured mid-transit times from the exact periodicity. One of the most interesting causes of TTVs is the gravitational interaction between planets. Here we consider a case of two planets in a mean motion resonance (orbital periods in a ratio of small integers). This case is important because the resonant interaction can amplify the TTV effect and allow planets to be detected more easily. We develop an analytic model of the resonant dynamics valid for small orbital eccentricities and use it to derive the principal TTV terms. We find that a resonant system should show TTV terms with two basic periods (and their harmonics). The resonant TTV period is proportional (m/M_*)^(-2/3), where m and M_* are the planetary and stellar masses. For m=10^(-4) M_*, for example, the TTV period exceeds the orbital period by ~2 orders of magnitude. The amplitude of the resonant TTV terms scales linearly with the libration amplitude. The ratio of the TTV amplitudes of two resonant planets is inversely proportional to the ratio of their masses. These and other relationships discussed in the main text can be used to aid the interpretation of TTV observations.


Introduction
Photometric observation of transits is one of the most powerful methods of planet detection. This method relies on a possibility that, if the planet's orbit is viewed nearly edge-on, the planet may repeatedly transit over the disk of its host star and periodically block a small fraction of the starlight. Thus, by monitoring the host star's brightness, the planet's presence can be revealed by a small dip in the photometric lightcurve. The main properties of the planet, such as its physical radius and orbital period, can be inferred from transit observations.
The spacing of transit lightcurves would be exactly the same over the course of observations if a planet moved on a strictly Keplerian orbit. Several dynamical effects, however, can produce deviations from the Keplerian case and induce the Transit Timing Variations (TTVs). TTVs were originally proposed as a non-transiting planet detection method (Miralda-Escudé 2002, Agol et al. 2005, Holman & Murray 2005, but have found more use in validating the transiting planet candidates from NASA's Kepler mission (e.g., Holman et al. 2010;Lissauer et al. 2011). Only a handful of non-transiting planets have been so far detected and characterized from TTVs (e.g., Nesvorný et al. 2012Nesvorný et al. , 2013. A significant progress has been made in the theoretical understanding of various dynamical causes of TTVs. These efforts were pioneered by Agol et al. (2005). Heyl & Gladman (2007) focused on a long-period interaction between planets and showed that the apsidal precession of their orbits can be detected only with a long TTV baseline. Nesvorný & Morbidelli (2008), on the other hand, developed a general analytic model for short-period TTVs and showed that they can be used, under ideal circumstances, to uniquely determine the mass and orbital parameters of the interacting planets.
Explicit analytic formulas for short-period TTVs are now available for zeroth- (Agol et al. 2005, Deck & Agol 2015 and first-order terms  in planetary eccentricities. The important case of near-resonant TTVs was highlighted in Lithwick et al. (2012) for the first-order resonances and in  and Hadden & Lithwick (2016) for the second-order resonances. The near-resonant TTV signal is a special case of the short-periodic variations when one harmonic becomes amplified due to a proximity of the system to a mean motion resonance.  considered a case of co-orbital planets and showed that co-orbital TTVs are expected to have a characteristic saw-tooth profile (for horseshoe orbits).
While the analytic works cited above cover a wide range of dynamically plausible planetary configurations, none of them (except for Agol et al. 2005) considered the case of a fully resonant planetary system, where two (or more) orbits are inside a mean motion resonance (e.g., 2:1, 3:2). This case is important because planetary migration in a protoplanetary gas disk should bring early planets into resonances (e.g., Masset & Snellgrove 2001). Indeed, many known planetary systems are consistent with having resonant orbits (e.g., Winn & Fabrycky 2015). This motivates us to consider the resonant case. Our main goal is to under-stand how resonant TTVs arise, and how their period and amplitude scale with planetary parameters.
The resonant case has not received much attention in the TTV literature so far at least in part because the resonant interaction of planets is complex and not easily amenable to analytic calculations. The short-periodic or near-resonant TTV signals, for example, can be computed with the standard methods of perturbation theory, where the unperturbed (Keplerian) motion is inserted in the right-hand side of the dynamical equations, and a linear variation of the orbital elements is obtained by quadrature. This method fails for the fully resonant orbits mainly because the resonant dynamics is non-linear.
Here, we take advantage of recent advances in the theoretical understanding of resonant dynamics (e.g., Batygin & Morbidelli 2013a; hereafter BM13) and derive approximate formulas for TTVs of a resonant pair of exoplanets. In Section 2, we first give a summary of resonant TTVs. The goal of this section is to highlight the main results of this work. Sections 3 and 4 explain how we obtained these results. We first show how a fully analytic solution can be obtained for a first-order resonance (Section 3). We then proceed to expand the exact solution in the Fourier series and explicitly derive the periods and amplitudes of the leading TTV terms (Section 4). The application of these results to real planetary systems is left for future work.

Summary of Resonant TTVs
Consider a planar system of two planets with masses m 1 and m 2 orbiting a central star with mass M * . The two planets gravitationally interact to produce TTVs. The TTVs signal of the two planets, δt 1 and δt 2 , can be approximated by where n j is the mean orbital frequency of planet j, h j = e j sin ̟ j , λ j is the mean longitude, e j is eccentricity, and ̟ j is the longitude of periapsis. Observer is assumed to see the system edge-on and measures angles λ j and ̟ j relative to the line of sight. Expression (1) is valid for small orbital eccentricities and small variations of the orbital elements, δλ j and δh j (Nesvorný 2009).
Here we consider a case with two planets in a first-order mean motion resonance (such that n 1 /n 2 ≃ k/(k − 1) with integer k), and proceed by calculating the variations of orbital elements due to the resonant interaction (Sections 3 and 4). The final expressions are given as the Fourier series with harmonics of two basic frequencies. As for δλ j we have Here, ν = (3/2)[(k−1) 2 n 1 /Λ 1 +k 2 n 2 /Λ 2 ], Λ j = m j GM * a j , a j is the semimajor axis, G is the gravitational constant, P τ is the period of resonant librations in scaled time units (P τ ∼ 2-4 in the most cases of interest; Section 3.7), A Ψ is the amplitude of the resonant oscillations of action Ψ (ranging from zero for an exact resonance to >1, where the approximation used to derive Eq. (2) starts to break down), f is the frequency of resonant librations, and 0 ≤ ǫ < 1 encapsulates the emergence of higher-order harmonics of f . 1 The resonant frequency f is given by where f 1 and f 2 are the resonant coefficients of the Laplacian expansion of the perturbing function (Table 1). The scaling of the resonant period, P = 2π/f , with different planetary parameters is discussed in Section 3.7. For practical reasons, Eq. (2) have been truncated at the first order in ǫ. Higher orders in ǫ and higher harmonics of the libration frequency can be computed using the methods described in Section 4.
The second term in Eq. (1) is related to the variation of eccentricities and apsidal longitudes of the two planets. In the most basic approximation (Section 4.3), it can be written as where P j = 2π/n j are the orbital periods, θ 0 is the initial value of θ = kλ 2 − (k − 1)λ 1 , f θ is the frequency of θ defined in Eq. (65). In the libration regime, f θ ≪ f . The resonant TTVs can therefore be understood as a sum of librational variations given in Eq.
(2) and slower modulation of the TTV signal given in Eq. (5). In Eq. (5), C u,0 is a coefficient of the order of unity (Eq. (67) in Section 4.3). For small libration amplitudes, A Ψ C u,0 ≃ √ 2Ψ 0 , where Ψ 0 is the initial value of the resonant action Ψ defined in Eqs. (21) and (24). This means that, unlike in Eq. (2), TTVs from the variation of eccentricities and apsidal longitudes do not vanish when A Ψ = 0. The dependence of the TTV amplitude on planetary masses arises from A, B and η 1 in Eq. (5), Two main approximations were adopted to derive Eqs. (2) and (5). In the first approximation, we retained the lowest-order eccentricity terms in the resonant interaction of planets. In the second approximation, we assumed that the libration amplitudes are not very large, and expanded the exact solution in the Fourier series. Both these approximations are tested in Section 5. Here we just illustrate the validity of the Fourier expansion for small libration amplitudes (Figures 1 and 2).
Let us briefly consider an application of our results to the TTV analysis. We assume that the photometric transits are detected for planets 1 and 2, and that the orbital periods inferred from the transit ephemeris are such that orbital period ratio P 2 /P 1 ≃ k/(k − 1) with small integer k (indicating near-resonant or resonant orbits). Furthermore, TTVs are assumed to be detected for both planets. To be specific, let us consider a realistic case with the available TTV data spanning several years of observations, which is not long enough to resolve the frequencies related to the apsidal precession of orbits.
The first step of the TTV analysis is to apply the Fourier analysis to the TTV data. This may reveal that the TTV signal contains a basic period. In principle, this period can be one of following two periods: (1) the super period defined as P s = (k/P 2 − (k − 1)/P 1 ) −1 (e.g, Lithwick et al. 2012), or (2) the resonant libration period P = 2π/f with f defined in Eq. (3). In the first case, the system is not in the libration regime of the resonance and TTVs can therefore be interpreted using the expressions appropriate for the near-resonant dynamics (Lithwick et al. 2012. The results described in this work apply in the second case. Section 3.7 explains how the libration period can be used to constrain planetary masses. Specifically, P ∝ (m/M * ) −2/3 (Agol et al. 2005, Holman et al. 2010, and therefore larger planetary masses imply shorter TTV periods. For the 2:1 resonance with P 1 = 10 days and m 1 ≃ m 2 = 10 −4 M * , for example, the libration period is P ≃ 4.5 yr (Section 3.7).
In the next step, it can be useful to check if the Fourier analysis of the TTV data provides evidence for harmonics of the basic period. If that's the case, this can indicate that the libration amplitude is relatively large. A comparison of the amplitudes of different harmonics can then be used to constrain the parameter ǫ in Eq. (2), which is related to the libration amplitude via equations reported in Appendix A.
We then proceed by comparing the TTV amplitudes of the two planets. From Eq.
(2) we have that δλ 1 /δλ 2 ≃ −[(k − 1)/k] 2/3 m 2 /m 1 . The TTV amplitude ratio therefore constrains the ratio of planetary masses (Agol et al. 2005). Figure 3 shows how the TTV amplitudes depend on planetary masses. For m 1 ≪ m 2 , we obtain from Eq. (2) that the TTV amplitudes of the inner and outer planets, A 1 and A 2 , are where we denoted α = a 1 /a 2 . This means that the TTV amplitude of the inner planet is independent of masses, while that of the outer planet constrains m 1 /m 2 . For P 1 = 10 day, P τ = 3 and A Ψ = 1, A 1 = 1.52/(k − 1) days for a k:(k − 1) resonance. A similar analysis can be performed for m 1 ≫ m 2 .
In addition to the dependence on masses, the TTV amplitudes (linearly) depend on the libration amplitude A Ψ . We therefore expect that some degeneracy should exist between the planetary mass and libration amplitude, with large TTV amplitudes being produced either by large masses or large libration amplitudes. This degeneracy can be broken if the libration frequency harmonics are detected, providing constraints on the libration amplitude, or if the measured TTV signal also contains short-period (chopping) terms (e.g., Nesvorný et al. 2013). A detailed analysis of this problem is left for future work.

Analytic Model of Resonant Dynamics
Here we discuss a Hamiltonian model of the resonant dynamics. Our approach closely follows the work of BM13. We take several shortcuts to simplify the reduction of the Hamiltonian to an integrable system. Then, in Section 3.6, we present an exact analytic solution. This solution is used in Section 4 to derive Eq. (2) and (5).

Hamiltonian Formulation of the Problem
The Poincaré canonical variables of two planets orbiting their host star are denoted by (r 0 , r 1 , r 2 ; p 0 , p 1 , p 2 ). The coordinate vector r 0 defines the host star's position with respect to the system's barycenter. Vectors r 1 and r 2 are the position vectors of the two planets relative to their host star. Momentum p 0 is the total linear momentum of the system (p 0 = 0 in the barycentric inertial frame). Momenta p 1 and p 2 are the linear momenta of the two planets in the barycentric inertial frame. The Poincaré variables are canonical, which can be demonstrated by calculating their Poisson brackets.
Using Poincaré variables the differential equations governing dynamics of the two planets can be conveniently written in a Hamiltonian form, where the total Hamiltonian is a sum of the Keplerian and perturbation parts, H = H K + H per , with and Here we denoted M j = m j + M * and the reduced masses µ j = m j M * /M j , where j = 1 and 2 stand for the inner and outer planet, respectively.
We assume that the planets are near or in a first-order mean motion resonance such that the ratio of their orbital periods is P 2 /P 1 ≃ k/(k − 1) for some integer k ≥ 2. In terms of the osculating orbital elements we have where a 1 and a 2 are the semimajor axes of planets, and where e 1 and e 2 are the orbital eccentricities, λ 1 and λ 2 are the mean longitudes, and ̟ 1 and ̟ 2 are the longitudes of pericenter. In Eq. (10) we only retained two most important resonant terms, and the lowest (first-order) eccentricity power. This expression is thus valid only near a specific resonance and for low orbital eccentricities of both planets. The planetary orbits are assumed to be in the same plane such that all inclination-dependent terms vanish.
The coefficients f 1 and f 2 are functions of the semimajor axis ratio α = a 1 /a 2 < 1, and can be written as where b (k) 1/2 (α) are the Laplace coefficients, D = d/dα, and δ k,2 is the Kronecker symbol (e.g., Brouwer & Clemence 1961). The last term in the expression for f 2 only appears if k = 2 (2:1 resonance). The Laplace coefficients were computed using the recurrences described in Brouwer & Clemence (1961, Secs. 15.7 and 15.8). Their values are reported in Table 1.

Resonant Variables
The resonant Hamiltonian discussed in the previous section needs to be written in canonical variables. The standard choice is the Delaunay elements Note that Γ j ≃ (1/2)Λ j e 2 j for small e j . The elements (Λ j ,Γ j ) are the canonical momenta, and (λ j , γ j ) are the conjugated canonical coordinates.
The resonant Hamiltonian H written in terms of the Delaunay elements has four degrees of freedom (DOF). Using canonical transformations we will reduce it to an integrable one-DOF system. We first perform a canonical transformation to the resonant canonical variables defined as When the new variables are inserted in Eqs. (9) and (10), it becomes clear that the resonant Hamiltonian depends on σ 1 and σ 2 , but not on λ 1 and λ 2 . Therefore, dK j /dt = −∂H/∂λ j = 0 and both momenta K 1 and K 2 are the new constants of motion. By summing them, we find that K 1 + K 2 = Λ 1 + Λ 2 − Γ 1 − Γ 2 is the total angular momentum of the system. Also, kK 1 + (k − 1)K 2 = kΛ 1 + (k − 1)Λ 2 = const. implies that any small changes of the semimajor axes, δa 1 and δa 2 , must be anti-correlated, and have relative amplitudes such that δa 1 /δa 2 ≃ −α 1/2 res (k − 1)m 2 /km 1 , where we denoted α res = [(k − 1)/k] 2/3 .

Approximation for Small Semimajor Axis Changes
The transformation to the resonant canonical variables produced a 2-DOF Hamiltonian H = H(σ 1 , σ 2 ; Γ 1 , Γ 2 ). In the next step, we assume that any changes of the semimajor axes of planets are small. Specifically, we write a j = a * j + δa j , where a * j is some reference value and δa j ≪ a * j , insert this expression into the Keplerian part of the Hamiltonian (Eq. 9), and expand it in powers of δa j . The first and second-order terms in δa j are retained. We then use δa j /a * j = 2δΛ j /Λ * j and rewrite all expressions in terms of Λ * j and δΛ j . 2 Finally, we substitute δΛ j → Λ j − Λ * j and drop all (dynamically unimportant) constant terms. This substitution is useful because it allows us to work with Λ 1 and Λ 2 instead of their variations δΛ 1 and δΛ 2 . Finally, we express H K in terms of the canonical variables defined in Eq. (13). This leads to with Here we denoted n j = GM j /a * 3 j and Λ * j = µ j GM j a * j . The quantities n j , n s and ν are constant parameters. Note that n s is related to the so-called super frequency, which is the expected frequency of the TTV signal for a pair of near-resonant planets (e.g., Lithwick et al. 2012). In addition to the usual term, kn 2 − (k − 1)n 1 , here n s also includes a correction that is a second order in the eccentricity (through its dependence on K 1 and K 2 ).
As for H per in Eq. (10), we have that e j = 2Γ j /Λ j for small eccentricity. In addition, because H per is already small, we do not need to retain terms proportional to δa j (thus Λ j → Λ * j ). The perturbation function then admits the following form where we introduced constants A = f 1 / Λ * 1 and B = f 2 / Λ * 2 .
2 A general result can also be obtained by directly performing the Taylor expansion in Λ j = Λ * j + δΛ j .
The momenta Φ 1 and Φ 2 defined by the transformations discussed above can be expressed in a compact form where z j = x j + ıy j = 2Γ j exp ıσ j . Condition Φ 2 = const. thus defines a circle in the plane of complex variables z 1 and z 2 and requires that Bz 1 − Az 2 lies on the circle at any time. Also, φ 1 = arg(Az 1 + Bz 2 ) and φ 2 = arg(Bz 1 − Az 2 ), where arg(z) denotes the argument of z.

Final Scaling
The Hamiltonian (20) depends on parameters n s , ν and C. We rescale Φ 1 and time to bring the Hamiltonian to a simple form where the parametric dependence is expressed by a single parameter The Hamiltonian equations are where τ relates to normal time t by with η 2 = νC 2 . When a solution of Eqs. (25) and (26) is found, the scaling parameters η 1 and η 2 can be used to map the solution back to the original variables.
The Hamiltonian (22) and the corresponding Eqs. (25) and (26) have been extensively studied in the past. They are equivalent to the second fundamental model of resonance (Henrard & Lemaître 1983), and to the Andoyer model discussed in Ferraz-Mello (2007). Here we first consider the dynamical flow arising from this Hamiltonian and its dependence on δ. In the next section, we show that Eqs. (25) and (26) have an exact analytic solution in terms of the Weierstrass elliptic functions.
The equilibrium points of Eqs. (25) and (26) control the general structure of the dynamical flow. Since dΨ/dτ = 0 implies that sin ψ = 0 in (26), the equilibrium points occur for ψ = 0 or π. The equilibrium values of Ψ are obtained from dψ/dτ = 0, which leads to a problem of finding roots of a cubic equation Ψ 3 − 2δΨ 2 + δ 2 Ψ − 1/8 = 0. There is only one (stable) equilibrium point for δ ≤ δ * = (27/32) 1/3 ≃ 0.945. This equilibrium point is located at ψ = π and 0 < Ψ < 1.26 ( Figure 4). The dynamical flow around the equilibrium point is simple. When projected to the (Ψ cos ψ, Ψ sin ψ) plane, 3 the trajectories are concentric, slightly deformed circles centered on the equilibrium point ( Figure 5a). Ψ changes only slightly during each cycle, and ψ circulates in the clockwise direction for most initial conditions except for the ones located very close to the equilibrium point, where ψ oscillates around π.
For δ = δ * , the stable equilibrium point is already substantially displaced from the origin. The dynamical transition is heralded by the appearance of a cusp trajectory (shown by a thin line in Figure 5b). Three equilibrium points appear for δ > δ * , two of which are stable and one is unstable. The unstable equilibrium point is traversed by a separatrix (shown by thick lines in Figures 5c and 5d), which surrounds the truly resonant trajectories for which ψ librates around π. The stable equilibrium point in the center of the resonant island is a smooth extension of the equilibrium point from δ < δ * to δ > δ * . It is located at ψ = π and Ψ = Ψ eq , where Ψ eq increases with δ ( Figure 4). For large δ, Ψ eq ≃ δ (Henrard & Lemaître 1983).

Exact Analytic Solution
A rearrangement of Eqs. (25) and (26) shows that the momentum Ψ satisfies with the right-hand side being a quartic polynomial The coefficients are where H (Ψ 0 , ψ 0 ) = h 0 is the conserved energy, and (Ψ 0 , ψ 0 ) is the initial condition. Equation (28) admits a general analytic solution (see Whittaker & Watson 1920, Sec. 20.4) is the n-th derivative of f at Ψ = Ψ 0 (n = 1, . . . , 4), and S = sgn[sin ψ 0 ] is the sign function (1 for sin ψ 0 > 0 and −1 for sin ψ 0 < 0). The time dependence of the solution is given in terms of the Weierstrass elliptic function ℘(z), and ℘ ′ (z) = d℘(z)/dz, whose invariants are Interestingly, as far as we know, the general solution (31) has not been discussed in the literature. Ferraz-Mello (2007) mentioned a particular solution valid for sin ψ 0 = 0, which implies that f 0 = 0 (see also Shinkin 1995).

Resonant Period
The solution (31) is periodic with a period for ∆ > 0, and for ∆ < 0. Here, K(k) is the complete elliptic integral of the first kind. Figure 7 shows the period of small amplitudes librations around the stable equilibrium points. Away from the resonance and if the eccentricities are small, n s /(2ν) ≫ Φ 2 , and the first bracketed term in Eq. (23) outweighs the second. In this case, δ is related to the super frequency, and increases, in the absolute value, when the system moves away from the resonance. In this sense, δ is a measure of the distance from the resonance. 4 Negative values of δ imply that the planetary orbits are spaced more widely than the actual resonance (P 2 /P 1 > k/(k − 1)), while δ > 0 means that the orbits are packed more tightly (P 2 /P 1 < k/(k − 1)). Also, δ ≃ η −1/3 2 n s /2 away from the resonance, and therefore the scaled period P τ ≃ π/|δ|. This represents a very good approximation of the period if δ < −2 or δ > 3 (see Figure 7).
The period P τ for different libration amplitudes is plotted in Fig. 8. The period is the shortest near the equilibrium point and increases with the libration amplitude. It becomes infinite on the separatrices. This corresponds to a situation when the amplitude of angle ψ becomes full π, and it takes infinitely long to reach the unstable equilibrium point at ψ = 0. Except for trajectories near the separatrix, the period values inside the libration island are P τ ≃ 2.5-4. They become shorter for δ > 5 (i.e., for higher eccentricities).
When rescaled according to Eq. (27), P t = η −1/3 2 P τ is the period of resonant librations in the normal time units (e.g., Julian days). It has the same dependence of P τ on δ and amplitude that we discussed above, and, in addition, contains an explicit dependence on the orbital period and planetary masses through the scaling parameter η 2 = νC 2 . If m 1 ≪ m 2 , then If, on the other hand, m 1 ≫ m 2 , then The multiplication coefficients in the square brackets are ≃1 for k = 2 (2:1 resonance) and decrease with increasing k. This expresses a stronger interaction of orbits that are more tightly packed for larger values of k. According to Eqs. (38) and (39), the orbital periods P 1 and P 2 set the basic time unit for P t .
The libration period scales with the mass of the more massive planet as (m/M * ) −2/3 , and is therefore shorter for a larger mass. It is insensitive to the mass of the lighter planet. If the masses of the two planets are comparable, then P t will depend on their combination via the scaling parameter η 2 . Figure 9 shows the contour plots of P t as a function of m 1 /M * and m 2 /M * . In the α → 1 limit, the period is proportional to [(m 1 + m 2 )/M * ] −2/3 . A reasonable approximation of the period is then These considerations have important implications for the TTV period and its scaling with various parameters. For example, for two planets in the 2:1 resonance with m 2 /M * = 10 −4 and m 1 ≪ m 2 , the expected TTV period is P t ∼ 78P 1 P τ . Thus, if P 1 = 10 days, and assuming that P τ ∼ 3 in the libration regime, the TTV period will be P t ∼ 6.4 years. We can therefore very roughly write in this case P t (2 : 1) ∼ 6.4 yr × P 1 10 d If, on the other hand, m 1 /M * = 10 −4 and m 1 ≫ m 2 , then the period will be approximately twice as long, because of having P 2 in Eq. (39) (instead of P 1 in Eq. (38)). Finally, if m 1 /M * = m 2 /M * = 10 −4 , then P t (2 : 1) ∼ 4.5 yr. The periods in the k > 2 resonances are expected to be shorter (Figure 9b).
This Hamiltonian has a simple solution (e.g., Batygin & Morbidelli 2013b). Defining complex variables z j = x j + ıy j , the solution can be written as where C j are integration constants related to the initial conditions. Now, since σ j = θ − ̟ j , we can define L12's 'complex eccentricities' w j = e j exp(ı̟ j ) and recast the approximate solution (44) as whereC j is the complex conjugate of C j . Here we assumed that θ = θ 0 + n s t and replaced Λ j → Λ * j . These assumptions are equivalent to those of L12 where the unperturbed (Keplerian) solution was inserted into the right-hand sides of the Lagrange equations, and the linearized solution was found by quadrature.
Eq. (45) can be compared to Eq. (A15) in L12. The first term in Eq. (45) is constant and corresponds to the 'free' term in Eq. (A15). The second term in Eq. (45) is identical to the second term in Eq. (A15) (this can be trivially shown by resolving notation differences). L12 proceeded by using the constants K 1 and K 2 (Eq. 13) to derive expressions for a j (t), which were then used to obtain λ j = n j t with n j = GM * /a 3 j (t). 5 Finally, L12 used a formula equivalent to Eq. (1) to compute TTVs. None of these steps requires a special clarification.
The main assumption of L12 was therefore to (effectively) neglect the term ν(Γ 1 +Γ 2 ) 2 in Eq. (42). For this to be valid |n s | ≫ ν(Γ 1 + Γ 2 ). Assuming that the masses and eccentricities of the two planets are comparable, and neglecting all factors of the order of unity, this condition can be written as e 2 j ≪ |P 2 /P 1 − k/(k − 1)|. This shows that the eccentricities cannot exceed certain threshold for the L12 formulas to be valid. This threshold is not 5 The method of L12 ignores the contribution to λ j from the derivatives of the Laplace coefficients. This is a correct assumption for the near-resonant orbits, because these terms have n s in the denominator, while the terms from n j (t) have n 2 s in the denominator. The latter terms are thus amplified near a resonance where n s is small. excessively restrictive. For example, for planetary orbits just outside the 2:1 resonance with P 2 /P 1 − 2 = 0.05, the eccentricities need to be e j ≪ √ 0.05 = 0.22.
In addition to the above eccentricity condition, the planetary orbits cannot be too close to a resonance where the non-linear effects become important even for negligible eccentricities. We tested this condition assuming initial orbits with small eccentricities and found that the L12 model is perfectly valid for δ < −4. A small, ∼10% discrepancy in both the amplitude and frequency of w j appears for δ ≃ −2. Further increase of δ leads to a situation where Eq. (45) is no longer an adequate representation of the resonant dynamics. According to our tests with e 1 ≃ e 2 0.01, the L12 model cannot usually be trusted for δ > −1. Figure 10 shows these thresholds as a function of the planetary mass and orbit period ratio.

Transit Timing Variations
If the variation of orbital elements is small, TTVs of two planets, δt 1 and δt 2 , can be computed from − n j δt j = δλ j + 2(δk j sin λ (see, for example, Nesvorný & Morbidelli 2008). Here, k j = e j cos ̟ j , h j = e j sin ̟ j , and λ j = λ * j + δλ j , where λ * j = n j (t − t 0 ) with constant n j . If the reference frame is chosen such that the orbital angles are measured with respect to the line of sight, transits occur when λ * j ≃ 0 (assuming small eccentricities). We therefore have −n j δt j = δλ j − 2δh j + O(e). The first-order eccentricity terms, which can be used to improve the validity for higher eccentricities, were given in Nesvorný (2009). TTVs will thus have a contribution from the mean longitude variation, δλ j , and another contribution from the variation of eccentricity and apsidal longitude, δh j = δ(e j sin ̟ j ). Defining Y j = − 2Γ j sin ̟ j , we have for small eccentricities that Γ j ≃ 1 2 Λ j e 2 j and h j = −Y j / Λ j . We then write Λ j = Λ * j + δΛ j with constant Λ * j , and obtain where we retained only the first-order terms in small variations.
It remains to compute δλ j and δY j . As for δλ j , we have dλ j /dt = ∂H K /∂K j where H K is given in Eq. (14) (note that K j appears in the n s term in (15)). Substituting Γ 1 + Γ 2 → Φ 1 + Φ 2 in Eq. (14), taking the derivative with respect to K j , and integrating with respect to t, we obtain Using the scaling relationship from Eqs. (24) and (27), we have that where Ψ(τ ) is given in Eq. (31). The first terms in (48) describe a uniform circulation of angles λ 1 and λ 2 (note that 4n j − 3K j (n j /Λ * j ) = n j for e j = 0, as expected). They will not contribute to TTVs. Instead, TTVs will arise from the integral term. As the sign in front of the integral term is positive for λ 1 and negative for λ 2 , TTVs of the two planets are anti-correlated. The amplitudes of δt 1 and δt 2 , denoted here by A λ,1 and A λ,2 , satisfy Thus, apart from a coefficient of the order of unity, the ratio of the TTV amplitudes from the λ terms is expected to be equal to the inverse of the planetary mass ratio.
As for δY j , we have where u 1 and v 1 were defined in Eq. (18). They are related to the Φ 1 and φ 1 variables via Eq. (19). Expressions (51) are derived in Appendix B. As we show in Appendix B, V 2 is a constant of motion, does not contribute to variations, and does not need to be computed.
Since the coefficients f 1 and f 2 have opposite signs (Table 1), A and B in Eq. (51) will also have opposite signs, and TTVs of the two planets will be anti-correlated. TTVs arising from terms in (51) will have amplitudes, A h,1 and A h,2 , such that The ratio of TTV amplitudes from the h terms is inversely proportional to the ratio of planetary masses. This is the same mass dependence as in Eq. (50). Thus, while resonant TTV period discussed in Section 3.7 can be used to constrain m 1 /M * and/or m 2 /M * , the TTV amplitude ratio is sensitive to m 2 /m 1 .
Above we reduced the problem in hand to the evaluation of v 1 cos θ − u 1 sin θ. The θ terms are simple. From Eq. (48) we have that As for the u 1 and v 1 terms, we obtain This is all we need for the computation of TTVs. To summarize, Ψ(τ ) from Eq. (31) needs to be inserted in (49) and (54). The integral in (49) needs to be computed and substituted into Eqs. (48) and (53). Then, θ, u 1 and v 1 obtained from Eqs. (53) and (54) are substituted into (51). The constant terms in λ j and Y j can be neglected, because they do not contribute to TTVs. Finally, Eq. (47) is used to compute TTVs.
In principle, the method outlined above can be used to derive fully analytic expressions for δt 1 and δt 2 , which would have the same general validity for low eccentricities as the analytic solution (31). For that, however, we would need to compute the derivative (needed for v 1 ) and integral (needed for λ j and θ) of Ψ(τ ) from (31). The integral ends up producing complex expressions (Appendix C). Here we therefore opt for a different approach, where we seek to find an expression for resonant TTVs in terms of the Fourier series.
For librations around the equilibrium point, we have that a ≫ b (Figure 11). We therefore identify a small parameter ǫ = b/a ≪ 1 and expand Eq. (55) in the Taylor series in ǫ.
Retaining only the first-order terms in ǫ (approximation of low-amplitude librations), we obtain Including higher order terms in ǫ would improve the validity of the approximation for large libration amplitudes. Next, we express cn and cn 2 in the Fourier series and Here we have k ′ = √ 1 − k 2 and q = exp(−πK ′ /K) with K ′ = K(k ′ ). E(k) is the complete elliptic integral of the second kind. General expressions for cn m (u, k), which can become useful when higher order terms in (56) are accounted for, can be found in Kiper (1984).
Both these Fourier series converge very rapidly, and we can therefore afford to use the lowest harmonics. In practice, given that the evaluation of u 1 from Eq. (54) will require a multiplication of the Fourier series, here we consider only the first and second harmonics of u = 2 √ γ(τ − τ 0 ). After substituting these terms into Eq. (56), we obtain Here we denoted the frequency f τ = 2π/P τ = π √ γ/K and D = C 1 /4a. To simplify things, we set τ 0 = 0 and drop all multiplication terms appearing from Eqs. (57) and (58)  Note that D is a proxy for the libration amplitude of Ψ ( Figure 11). Eq. (59) can be used to trivially compute the integral Ψdτ appearing in Eq. (49), which is then inserted into Eqs. (48) and (53). In the following text, the two TTV contributions, δλ j and δY j , will be considered separately.

Contribution from Mean Longitude Variations
The calculation of δλ j is simple. Retaining the periodic terms in Eq. (48) we obtain with the coefficients C λ,1 = 1 + ǫ and C λ,2 = −ǫ/4, and frequency f = η These equations are the source of Eq.
(2) in Section 2, where we give the final expressions for TTVs arising from the variation of λ 1 and λ 2 . In Eq. (2), we have taken the liberty to drop the star from Λ * j , but it is understood that n j , Λ j , and other quantities depending on these parameters in Eq. (2) are considered to be constant. Also, since D is a good proxy for the amplitude A Ψ , we replaced D → A Ψ in Eq. (2). Note that A Ψ is positive for Ψ 0 > Ψ eq and negative for Ψ 0 < Ψ eq .

Contribution from Eccentricity and Apsidal Longitude
We need to compute v 1 cos θ − u 1 sin θ and insert it in Eq. (51). As for the terms including θ, we obtain from Eqs. (49), (53) and (59) where We then use the following expansions to obtain cos θ and sin θ where J n (x) are the Bessel functions. Both these Fourier series converge rapidly for x ≪ 1. We therefore retain only the lowest-order harmonics. This leads to Here we neglected the coefficients C θ,2 that are of the order of ǫ (note that ǫ ≪ D; Fig. 11).
The harmonics with f in Eq. (64) appear from the resonant librations (Section 3.7). Since θ = φ 1 − ζ 1 (Appendix B), and φ 1 oscillates around π in the libration island, θ has the same circulation frequency as ζ 1 . Now, given that ζ 1 is defined from the apsidal longitudes ̟ 1 and ̟ 2 (Appendix B), the interpretation of f θ is that it is the mean precession frequency of the longitudes of periapsis. Neglecting terms O(ǫ) and using the definition of δ in Eq. (23) we find from (62) that where substituted Ψ 0 + D → Ψ eq (see Figure 11). The f θ frequency therefore scales in the same manner with planetary parameters as f = η 1/3 2 f τ (as expressed by the η 1/3 2 factor). Unlike f , which derives from f τ = 2π/P τ ∼ 3 in the libration island, f θ contains the factor 2(δ − Ψ eq ). This factor is ≪ 1 (compare the dotted and solid 'stable 1' lines in Figure 4). Therefore, f θ is substantially smaller than f , which shows that the variations from θ are expected to occur on a long timescale.
Specifically, the period P θ = 2π/f θ = π/(δ − Ψ eq ) is equal to 10.1 for δ = 1, 13.3 for δ = 2, and 18.1 for δ = 4. The longer periods for larger δ values are a consequence of Ψ eq approaching δ for increasing values of δ (Section 3.5). Note that P θ , at least in the approximation adopted here, is independent of the libration amplitude. Also, given that δ − Ψ eq < 0 in the libration island, the f θ frequency is negative as well, meaning that the circulation of θ is retrograde (implying retrograde rotation of ̟ 1 and ̟ 2 ).
The terms in Eq. (64) containing frequencies f θ and f could be combined together to produce harmonics with frequencies f θ ± f and f θ ± 2f . Given that, as we discussed above, the characteristic periods of these terms are largely different, we prefer to leave them multiplying each other in Eq. (64). Accordingly, Eq. (64) is interpreted as the resonant variations around the mean value that is slowly modulated with frequency f θ .
The expressions for u 1 and v 1 are derived from Eq. (54), after substituting Ψ from Eq. (59). After some algebra we obtain with coefficients In the above expression for the coefficient C u,0 , u 0 denotes the initial value  25) and (26). It shows that the approximation (66) is excellent for small libration amplitudes but loses precision for large libration amplitudes. This happens mainly because the terms O(ǫ 2 ) were neglected in Eq. (56). In principle, it should be possible to include these and higher order terms and improve the validity of the Fourier approximation by producing more general expressions. We leave this for future work. Now we should combine Eqs. (64) and (66) together. Unfortunately, this generates a very long expression for δY j . We do not explicitly give this equation here. The full expression was coded in a program and used to generate Figs. 1 and 2. We find that the TTV terms from v 1 cos θ − u 1 sin θ with the frequency f have amplitudes that are generally much smaller than the TTV amplitudes arising from the δλ j terms (Eq. 60). Here we therefore explicitly report only the most important harmonic with frequency f θ . These terms do not have a counterpart in Eq. (60). They are important for the long-term modulation of the TTV signal. Specifically, we find that This is the source of Eq. (5). Since J 0 (C θ,1 ) ∼ 1, we do not list this term in Eq. (5), where we also substitute D → A Ψ .

The Domain of Validity
We adopted several approximations in this work: I. The Laplacian expansion of the perturbing function used in Section 3.1 is convergent only if the planetary eccentricities are small enough (Sundman 1912). For a planet on a circular orbit, this limits the validity of the expansion to e < 0.25 (e < 0.2) for orbits near its inner (outer) 2:1 resonance, and to e < 0.15 (e < 0.12) for orbits near its inner (outer) 3:2 resonance. The analytic results derived here are not valid above these limits. See Nesvorný & Morbidelli (2008) for a discussion of Sundman's criterion.
II. All non-resonant terms were neglected in Section 3.1. The short-periodic terms with non-resonant frequencies produce short-periodic TTVs that can be calculated by the method described in Nesvorný & Morbidelli (2008). These terms can be linearly added to the expressions obtained here for resonant TTVs. The secular terms are second and higher orders in planetary eccentricities and contribute by only a small correction to the precession of orbits if the eccentricities are small.
III. The second-and higher-order resonant terms in planetary eccentricities were neglected in the perturbing function. This is an important approximation that limits the validity of the results to small eccentricities. The same assumption was adopted when writing δt j as a variation of orbital elements.
IV. The amplitude of the semimajor axis variations was assumed to be small (this allowed us to simplify the Keplerian Hamiltonian in Section 3.3). The same assumption was adopted to compute δt j from Eq. (1), where we also neglected all second and higher order terms in small variations of the orbital elements. We find that these approximations are generally valid and do not impose any meaningful limits on the range of planetary parameters where our analytic results are valid.
V. The exact solution of the second fundamental model of resonance was expanded in the Taylor series in ǫ and only the terms O(ǫ) were retained (Section 4.1). In addition, the Jacobi elliptic functions were written as the Fourier series and only the lowest harmonics were retained. Both these approximations limit the validity of our analytic TTV model to relatively small libration amplitudes. [In principle, the methods described in Section 4 can be used to obtain more general expressions.] Here we perform tests of these assumptions to establish the domain of validity of our analytic model. To this end we developed several codes that compute the resonant TTV signal at various stages of approximation. They are: A. A full N-body integrator of Eq. (7) and (8) where the gravitational interaction of planets is taken into account exactly. We used the symplectic code known as Swift (Levison & Duncan 1994) with routines for an efficient and precise determination of TTVs (Nesvorný et al. 2013; see also Deck et al. 2014).

B.
A numerical integrator in orbital elements that uses the Laplacian expansion of the perturbing function. Various terms can be included or excluded in this integrator. In the most basic approximation, the code includes only the first-order resonant terms from Eq. (10). Optionally, it also accounts for the second-order secular and/or resonant terms. This code is used to test the approximations II and III listed above C. A code that maps the initial orbital elements onto (22) and numerically integrates the corresponding Eqs. (25) and (26). Another code uses the exact analytic solution (31). As expected, these two codes give exactly the same result, which shows that our implementation of Eq. (31) is working correctly.

D.
A TTV code based on the analytic formulas derived in Section 4. This code is subject to all approximations discussed above. It cannot produce accurate results if the orbital eccentricities and/or libration amplitudes exceed certain limits.
We first test the approximation V. To this end we compare the results obtained with code B with the analytic results from method D. In B, we include the two first-order resonant terms and neglect terms that are the second or higher order in planetary eccentricities. 7 The masses and initial orbits are chosen such that δ = 2 (Fig. 13) or δ = 4 (Fig. 14). The initial orbits are then varied to survey different libration amplitudes. These tests show that the analytic method produces very reliable results for A Ψ 1 (top panels in Figs. 13 and 14). For the libration amplitudes much larger than that, our analytic expressions for u 1 and v 1 in Eq. (66) become inaccurate (see Fig. 12). As a consequence, the analytic approximation of the TTV terms from δh j fails (bottom panels in Figs. 13 and 14).
Interestingly, however, the analytic approximation of the full TTV signal is reasonable even if A Ψ > 1. This happens because the TTV terms from δλ j increase with A Ψ and become dominant for large A Ψ . We are able to reproduce these terms correctly, because the analytic formula in Eq.
(2) has more a general validity than the one that requires a correct approximation of the boomerang-shaped trajectories in the (u 1 , v 1 ) plane (Eq. 66). We therefore conclude that the analytic TTV model can be used, with some caution, even if the libration amplitudes are relatively large.
We now turn our attention to the approximations I, II and III. We find that the omission of the higher-order resonant terms in III is the most restrictive assumption. To illustrate this, Figures 15 (inner planet) and 16 (outer planet) show a comparison of the analytic model with TTVs computed from the N-body code (method A above). Here we set different planetary eccentricities ranging from e 1 = e 2 = 0.001 (left panels in both figures) to e 1 = e 2 = 0.05 (right panels). We find that that the analytic model works well for e 1 = e 2 = 0.001.
Already for e 1 = e 2 = 0.01, a significant discrepancy appears (see middle panels in Figs.  15 and 16). An important part of the discrepancy, however, is not due to the assumption III, but is rather related to the choice of initial conditions. Recall that, in addition to the resonant terms, the exact computation of TTVs with method A also contains the short-periodic harmonics, while the analytic method D does not account for these terms. This presents a difficulty when choosing the initial conditions in A and D that are consistent with each other. If the same values are adopted in A and D, the initial semimajor axes in A generate slightly different values of the mean orbital frequencies than the same initial semimajor axes in D. This effect then propagates into a difference in the libration frequency f . To demostrate this, we surveyed a small neigborhood of the initial conditions and found that it is always possible, if the eccentricities are sufficiently small, to apply a small adjustment such that the difference between the analytical and numerical results vanishes (left and middle bottom panels in Figs. 15 and 16). Note that this initial value problem does not seriously limit the application of the analytic model to the real data, because it requires only a very small adjustment of a 1 or a 2 (or equivalently n 1 and n 2 ), which can easily be absorbed by other parameters.
Another more fundamental discrepancy appears for e 1 = e 2 = 0.05. In this case, the TTV frequency computed from the analytic model is nearly 40% higher than the actual frequency, and the TTV amplitudes are ≃25% smaller than their actual values (right panels in Figs. 15 and 16). In this case, it is not possible to adjust the initial conditions to cancel the difference. This shows that the assumption III starts to fail. We confirm this by method B, where it becomes apparent that including the second-order resonant terms improves model's precision. Still, for e 1 = e 2 = 0.05 the amplitude discrepancy is relatively minor and can be compensated, for example, by a relatively small correction of planetary masses. We therefore find that the analytic model is still useful in this case. Our additional tests show that the analytic expressions for TTVs are not reliable for eccentricities exceeding ∼0.1. [The validity domain in e should be slightly larger for distant resonances such as 2:1, and smaller for k ≥ 4.] The analytic model was developed under the assumption of exactly co-planar planetary orbits. This assumption was used in Section 3.1 to neglect all terms in the Laplacian expansion of the perturbing function that depend on inclinations. The model is therefore not expected to be valid if the mutual inclination between orbits, I mutual , is large. We performed various tests of this assumption and found that the analytic model is reasonably accurate for I mutual < 10 • , but fails to produce reliable results for I mutual 10 • . This should not be a severe limitation of the applicability of the analytic results to the multi-transiting planetary systems, because the orbits in these systems are expected to be nearly co-planar (e.g, Fang & Margot 2012).

Conclusions
In this work we developed an analytic model for TTVs of a pair of resonant planets, and discussed how the TTV period and amplitude constrain the masses and orbits of the two planets. The model is strictly valid only for small orbital eccentricities (e < 0.1). It was developed under the assumption of co-planar orbits but our tests show that it is valid even if the mutual inclination of orbits is not large (< 10 • ).
The resonant TTV signal is expected to contain the harmonics of two basic periods: the period of resonant librations and the period of apsidal precession of orbits. The latter is expected to be ∼5 times longer than the former, and may be difficult to detect with a short baseline of the TTV measurements. The libration period is relatively insensitive to the exact location of the system parameters in the resonant island, and scales with (m/M * ) −2/3 . Its determination from the TTV measurements can therefore help to constrain the planetary masses. This is an important difference with respect to the near-resonant case (Lithwick et al. 2012), where the TTV period is the super period, which is independent of mass.
The TTV amplitudes, on the other hand, can be used to constrain the ratio of planetary masses m 1 /m 2 . Since both the TTV period and amplitude depend on the resonant amplitude A Ψ , some mild degeneracies between the mass and orbital parameters are expected, but these degeneracies can be broken by a detection of higher-order resonant harmonics, which constrain A Ψ , and/or short-periodic (chopping) effects. A detailed analysis of this problem and the application of our analytic model to the resonant exoplanets (Winn & Fabrycky 2015) is left for future work.
The work of DN was supported by NASA's ADAP program. DV was supported by the Czech Grant Agency (grant P209-13-01308S). Slawek Breiter pointed to us the general solution Eq. (31) in Whittaker & Watson (1920). We thank Katherine Deck, Eric Agol and an anonymous reviewer for useful comments on the manuscript.
The coefficients C n in Eq. (31) are obtained from the derivatives of f (Ψ) with Ψ = Ψ 0 The expressions for C 3 and C 4 are not needed if we set ψ 0 = π and thus f 0 = 0. The invariants in Eq. (32) can be written as where we denoted If we define then the three roots of the cubic equation can be obtained from In the equilibrium point, Ψ 0 = Ψ eq , we have that ∆ = 0 and thus F = g 1/3 3 /2. It follows that β = 0 and α = −F = −g 1/3 3 /2 = C 2 /24 < 0. Therefore, b = 0 and a = −C 2 /4, where a and b are defined in the main text. Figure 11 shows ǫ = b/a for δ = 3. It is zero at the equilibrium point and increases to ǫ ≃ 0.2 at the separatrix. Higher values of δ lead to smaller values of ǫ. For δ ≃ 1, on the other hand, ǫ can be as large as 0.6 near the separatrix.

B. Expressions for Y j
In Section 3 we omitted to explain one important issue that becomes apparent if the degree of freedom related to θ = kλ 2 − (k − 1)λ 1 is treated separately from those related to ̟ 1 and ̟ 2 (as it was done in BM13). To explain this issue we first define in direct correspondence to Eq. (17) where γ j = −̟ j . Second, we perform a transformation (see Section 3.4) to the new variables (V 1 , V 2 ; U 1 , U 2 ) And last, we introduce new polar variables (ζ 1 , ζ 2 ; Φ 1 , Φ 2 ) such that The inverse transformation to (B2) is With these definitions, it is straightforward to show that φ j = θ + ζ j , where φ j are the original angles defined in (19), ζ 1 = arg(AZ 1 + BZ 2 ) and ζ 2 = arg(BZ 1 − AZ 2 ), where Z j = X j + ıY j = 2Γ j exp ıγ j . We thus find that V 1 = v 1 cos θ − u 1 sin θ. When substituted into (B4), we obtain Eq. (51) in the main text. In addition, it can be shown that dζ 2 /dt = dφ 2 /dt − dθ/dt = 0. The angle ζ 2 is therefore constant. Consequently, since Φ 2 = const. as well, both U 2 and V 2 are constants of motion. This result is used in Section 4, where V 2 in Eq. (51) does not contribute to TTVs.
This preprint was prepared with the AAS L A T E X macros v5.2.
res.  Table 1: The coefficients f 1 and f 2 for different resonances. In the second column, we report the semimajor axis ratio for an exact resonance. The coefficient values are given for α = a 1 /a 2 = α res . Fig. 1.-A demonstration of the validity of the analytic TTV formulas obtained in this work. Here we set m 1 = m 2 = 10 −5 M * with M * = M Sun . The initial orbital elements were a 1 = 0.1 AU, a 2 = 0.13115 AU, e 1 = e 2 = 0.02, λ 1 = π, λ 2 = 0, ̟ 1 = 0, ̟ 2 = π. This orbital configuration corresponds to the libration regime in the 3:2 resonance (k = 3). The upper (lower) panel shows the results for the inner (outer) planet. The green and blue lines were computed by numerically integrating the differential equations corresponding to the resonant Hamiltonian (9) and (10). The black lines were obtained from the analytic TTV expressions (1), (2) and a generalization of (5) derived in Section 4.3. From left to right the panels show TTVs from δλ j and δh j , and their sum from Eq. (1). The validity of the analytic model is excellent in this case because the resonant amplitude is relatively small (A Ψ ≃ 0.65 with δ = 2.36 and Ψ 0 = 3.01; see Section 3).  Fig. 1 but with m 1 = m 2 = 3 × 10 −4 M * , M * = M Sun , a 1 = 0.1 AU, a 2 = 0.132 AU, e 1 = e 2 = 0.05, λ 1 = π, λ 2 = 0, ̟ 1 = 0, ̟ 2 = π. This corresponds to A Ψ ≃ −0.65, δ ≃ 1.38 and Ψ 0 ≃ 1.96 (Section 3). The periods are shorter in this plot than in Fig. 1, because the two planets were given larger masses.   (22). The solid line denoted 'stable 1' is the first stable equilibrium that exists for any value of δ. The second stable equilibrium, denoted by 'stable 2', appears only for δ > δ * ≃ 0.945. The dashed line is the unstable equilibrium. The gray area is the place where the resonant librations occur. The dynamical regime where the two orbits are just wide (narrow) of the resonance is labeled by α < α res (α > α res ). The dotted line is an approximation of the first equilibrium point, Ψ eq = δ, that becomes progressively better with increasing δ.    Fig. 4), respectively. The dotted lines are various approximations. The dotted lines denoted by 'a' is the super-period approximation with P τ = π/|δ|. The one denoted by 'b' is an oscillator approximation P τ = π(2/Ψ eq ) 1/4 , where Ψ eq of the libration point is computed exactly for each δ. The dotted line 'c' shows P τ = π(2/δ) 1/4 .  Here we assumed that P 1 = 10 days and P τ = 3 and computed P t = (νC 2 ) −1/3 P τ . Panels (a) and (b) show the results for the 2:1 and 3:2 resonances, respectively. This is the expected TTV period produced by the resonant librations in these resonances. The dashed lines show the approximation from Eq. (40), which becomes better for larger values of k (i.e., for α → 1). Fig. 10.-The validity domain of the L12 model. Here we assumed m 1 ≃ m 2 = m and small orbital eccentricities, and plotted the isolines of δ from Eq. (23) as a function of planetary mass and P 2 /P 1 − k/(k − 1). According to the discussion in the main text, the L12 model is valid for δ −2. This condition represents a combined constraint on the orbital period ratio and planetary masses. The parameter region below and to the right of the δ = −2 line is where L12's TTV formula is strictly valid. This region covers most of the period range shown here for planetary masses below that of Saturn. Fig. 11.-The dependence of various parameters of the analytic model on initial conditions. Here we set δ = 3, ψ 0 = π, and vary Ψ 0 . The solid line shows the small parameter ǫ = b/a. It is zero at the equilibrium point Ψ 0 = Ψ eq (labeled by the arrow) and increases to ǫ ≃ 0.2 near the separatrices. The dotted line shows the amplitude A Ψ = Ψ 0 − Ψ eq . The parameter D = C 1 /4a (Section 4.1), showed by the dashed line, is an excellent approximation of A Ψ for small amplitude librations.  The planetary masses and orbital parameters were chosen such that δ = 2 in all cases shown here. We fixed m 1 = m 2 = 3 × 10 −4 M * , M * = 1 M Sun , a 1 = 0.1 AU, and varied a 2 near the external 3:2 resonance with the inner planet (a 2 ≃ 0.131 AU). The eccentricities were adjusted to give δ = 2 and a desired initial value Ψ 0 . From top to bottom, the panels show the results for Ψ 0 = 2, Ψ 0 = 2.5 and Ψ 0 = 3.5. The stable equilibrium point is located at Ψ eq = 2.24 for δ = 2. The different cases shown here thus correspond to the libration amplitudes A Ψ = 0.25, 1.24 and 1.74 (from top to bottom). The plots show TTVs of the inner planet (the results for the outer planet are similar). The green dots were obtained by numerically integrating the differential equations corresponding to the resonant Hamiltonian (9) and (10). The black lines were obtained from the analytic TTV expressions (1), (2) and a generalization of (5) derived in Section 4.3. From left to right the different panels show the TTVs from δλ j and δh j , and their sum from Eq. (1). Fig. 14.-The same as Figure 13 but for δ = 4. In this case, the stable equilibrium point is at Ψ eq = 4.14. From top to bottom, different panels correspond to Ψ 0 = 3.5 (A Ψ = 0.67), Ψ 0 = 2.5 (A Ψ = 1.17) and Ψ 0 = 2 (A Ψ = 1.57). Fig. 15.-A comparison of TTVs obtained with the analytic model (black lines) and precise N-body integrator (green lines). From left to right, the orbital eccentricities of planets were increased to test the validity of the analytic model. We fixed m 1 = m 2 = 10 −5 M * , M * = 1 M Sun , a 1 = 0.1 AU, and a 2 near the external 3:2 resonance with the inner planet (a 2 = 0.13095 AU). TTVs of the inner planet are shown here. The case with e 1 = e 2 = 0.001 (left panels) corresponds to δ = 0.51 and Ψ 0 = 0.007. The case with e 1 = e 2 = 0.01 (middle panels) corresponds to δ = 1.25 and Ψ 0 = 0.75. The case with e 1 = e 2 = 0.05 corresponds to δ = 19.3 and Ψ 0 = 18.8. In the two bottom panels on the left, we illustrate how a small adjustment of the initial conditions improves the results (see the main text for a discussion). In the bottom-right panel, we rescaled time to show that a modest adjustment of the frequency can resolve the discrepancy when eccentricities are larger.