Exciting the transit timing variation phases of resonant sub-Neptunes

ABSTRACT

There are excesses of sub-Neptunes just wide of period commensurabilities like the 3:2 and 2:1, and corresponding deficits narrow of them. Any theory that explains this period ratio structure must also explain the strong transit timing variations (TTVs) observed near resonance. Besides an amplitude and a period, a sinusoidal TTV has a phase. Often overlooked, TTV phases are effectively integration constants, encoding information about initial conditions or the environment. Many TTVs near resonance exhibit non-zero phases. This observation is surprising because dissipative processes that capture planets into resonance also damp TTV phases to zero. We show how both the period ratio structure and the non-zero TTV phases can be reproduced if pairs of sub-Neptunes capture into resonance in a gas disc while accompanied by a third eccentric non-resonant body. Convergent migration and eccentricity damping by the disc drives pairs to orbital period ratios wide of commensurability; then, after the disc clears, secular forcing by the third body phase shifts the TTVs. The scenario predicts that resonant planets are apsidally aligned and possess eccentricities up to an order of magnitude larger than previously thought.

1 INTRODUCTION

Sub-Neptunes, planets with radii ≲ 4R_⊕, abound within an au of solar-type and low-mass stars (e.g. Fressin et al. 2013; Dressing & Charbonneau 2015; Petigura et al. 2018; Zhu et al. 2018; Sandford, Kipping & Collins 2019; Daylan et al. 2021; Otegi et al. 2021; Turtelboom et al. 2022). They are frequently members of multiplanet systems.

Fig. 1 shows the period ratios P₂/P₁ of neighbouring sub-Neptunes (subscript ‘1’ for the inner planet and ‘2’ for the outer), compiled from the NASA Exoplanet Archive. There appear to be excess numbers of pairs with commensurable periods, most obviously near P₂/P₁ = 3:2, and arguably also near 2:1 and 4:3. These ‘peaks’ in the histogram actually occur at P₂/P₁ values slightly greater than small-integer ratios. In terms of

which measures a pair’s fractional distance from a (q + 1):q commensurability (q is a positive integer), the observed 3:2, 2:1, and 4:3 peaks lie at |$\Delta \sim +1~{{\ \rm per\ cent}}$|⁠. Just adjacent to the peaks are population deficits or ‘troughs’ at |$\Delta \sim -1~{{\ \rm per\ cent}}$|⁠, with the clearest trough situated shortward of 2:1 (Lissauer et al. 2011; Fabrycky et al. 2014; Steffen & Hwang 2015). These peak-trough features can be discerned out to orbital periods of tens of days (Choksi & Chiang 2020).

The peaks and troughs are signatures of capture into mean-motion resonance (for an introduction, see Murray & Dermott 1999). Capture is effected by slow changes in a planet pair’s semimajor axes that bring period ratios nearer to commensurability. The migration is driven by some source of dissipation, e.g. the dissipation of tides raised on planets by their host stars (Lithwick & Wu 2012; Batygin & Morbidelli 2013; Millholland & Laughlin 2019), or the dissipation of waves excited by planets in their natal gas discs, a.k.a. dynamical friction (e.g. Goldreich & Tremaine 1980; Lee & Peale 2002). The 3:2 and 2:1 peak-trough statistics can be reproduced in a scenario, staged late in a protoplanetary disc’s life, whereby some sub-Neptune pairs migrate convergently into the peak (so that P₂/P₁ approaches commensurability from above), while others migrate divergently, out of the trough, because of eccentricity damping by dynamical friction (Choksi & Chiang 2020; MacDonald et al. 2020). The observed preference of first-order resonances for Δ > 0 is a consequence of dissipation driving systems to their fixed points, which are located intrinsically at Δ > 0 (Goldreich 1965; see also the introduction of Choksi & Chiang 2020).

The orbital periods in Fig. 1 are evaluated using the average time between a given planet’s transits, as measured from transit surveys like the Kepler and TESS space missions. The actual time between transits can vary from transit to transit – the transit timing variation (TTV) is the fluctuation about the average. The closer a pair of planets is to a period commensurability (i.e. the smaller is Δ), the larger are its TTVs (e.g. Agol et al. 2005; Holman & Murray 2005).¹ Many of the sub-Neptune pairs in the 3:2 and 2:1 resonant peaks exhibit simple sinusoidal TTV cycles; some examples from the Kepler TTV catalogue of Rowe et al. (2015) are illustrated in Fig. 2 (see also table 2 of Hadden & Lithwick 2017). Sinusoidal TTVs are expected for two bodies near commensurability (Lithwick, Xie & Wu 2012; Deck & Agol 2016; Hadden & Lithwick 2016; Nesvorný & Vokrouhlický 2016a), with the TTVs of the inner and outer members varying over the same TTV ‘super-period’, P_super ∼ P₁/|Δ|. For the examples shown in Fig. 2, the inner and outer planet TTVs are anticorrelated: When one rises, the other falls. The TTV amplitudes depend on planet masses and eccentricities (Lithwick et al. 2012; Wu & Lithwick 2013; Hadden & Lithwick 2014, 2016, 2017), and for the systems of interest here are on the order of an hour (Mazeh et al. 2013; Rowe et al. 2015).

A TTV cycle is described not only by its amplitude and period but also by its phase (Lithwick et al. 2012; Deck & Agol 2016). A simple working definition for the TTV phase starts by identifying a TC. A TC is a conjunction (occurring when the star and the two planets lie along a line) that coincides with a transit. Equivalently, TCs occur when both planets transit simultaneously. Such times are marked ‘TC’ in Fig. 2; they occur once per super-period. The phase Φ_TTV is the phase interval between a TC and the nearest TTV zero-crossing. We use the ‘descending’ zero-crossing (when the TTV crosses zero from above) when evaluating the phase Φ_{TTV, 1} of the inner member of a resonant pair, and the ‘ascending’ zero-crossing for the outer phase Φ_{TTV, 2}. For anticorrelated TTV pairs like those in Fig. 2, this convention implies Φ_{TTV, 1} = Φ_{TTV, 2},² in which case a single phase Φ_TTV uniquely describes the system. The resonant pair featured in the middle column of Fig. 2 has Φ_TTV = 0 for both planets because the TCs coincide with the TTV zero-crossings. If the TC precedes the nearest zero-crossing, then Φ_TTV < 0 (left column of Fig. 2); otherwise Φ_TTV > 0 (right column).

The TTV phase can be interpreted geometrically. Fig. 3 (left-hand panel) illustrates that Φ_TTV = 0 corresponds to a TC where the inner planet is at periapse (‘peri’) while the outer planet is at apoapse (‘apo’). The symmetry of this orbital configuration about both the line of conjunctions and the line of sight leads to a TTV that is momentarily zero. This symmetry is broken for Φ_TTV ≠ 0. For the case of non-zero phase, during the TC, the inner planet is displaced from its peri, and the outer planet is displaced from its apo (right-hand panel of Fig. 3).

Interpreted geometrically as such, Φ_TTV has physical significance. The symmetric peri/apo configuration at conjunction (leading to Φ_TTV = 0) represents a fixed point of a first-order resonance (e.g. Peale 1986). As we mentioned earlier, dissipation leading to resonance capture drives the system toward this fixed point at Δ > 0. Thus the various dissipative mechanisms – disc migration, disc eccentricity damping, and tides – that reproduce the observed period ratio statistics (Fig. 1) predict Φ_TTV = 0.³

Is this expectation of zero TTV phase confirmed by the observations? The answer is no – Fig. (4a) plots Φ_{TTV, 1} and Φ_{TTV, 2} for a sample of sub-Neptune pairs near 2:1 and 3:2 resonances. Histograms of these phases are plotted in Fig. (4b), which does not bother to distinguish Φ_{TTV, 1} from Φ_{TTV, 2} since Fig. (4a) shows they are nearly equal. The observed distribution of Φ_TTV runs the gamut from −180° to +180°, in apparent contradiction to scenarios of dissipative, migration-induced resonance capture that predict Φ_TTV = 0. This behaviour does not seem to change with distance from the host star, as shown in Fig. (4c).

Our goal in this paper is to resolve this discrepancy – to find a formation history for sub-Neptunes that can explain both their period ratio statistics (Fig. 1) and their non-zero TTV phases (Fig. 4). We will continue to adhere to the picture in which a residual protoplanetary disc drives a pair of sub-Neptunes into resonance while damping their eccentricities (Choksi & Chiang 2020; MacDonald et al. 2020). To this scenario we will add a third, non-resonant planet. The extra body might be another sub-Neptune. Or the extra body could be a giant planet – a ‘cold Jupiter’ or ‘sub-Saturn’ – that radial velocity surveys have found in ∼40  per cent of systems hosting a transiting sub-Neptune (Bryan et al. 2016; Zhu & Wu 2018; Mills et al. 2019; Rosenthal et al. 2021).

We will show that a third, non-resonant planet can excite the TTV phase of a resonant pair without interfering with P₂/P₁. The excited phase is not constant, but varies as the third body drives secular changes in the pair over a time-scale P_sec  ≳ P₁/μ₃ ∼ (10³ − 10⁵)P₁, where μ₃ = m₃/M_⋆ is the third planet’s mass ratio with the star. This secular time is much longer than the TTV cycle super-period, P_super ∼ 10²P₁. Over a typical observing window spanning just a few cycles, the TTVs of a resonant pair masquerade as a standard two-planet sinusoid with constant phase.

This paper is organized as follows. Section 2 studies a representative three-planet system. Section 3 explores how TTV phases scale with the third planet’s mass, eccentricity, and orbital period. Section 4 summarizes and discusses. Technical details including the equations solved in this paper are contained in the appendices.

2 RESONANT DYNAMICS INCLUDING A NON-RESONANT COMPANION

We consider a coplanar pair of sub-Neptunes (labelled ‘1’ for the inner body and ‘2’ for the outer) that start outside the 2:1 mean-motion resonance. Exterior to this pair orbits another coplanar planet (‘3’), situated far from any first-order resonance with the pair. At the start of our calculations, we apply dissipative forces to the interior pair to smoothly damp their eccentricities and have their semimajor axes converge. Once the pair captures into resonance and equilibrates, we shut off these external forces and continue evolving all three planets without dissipation.

We numerically integrate Lagrange’s equations (A1)–(A18) for the mean motions n(t), eccentricities e(t), and longitudes of periapse ϖ(t) of the three planets. We feed these numerical solutions into the analytic formulae of Lithwick et al. (2012) – their equations (1)–(13) – to compute the TTVs of planets 1 and 2, including TTV phases. Details are provided in Appendix  A, which also checks our semi-analytic approach against an N-body simulation.

2.1 Results for period ratios and TTV phases

The left column of Fig. 5 reviews how the inner pair evolves in the absence of the third body (e.g. Lee & Peale 2002; Choksi & Chiang 2020). Two sub-Neptunes with masses |$m_1 = m_2 = 7.5\, \text{M}_{\hbox{$\oplus $}}$| convergently migrate toward the 2:1 resonance. They are captured into resonance, with both resonant arguments ϕ₁ and ϕ₂ librating (panel e), and Δ equilibrating to  +0.9  per cent (panel a). The TTV phases rapidly damp to ≲10° from the applied dissipation (panel b). The eccentricities level out as the imposed damping balances excitation by resonant migration (panel c). After the dissipation shuts off, the eccentricities oscillate a bit more, but overall the system properties do not change much.

The right column of Fig. 5 shows how the same pair evolves when accompanied by a third planet of mass m₃ = 1M_J, period P₃ = 20P₁, and eccentricity e₃ = 0.3. While disc torques are still being applied, the evolution of planets 1 and 2 is not too different from the no-3 case. But after disc torques are removed, the resonant pair is free to respond in full to planet 3. Eccentricities e₁ and e₂ are secularly driven by planet 3 (panel h), which also forces apsidal alignment, ϖ₁ ≃ ϖ₂ (panel i). The resonant locks are broken (panel j). Despite all this, the period ratio P₂/P₁ is left at the same equilibrium value as in the two-planet case; compare panels a and f. The main effect of the third planet on Δ is to induce modest, ∼5  per cent oscillations. The oscillations have both a short period equal to the TTV super-period |$P_{\rm super} \approx 500 n_1^{-1}$| (panel f inset) and a long period corresponding to the secular forcing period |$P_{\rm sec}~\approx 10^6\, n_1^{-1} \approx 2 \times 10^3 P_{\rm super}$|⁠. It is not surprising that Δ is relatively unaffected by planet 3, since secular interactions do not change semimajor axes.

By contrast to Δ that is modulated only mildly by the third body, the TTV phases are dramatically affected. Over P_sec, both phases Φ_{TTV, 1} and Φ_{TTV, 2} are freed from zero and cycle between −180° and 180° (panel g).

2.2 Analysis of TTV phases

The TTV phases can be dissected using the analytic formulae derived by Lithwick et al. (2012) for a pair of planets near a first-order commensurability⁴:

Expressions for the resonantly forced eccentricities e_{forced, 1} and e_{forced, 2} and the order-unity coefficients f₁ and f₂ are given in Appendix  A; we note here that to order-of-magnitude, e_{forced, 1} ∼ μ₂/Δ, and vice-versa. The angle λ^q = (q + 1)λ₂ − qλ₁, where λ is the mean longitude, measures the longitude of the line of conjunctions between planets 1 and 2. It sweeps at constant rate dλ^q/dt = −Δqn₁. Our convention is that the observer sits along the positive x-axisso that λ^q = 0 during a TC. If the forced term proportional to e_{forced, 1}sin λ^q in equation (2) were dominant, then Φ_{TTV, 1} = 0 (since every λ^q = 0 TC would give TTV₁ = 0). This forced portion of the TTV arises because the orbit precesses apsidally, at a rate dϖ₁/dt = dλ^q/dt = −Δqn₁ when the resonant argument ϕ₁ = λ^q − ϖ₁ is fixed at 0 in resonance lock. Analogous statements apply to the purely forced term e_{forced, 2}sin (λ^q − π) in equation (3).

Adding the terms proportional to |$\mathcal {Z}$| leads to Φ_TTV ≠ 0. These terms depend on the free component of each planet’s eccentricity, i.e. the portion of the eccentricity beyond the resonant fixed point value. See the top panel of Fig. 6 for how e_free and the related angle ϖ_free are extracted from the osculating e and ϖ. The free elements e_free and ϖ_free define a free eccentricity vector e_free which is constant in time; it is a kind of integration constant determined by initial conditions. The forced eccentricity vector’s magnitude is also constant, but its direction varies with time at rate dλ^q/dt. Thus if the free eccentricity is non-zero, then the osculating eccentricity, given by the vector sum of the forced and free eccentricities, oscillates about the forced (fixed point) value (see also Lithwick et al. 2012, fig. 1). Such eccentricity variations drive mean-motion variations (via the Brouwer integrals associated with the resonance; equation A18 of Lithwick et al. 2012), which in turn contribute to the TTV.

These mean-motion TTVs due to free eccentricities are generally phase-shifted relative to the aforementioned apsidal precession TTVs. The bottom panel of Fig. 6 illustrates this phase shift geometrically. From equation (4), we construct the difference vector e_{free, 1} − |f₂/f₁|e_{free, 2}. The direction of this vector gives Φ_TTV, in the limit that the magnitude of this vector ≫Δ so that the mean-motion TTVs dominate.

Returning to our disc-driven capture scenario, after the disc clears, the third body perturbs the inner pair off the fixed point, secularly exciting free eccentricities (panel h of Fig. 5), which lead to mean-motion changes (panel f), and by extension non-zero TTV phases (panel g). When total eccentricities e₁ and e₂ are excited up to several per cent – an order of magnitude larger than their forced values (panel c) – they are dominated by their free components, i.e. e_free ≃ e and ϖ_free ≃ ϖ. In addition, the two resonant planets apsidally align under the influence of the third body, with ϖ₁ ≃ ϖ₂ varying about ϖ₃ over a secular time-scale (panel i; see also Beaugé, Michtchenko & Ferraz-Mello 2006 and Laune, Rodet & Lai 2022).⁵

The conditions e_free ≫ e_forced and ϖ₁ = ϖ₂ inserted into equations (2)–(4) imply that each planet’s TTV phase lies between 0 and ϖ. Where the phase lands in this range is controlled by the ratio

obtained by comparing the magnitudes of the terms in the square brackets in equations (2)–(3), i.e. the magnitude of the phase-shifted mean-motion TTV to that of the zero-phase precessional TTV. As Fig. 7 shows, when |$\mathcal {D}/\Delta \ll 1$|⁠, both planets have zero TTV phase, whereas when |$\mathcal {D}/\Delta \gg 1$|⁠, Φ_TTV = ϖ. For our sample evolution, |$\mathcal {D}/\Delta$| stays above unity over the secular cycle (Fig. 8a), and the TTV phases approximately track ϖ (Fig. 8b).

Externally amplifying the eccentricities of a resonant pair of planets excites Φ_TTV more readily than it changes Δ. Whereas Φ_TTV/ϖ is of order |$\mathcal {D}/\Delta$| (for |$\mathcal {D}/\Delta \lesssim 1$|⁠), the fractional change in Δ is of order |$\sim 10e_{\rm forced} \mathcal {D}/\Delta$|⁠, by virtue of the Brouwer constants of the motion associated with the resonance (equation A19 of Lithwick et al. 2012). The relative insensitivity of Δ means that the fits obtained to observed period ratios within the scenario of disc-driven dissipation and migration (Choksi & Chiang 2020) should not change much when eccentricity forcing by an external companion is added.

3 PARAMETER DEPENDENCE

We explore how the non-resonant third planet’s eccentricity e₃, mass m₃, and orbital period P₃ affect the TTV phases of the resonant pair. As we saw in the previous section, the larger is |$\mathcal {D}/\Delta$| [equation (5)], the more the TTV phases can deviate from zero. When |$\mathcal {D}/\Delta \gt 1$|⁠, the TTV phases approximately track ϖ₁ or ϖ₂ – these apsidal longitudes are nearly the same under secular driving by an eccentric third body. The pair’s apses in turn follow ϖ₃, which is free to take any value between 0 and 2π relative to our line of sight. Thus as a proxy for the range of accessible TTV phases, we measure how |$\mathrm{max}(\mathcal {D}/\Delta)$| depends on the properties of the third planet. To build intuition, we first study the test particle limits in which the resonant planets have very different masses, m₁ ≪ m₂ ≲ m₃ or m₂ ≪ m₁ ≲ m₃. Then we see how our results change when m₁ is comparable to m₂, keeping the total mass fixed to |$m_1 + m_2 = 15\, \text{M}_{\hbox{$\oplus $}}$|⁠.

Fig.9 plots |$\mathrm{max}(\mathcal {D}/\Delta)$| as measured from our numerical integrations varying one parameter at a time around a fiducial set |$\lbrace e_3,\, m_3,\, P_3\rbrace = \lbrace 0.3,\, 1\, M_{\rm J},\, 20P_1\rbrace$|⁠. In both test particle limits (purple squares and green crosses) near the 2:1 resonance (top row), |$\mathrm{max}(\mathcal {D}/\Delta)$| spans a factor of ∼100 across our explored parameter space. Most of this variation stems from changes in |$\mathcal {D}$| (and not Δ; see the end of Section 2.2), reflecting how strongly the third body secularly amplifies e₁ and e₂. Reducing the perturber’s eccentricity (left-hand panels of Fig. 9) or moving it farther from the pair (right-hand panels) weakens this amplification. At the same time, the perturber’s mass seems to hardly matter (middle panels); in the test particle limit, perturbers weighing anywhere from 15 M_⊕ to a few M_J produce nearly the same |$\mathrm{max}(\mathcal {D}/\Delta)$|⁠. Pairs near the 3:2 resonance (bottom row) follow all the same trends, but their |$\mathrm{max}(\mathcal {D}/\Delta)$| values are systematically lower than for the 2:1. The TTV phases of 3:2 pairs are harder to excite by an external perturber than those of 2:1 pairs, a consequence of 3:2 pairs being closer together and enjoying a stronger resonant interaction that better isolates them from external influences. More details about this difference between the 3:2 and 2:1 are given in Appendix  B.

Our numerical results in the test particle limits suggest power-law scalings that can be reproduced as follows. When m₁ ≪ m₂, the eccentricity e₂ is forced by m₃ and not m₁. If the perturber is distant (P₃ ≫ P₂) and massive (m₃ ≫ m₂), it secularly excites e₂ up to max(e₂) = (5/2)(P₂/P₃)^2/3e₃ (Murray & Dermott 1999). If we further assume that Δ is constant and that e₁ tracks e₂, then

scalings that are consistent with our numerical data in the appropriate limits for both the 2:1 and 3:2. The same scalings apply when m₂ ≪ m₁ ≪ m₃.

In the more realistic case where m₁ and m₂ are comparable (blue circles in Fig. 9), |$\mathrm{max}(\mathcal {D}/\Delta)$| and by extension TTV phases drop relative to the m₁ ≪ m₂ test particle case, precipitously in some regions of parameter space. Allowing for m₁ ∼ m₂ strengthens the mutual interaction of the resonant pair and shields them from perturbations by an external third body – see also Fig. 10, which underscores the sensitivity of |$\mathrm{max}(\mathcal {D}/\Delta)$| to perturber mass when m₁ ∼ m₂. To keep |$\mathrm{max}(\mathcal {D}/\Delta)\gt 1$| for the 3:2 resonance – which automatically keeps |$\mathrm{max}(\mathcal {D}/\Delta)\gt 1$| for the 2:1 – it suffices to have an eccentric Jupiter-mass perturber at P₃ ∼ 15P₁; or an eccentric Saturn-mass perturber with P₃ ∼ 10P₁; or an eccentric Neptune-mass perturber with P₃ ∼ 5P₁. Interestingly, in some regions of parameter space for the 3:2, |$\mathrm{max}(\mathcal {D}/\Delta)$| actually increases for m₁ = m₂ relative to the m₂ ≪ m₁ test particle case (Fig. 9), but the enhancements are less than factors of 2–3.

So far in this section, we have used |$\mathrm{max}(\mathcal {D}/\Delta)$| as a proxy for whether TTV phases can be large or must remain small. In reality |$\mathcal {D}/\Delta$| and by extension Φ_TTV/ϖ vary with time over the secular cycle (Figs 5 and 8). Time-sampled distributions of Φ_TTV generated from ensembles of three-planet numerical integrations are plotted in Fig. 11. For each integration, set up the same way as in Section 2 and with an initial ϖ₃ drawn randomly from a uniform distribution between 0 and 2π (reflecting the isotropy of space), we record Φ_{TTV, 1} and Φ_{TTV, 2} at evenly spaced time intervals over one secular cycle post-dissipation. Each panel in Fig. 11 compiles time samples from 30 such integrations, each with a different starting ϖ₃, and combining Φ_{TTV, 1} and Φ_{TTV, 2} because they are similar. The resultant distributions of Φ_TTV concentrate at zero, more strongly for the 3:2 since it is less sensitive to secular forcing by the third body. The system parameters for the left- and right-hand panels are chosen to give practically identical Φ_TTV distributions, illustrating the degeneracy between e₃, m₃, and P₃.

4 SUMMARY AND DISCUSSION

Many mean-motion resonant sub-Neptunes exhibit sinusoidal TTVs whose phases Φ_TTV are non-zero, consistent with circulation or large-amplitude libration about their fixed points. This observation is surprising because dissipative processes that capture pairs into resonance also drive them to their fixed points, damping Φ_TTV to zero. In this paper, we showed how secular forcing of a resonant pair by an eccentric third body can phase-shift TTVs. The third body can force the pair into apsidal alignment, and Φ_TTV can track their common longitude of pericentre. While the non-resonant third body knocks the pair off the resonant fixed point, it does not much alter their period ratio, P₂/P₁. Thus scenarios that do alter P₂/P₁ to reproduce the observed ‘peaks’ and ‘troughs’ in the period ratio histogram can be simply augmented by a third-body secular perturber to account for non-zero TTV phases. In particular, orbital migration and eccentricity damping in a residual protoplanetary disc can drive a fraction of sub-Neptunes into resonance with P₂/P₁ just wide of perfect commensurability, as observed (Choksi & Chiang 2020); and if such capture plays out in the presence of a sufficiently massive and eccentric non-resonant third body, the TTVs can have non-zero phases after the disc clears, also as observed.

Fig.12 summarizes the parameter space favouring the secular excitation of TTV phases and shows where known companions to TTV pairs lie in this space. Two such companions have had their eccentricities e₃ measured from radial velocities (RVs). One of them, a sub-Saturn with |$e_3 = 0.13^{+0.13}_{-0.09}$|⁠, looks marginally capable of exciting Φ_TTV, which is measured to be ∼20°. The other, a super-Jupiter with |$e_3 = 0.20^{+0.01}_{-0.01}$|⁠, seems too distant to generate the ∼100° phase of its associated TTV pair. Most of the remaining companions in Fig. 12 are other transiting sub-Neptunes lacking measured eccentricities. We see that their eccentricities e₃ would have to be ≳0.05–0.15 to secularly excite TTV phases. How such eccentricities were acquired and maintained against disc damping would need to be explained.

Fig. 12 also shows that many of these smaller-mass third-body perturbers lie near period commensurabilities with their resonant pairs. The perturbers have Δ’s of ∼5–20 per cent, larger than the Δ’s of ∼1  per cent exhibited by the pairs. In preliminary N-body integrations of these systems, we have not found near-resonant forcing by the third body to affect the TTV phase of the resonant pair. We do find that when the third body’s Δ ≲ the pair’s Δ, the TTVs of the pair are no longer sinusoidal. Non-sinusoidal TTVs are observed for extensive resonant chains like TRAPPIST−1, Kepler-60, and TOI-1136 (e.g. fig. 2 of Agol et al. 2021; fig. 10 of Jontof-Hutter et al. 2016; fig. 5 of Dai et al. 2023). Although a phase seems challenging to define for a non-sinusoidal signal, one can always mark when two (or more) planets transit the star at the same time (as we have done in Fig. 2). Future work could explore how the relative timing of TCs in resonant chains depends on formation history, including how much eccentricity and semimajor axis damping each planet experienced.

TTVs can be combined with transit duration and/or RV data to place tighter dynamical constraints on systems. Radial velocities can pin down planet masses, and transit durations can help measure individual planet eccentricities, breaking the ‘mass–eccentricity’ and ‘eccentricity–eccentricity’ degeneracies inherent to using TTVs alone (Lithwick et al. 2012). For example, Dawson et al. (2021) combined RV and transit data to measure a 60° ± 2° libration amplitude for the resonant argument ϕ of the sub-Neptune TOI-216b in 2:1 resonance with the gas giant TOI-216c. From our study, we have learned that a third body can secularly excite resonant libration amplitudes, all the way up to circulation, without materially changing the period ratio at the time of resonance capture. For the resonant planet pairs Kepler-90gh and K2-19bc, Liang, Robnik & Seljak (2021) and Petigura et al. (2020) measured aligned apsides, which our study has shown is possible if these pairs are secularly forced by third bodies. The apsidal alignment in Kepler-90 might be enforced by any of the six other transiting planets in the system. For K2-19, we find that an eccentric |$30\, \text{M}_{\hbox{$\oplus $}}$| planet orbiting with P₃ = 5P₁ ≈ 40 d could explain the observed alignment. With an RV semi-amplitude of 7 m s⁻¹, such a body could hide below the current noise threshold of 13 m s⁻¹. Continued RV monitoring will decide whether an eccentric perturber exists, or if the observed alignment requires another explanation such as eccentricity pumping by the natal disc (Laune et al. 2022).

Wu & Lithwick (2013) and Hadden & Lithwick (2014) used the observed distribution of TTV phases to estimate the free eccentricities of resonant sub-Neptunes, assuming that the planets in a given resonant pair have uncorrelated apsidal orientations and eccentricities. This assumption does not hold in our scenario where both planets are secularly forced by a third body into apsidal alignment and to have a particular ratio of eccentricities. Whereas they inferred free eccentricities of 1–2  per cent, we would derive, in the context of our model, free eccentricities as large as the planets’ total eccentricities, which can be secularly forced by an eccentric third body up to values of ∼0.2. Do severely underestimated free eccentricities affect the distribution of planet masses they derived using TTV amplitudes – a distribution known to be consistent with independent mass measurements from RVs? Fortunately, no – the TTV amplitudes depend on free eccentricities only through the particular linear combination of free eccentricity vectors shown in the bottom panel of Fig. 6, and this vector combination is reliably measured from the observed TTV phase.

We have assumed in this work that gas disc torques always drive planet pairs to their resonant fixed points. If the disc is turbulent, density fluctuations add a random walk component to a planet’s orbital evolution (Laughlin, Steinacker & Adams 2004; Adams, Laughlin & Bloch 2008; Rein & Papaloizou 2009; Batygin & Adams 2017). Planetesimal discs also drive migration having both smooth and stochastic contributions (Murray-Clay & Chiang 2006; Ormel, Ida & Tanaka 2012; Nesvorný & Vokrouhlický 2016b). On the one hand stochasticity can excite free eccentricities and TTV phases (Goldberg & Batygin 2022). But too much stochasticity can wipe out the peak-trough asymmetries seen in the Kepler period ratios (Fig. 1; Lissauer et al. 2011; Fabrycky et al. 2014; Steffen & Hwang 2015; Choksi & Chiang 2020). Rein (2012) argues that a suitably tuned mix of stochastic and smooth migration torques can reproduce the observed period ratio distributions. This study’s reproduction of resonant peak-trough structures would be better assessed by replacing their cumulative distributions with differential ones – the modeled troughs short of commensurability might actually be too filled in. What the stochasticity and migration parameters advocated by Rein (2012) predict for TTVs and in particular TTV phases is unknown; it is also not clear that the implied period distribution (distinct from the period ratio distribution) satisfies observations (cf. Lee & Chiang 2017). Goldberg & Batygin (2022) neglect the problem of the peak-trough asymmetries in the period ratios, as their plots record |Δ| instead of Δ.

ACKNOWLEDGEMENTS

We are indebted to Yoram Lithwick for helping to launch this work. We also thank Eric Agol, Konstantin Batygin, Lister Chen, Rebekah Dawson, Dan Fabrycky, Eric Ford, Max Goldberg, Sam Hadden, Matt Holman, Shuo Huang, Sarah Millholland, Erik Petigura, Hanno Rein, Jason Rowe, Andrew Vanderburg, Shreyas Vissapragada, and Jack Wisdom for useful exchanges. An anonymous referee provided a constructive report. Simulations were run on the Savio cluster provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Financial support was provided by NSF AST grant 2205500, and an NSF Graduate Research Fellowship awarded to NC.

DATA AVAILABILITY

No new data were collected as part of this work.

Footnotes

References

APPENDIX A: EQUATIONS SOLVED

The results in this paper were obtained by integrating Lagrange’s equations of motion, using a disturbing function that includes resonant and secular terms to leading order. We list these equations in Section A1 and check their validity in Section A2.

A1 Equations of motion

We consider two planets lying near a (q + 1):q mean-motion resonance (where q is a positive integer) and accompanied by a third non-resonant planet. The equations below are for the case where the third planet orbits exterior to the near-resonant pair. We also modeled third bodies lying interior to the pair, but do not list the equations for that case; they may be derived in a straightforward way from the exterior case (Murray & Dermott 1999).

The inner resonant planet evolves according to:

where n is the mean motion, e is the eccentricity, ϖ is the longitude of periapse, ϕ_i = (q + 1)λ₂ − qλ₁ − ϖ_i is the resonant argument describing where conjunctions happen in orbital phase (i.e. relative to periapse), λ is the mean longitude, μ is the planet-to-star mass ratio, and subscripts 1, 2, and 3 refer respectively to the inner resonant planet, the outer resonant planet, and the third exterior non-resonant planet. The coefficients f₁, f₂, C_A, and C_B depend on the semimajor axis ratio α_ij = a_i/a_j = (P_i/P_j)^2/3:

(Murray & Dermott 1999). The term δ_{q, 1} in (A6) is the Kronecker delta and accounts for the contribution of the indirect potential of the 2:1 resonance; we will return to this term in Appendix  B. We keep the coefficients (A5)–(A8) fixed at their initial values (Table A1) since they vary negligibly over the course of our integrations.

Terms involving t_e and t_a describe eccentricity damping and orbital migration, respectively, due to disc torques. We assume that eccentricity damping alone conserves the planet’s orbital angular momentum by setting the coefficient p = 3 (see also section 2.1 of Goldreich & Schlichting 2014). To ensure convergent migration, we set t_a < 0 for planet 1, and do not allow planet 2 to migrate. In reality both planets probably migrate, but their mutual interaction is controlled by the relative migration rate.

The outer resonant planet obeys an analogous set of equations:

As mentioned above, we have removed the migration term involving t_a for planet 2 since outcomes depend only on the relative migration rate. For simplicity, the eccentricity damping time t_e is assumed equal for planets 1 and 2.

We set |$t_a = -10^7\, n_1^{-1}$| and remove all terms due to disc torques once the pair has equilibrated in resonance. For our simplifying assumptions, the final equilibrium state depends on the ratio t_e/|t_a| (Goldreich & Schlichting 2014; Terquem & Papaloizou 2019; Choksi & Chiang 2020). A common estimate is t_e/|t_a| ∼ (h/a)² ∼ 10⁻³ for disc aspect ratios h/a ∼ 0.03 at a ∼ 0.3 au (Tanaka, Takeuchi & Ward 2002; Tanaka & Ward 2004; Cresswell & Nelson 2008). This estimate for t_e/|t_a| derives from idealized discs with power-law surface density profiles and constant temperatures. Accounting for more realistic disc properties, including gaps, horseshoe substructures, and circumplanetary material (e.g. D’Angelo, Kley & Henning 2003; D’Angelo, Bate & Lubow 2005; Masset, D’Angelo & Kley 2006; Duffell & Chiang 2015), and non-isothermal equations of state (Kley & Nelson 2012, their fig. 3), can change both t_e and t_a.

Given these complications, we do not rely on the disc’s aspect ratio h/a to set t_e/|t_a|. Instead, we use equations (49)–(51) of Terquem & Papaloizou (2019) to choose, for each integration, the value of t_e/|t_a| that yields an equilibrium Δ (fractional distance from period commensurability) of 1 per cent, close to the values observed for Kepler sub-Neptunes. The precise value depends on the individual planet masses and the resonance being considered. For reference, when μ₁ = μ₂ = 3 × 10⁻⁵, we calculate that t_e/|t_a| = 1 × 10⁻⁴ for the 3:2 resonance and 4 × 10⁻⁵ for the 2:1. These are similar to the values advocated by Huang & Ormel (2023, fig. 7). Our results for TTV phases do not seem especially sensitive to t_e/|t_a|. We find for the example model in Fig. 5 that changing t_e/|t_a| from its nominal value by a factor of 10 changes the TTV phase proxy |$\mathrm{max}(\mathcal {D}/\Delta)$| by a factor of 2.

The resonantly forced eccentricities used in equations (2) and (3) to determine TTVs analytically are

obtained by taking |$\dot{\phi }_1 = \dot{\phi }_2 = 0$| while neglecting all but the resonant interaction terms.

Finally, we account for the secular back-reaction of planets 1 and 2 on to planet 3:

Note that we neglect eccentricity damping and migration of planet 3. To generate non-zero TTV phases for planets 1 and 2, we need planet 3’s orbit to be sufficiently eccentric, but do not model explicitly how planet 3 acquired or retained its eccentricity.

We numerically integrate the 10 coupled differential equations in (A1)–(A18) using the lsoda algorithm implemented in the solve-ivp module of scipy (Virtanen et al. 2020). The fractional tolerance of the solutions is set to 10⁻⁹.

A2 Checks

Equations (2)–(4) for the TTVs of a near-resonant pair of planets are derived from Lithwick et al. (2012) who assumed that the planets’ mean longitudes do not deviate from linear ephemerides by more than ∼1 rad (see their equations A7 and A22, and their appendix A.2). This requires Δ, the pair’s distance from a period commensurability, to be not too small. In our paper, we have focussed on resonant pairs of planets forced by third bodies to have free eccentricities ≫ forced eccentricities, and to be apsidally aligned with ϖ₁ ≃ ϖ₂; under these conditions, the assumption of Lithwick et al. (2012) translates to |$\Delta \gtrsim \sqrt{\mu |e_1 -|f_2/f_1|e_2|}$|⁠, where μ is either μ₁ or μ₂. Our modeled sub-Neptunes have |$\Delta \sim 1 {{\ \rm per\ cent}}$|⁠, e ≲ 0.1, and μ ∼ 3 × 10⁻⁵, and thus satisfy this requirement.

As an added check, we integrate the same representative system shown in Fig. 5 using the REBOUND N-body code (Rein & Liu 2012). We use the WHFast (Rein & Tamayo 2015) implementation of the Wisdom & Holman (1991) algorithm with a timestep of 0.02P₁, where P₁ is the orbital period of the innermost planet. Migration and eccentricity damping forces are included using the modify-orbits-forces package in REBOUNDx (Tamayo et al. 2020). Simulation outputs are recorded every 500P₁ using a SimulationArchive (Rein & Tamayo 2017). We compute TTV phases by rerunning segments of our simulation with finer timestepping. For each snapshot in the SimulationArchive we integrate forward by 300P₁ and record transit times to a precision of 10⁻⁷P₁ using the bisection algorithm suggested in the REBOUND documentation.⁶ TCs are identified using the same method applied to the observations in Fig. 2, and used to compute TTV phases. If no TCs occur in the given segment, we move on to the next simulation snapshot. We have checked that we obtain nearly identical results for the TTV phases if instead of using TCs as our reference we use the times when the two planets have true longitudes that differ from zero by less than 0.1 rad.

Fig.A1 compares solutions obtained from our standard method of integrating Lagrange’s equations with those from REBOUND.

APPENDIX B: DIFFERENCES BETWEEN THE 3:2 AND 2:1 RESONANCES

In Section 3, we showed that TTV phases are easier to excite for planets near the 2:1 than near the 3:2 (Fig. 9). How the degree of phase excitation depends on mass ratio m₁/m₂ also differs between the two resonances (Fig. 10). Here, we explore the reasons for these behaviours.

As discussed in Section 2.2, a non-zero TTV phase shift depends on a time-varying mean motion. After disc torques are removed, and the third body enforces ϖ₁ = ϖ₂ so that ϕ₁ = ϕ₂, the only way to generate mean-motion variations in planets 1 and 2 is for their eccentricities to deviate from the ratio (e₂/e₁)_res = |f₁/f₂| [see equations (A1) and (A10)]. Deviations from (e₂/e₁)_res are easier to achieve at larger Δ, i.e. for systems that behave more secularly than resonantly. In the secular limit, planets 1 and 2 are driven by disc eccentricity damping into a secular eigenmode – a fixed point in the secular theory, where the eccentricities do not change and the apses precess at the same rate (for given α₁₂), distinct from the resonant fixed point discussed elsewhere in this paper. The relevant mode has ϖ₁ = ϖ₂ and is easiest to write down in the test particle limit: When m₁ ≪ m₂, (e₂/e₁)_sec = |2C_A/C_B|, and when m₂ ≪ m₁, (e₂/e₁)_sec = |C_B/2C_A|.

Fig. B1 confirms that post-disc damping, e₂/e₁ varies between the resonance-dominated limit (e₂/e₁)_res at small Δ, and the secular-dominated limit (e₂/e₁)_sec at larger Δ. The figure derives from simulations like the one in Fig. 5, but performed in the m₁ ≪ m₂ and m₂ ≪ m₁ test particle limits, and with varying t_e/|t_a| to generate different equilibrium Δ values. The values of e₂/e₁ plotted are measured post-disc, at the peak of the secular cycle driven by planet 3; we found that this ratio hardly changed throughout the cycle.

We see from Fig. B1 that 3:2 pairs are confined to a much narrower range of e₂/e₁ than 2:1 pairs. The restriction is particularly severe for 3:2 resonance in the m₂ ≪ m₁ limit (bottom right panel). The smaller the range of permitted values, the closer e₂/e₁ is to (e₂/e₁)_res, and the harder it becomes to excite TTV phases. Thus, we can make some sense of the trends seen in Figs 9 and 10. The 2:1 resonance is more susceptible to secular forcing than the 3:2, in part because the 2:1 is weakened (|f₂| is made smaller) by the indirect potential.