Evidence for a Nondichotomous Solution to the Kepler Dichotomy: Mutual Inclinations of Kepler Planetary Systems from Transit Duration Variations

Sarah C. Millholland; Matthias Y. He; Eric B. Ford; Darin Ragozzine; Daniel Fabrycky; Joshua N. Winn

doi:10.3847/1538-3881/ac0f7a

1. Introduction

One of the fundamental features of our solar system is the near coplanarity of the planetary orbits. The orbital inclinations of the solar system planets relative to the invariable plane are a few degrees or less, with the exception of Mercury, which is inclined by 6°. This coplanarity was one of the key observables leading to the nebular hypothesis (Kant 1755; Laplace 1796) and the modern-day paradigm of planet formation in thin gaseous disks (e.g., Armitage 2011). It is unclear whether most planetary systems retain their primordially small inclinations (like the solar system) or develop larger misalignments. Accordingly, constraining the mutual inclination distribution of exoplanetary systems is a large-scale objective with fundamental implications for planet formation theory (e.g., Hansen & Murray 2013) and population statistics (e.g., Tremaine & Dong 2012).

Although the transit and Doppler methods do not generally reveal the mutual inclinations of multiplanet systems without further assumptions, there have still been several advancements in our knowledge of the mutual inclinations of the population of short-period (P ≲ 1 yr), tightly packed systems discovered by the Kepler mission (Borucki et al. 2010; Lissauer et al. 2011; Batalha et al. 2013; Fabrycky et al. 2014). (See Zhu & Dong 2021 for a recent review.) Roughly half of Sun-like stars host these short-period planets (Winn & Fabrycky 2015; Zhu et al. 2018; He et al. 2019; Yang et al. 2020), most of which are likely in multiplanet (but not necessarily multiple-transiting) systems (Mulders et al. 2018; He et al. 2020). Such systems will be our primary focus in this paper.

Kepler systems of multiple transiting planets ("Kepler multis") are nearly coplanar on average (e.g., Lissauer et al. 2011).¹⁰ This conclusion has been reached by analyzing, for instance, the distribution of ratios of transit chord lengths between pairs of planets within the same system (Steffen et al. 2010; Fang & Margot 2012; Fabrycky et al. 2014). According to these studies, at least half of Kepler systems have mutual inclinations that are consistent with a Rayleigh distribution with scale parameter σ_i ∼ 1°–2°.

Within the Kepler sample, there is also evidence for a subset of systems with larger mutual inclinations, based on analyses of the observed transiting multiplicity distribution. Inferences from this distribution are difficult, though, because it depends on both the intrinsic multiplicity distribution and the mutual inclination distribution (e.g., Lissauer et al. 2011; Tremaine & Dong 2012). To break the degeneracy, additional assumptions or observations are required, such as inputs from radial velocity (RV) surveys or the incidence of transit timing variations (TTVs).

Several studies used statistical planet population models to fit the observed transiting multiplicity distribution of Kepler systems (e.g., Lissauer et al. 2011; Fang & Margot 2012; Tremaine & Dong 2012; Ballard & Johnson 2016). These studies identified an apparent excess of single-transiting systems ("singles"), compared to the number of singles one would expect if all planetary systems have the same ∼1°–2° mutual inclination dispersion (and intrinsic multiplicity distribution) of the transit multis. This apparent discrepancy between the observed and expected number of singles is known as the "Kepler dichotomy," and it is sometimes interpreted as evidence for two (or more) planet populations with distinct architectures (e.g., Lissauer et al. 2011; Johansen et al. 2012; Hansen & Murray 2013; Ballard & Johnson 2016). Biases due to detection order in multiple-transiting systems likely contribute to the Kepler dichotomy (Zink et al. 2019) but cannot entirely explain it (He et al. 2019).

One proposed resolution to the Kepler dichotomy is to invoke a substantial fraction of intrinsically single-planet systems (Fang & Margot 2012; Sandford et al. 2019). While such a model can fit the observed multiplicity distribution, He et al. (2019) showed that it is an unlikely explanation because it requires too many stars to have only one planet. Specifically, after assigning nearly coplanar multiplanet systems to target stars with an occurrence derived from the transiting multiplanet systems, the remaining proportion of stars is not large enough to account for the observed number of transit singles. Moreover, the observed incidence of TTVs of Kepler planets does not provide evidence for a large population of intrinsic singles (Ford et al. 2011; Holczer et al. 2016).

The leading explanation of the Kepler dichotomy is that there is a distribution of mutual inclinations that includes a substantial fraction of systems with low mutual inclinations and a separate fraction of systems with significantly larger mutual inclinations. For example, there might be two (or more) subpopulations, one with 1°–2° mutual inclinations and one with >30° mutual inclinations. An alternative is a smooth distribution, such as a model in which the inclinations are drawn from a Rayleigh distribution with a width parameter that is itself a variable drawn from a smooth distribution (e.g., Lissauer et al. 2011) or that depends on the number of planets in the system (e.g., Moriarty & Ballard 2016; Zhu et al. 2018). However, both the statistical properties and physical origins of these variations are uncertain.

Recently, He et al. (2019, 2020) built a forward modeling framework capable of fitting a statistical description of the underlying planet population to the survey data by "observing" simulated planetary systems with the Kepler detection pipeline. He et al. (2019, hereafter H19) showed that the observed multiplicity distribution (together with several other aspects of the Kepler population, including the distributions of orbital periods, period ratios, transit depths, depth ratios, and transit durations) can be described by two populations consisting of a low and high mutual inclination component, with σ_i,low ≈1°–2°, σ_i,high ≈ 30°–65°, and ∼40% of systems having high mutual inclinations. H19's findings are similar to those of Mulders et al. (2018), who used an analogous forward modeling framework and also found the data to be compatible with a dichotomous two-population model.

A more recent model by He et al. (2020, hereafter H20) showed that the observations can also be reproduced when the mutual inclinations (and eccentricities) are distributed according to a stability limit dictated by the system's total angular momentum deficit (AMD; Laskar & Petit 2017). An interesting prediction of this model is an inverse correlation between the mutual inclination dispersion and the intrinsic multiplicity that is well described by a power law, σ_i,n ∝ n^α (with n the multiplicity and α < 0). This feature is qualitatively similar to the model of Zhu et al. (2018), and both lead to the same key conclusion that the mutual inclination distribution does not necessarily have to be dichotomous, but rather it can be characterized by a broad and multiplicity-dependent distribution.

The H19 and H20 models provide comparably good fits to the observations, which again reflects the fundamental degeneracy between the intrinsic multiplicity and mutual inclination distributions. Breaking this degeneracy requires an extra source of information, such as data from RVs (Figueira et al. 2012; Tremaine & Dong 2012) or TTVs (Zhu et al. 2018). In this work, we consider transit duration variations (TDVs) as a hitherto-unexploited source of extra information that is highly sensitive to the mutual inclination distribution. Secular (long-term average) planet–planet interactions lead to apsidal and nodal precession of the orbits on a timescale of 10²–10³ yr for typical Kepler planets (Millholland & Laughlin 2019). As first pointed out by Miralda-Escudé (2002), this orbital precession leads to variations in the transit duration of transiting planets that manifest as a slow drift on observable timescales (see Figure 1). The drift timescale is sensitive to mutual inclinations because the signal goes as ${\dot{T}}_{\mathrm{dur}}\propto \dot{{\rm{\Omega }}}P\sin i$ (Miralda-Escudé 2002), where i is the transiting planet's inclination with respect to the invariable plane and Ω is the longitude of the ascending node. Moreover, TDVs of planets in single-transiting systems encode information about inclined, nontransiting companion planets.

**Figure 1.** Schematic illustration of the physical setup and concept of this study. A transiting planet experiences orbital precession induced by secular perturbations from a companion planet, which is nontransiting in this illustration. The orbital precession leads to changes in the transiting planet's inclination with respect to the line of sight (but not to the invariable plane) and the transit duration, T_dur, here defined as the duration during which the center of the planet is projected in front of the stellar disk. Over long timescales (10²–10³ yr), these TDVs manifest as an oscillation at the orbital precession period, as indicated with the example time evolution plot on the bottom left. However, over short timescales, such as the ∼4 yr baseline of the Kepler prime mission, the TDVs manifest as a linear trend with slope ${\dot{T}}_{\mathrm{dur}}$ . This is indicated with the zoom-in on the bottom right. The TDV slope depends on the mutual orbital inclinations. Examining the TDVs of an ensemble of planets thus allows us to place constraints on the distribution of mutual inclinations.
Download figure:
Standard image High-resolution image

**Figure 1.** Schematic illustration of the physical setup and concept of this study. A transiting planet experiences orbital precession induced by secular perturbations from a companion planet, which is nontransiting in this illustration. The orbital precession leads to changes in the transiting planet's inclination with respect to the line of sight (but not to the invariable plane) and the transit duration, T_dur, here defined as the duration during which the center of the planet is projected in front of the stellar disk. Over long timescales (10²–10³ yr), these TDVs manifest as an oscillation at the orbital precession period, as indicated with the example time evolution plot on the bottom left. However, over short timescales, such as the ∼4 yr baseline of the Kepler prime mission, the TDVs manifest as a linear trend with slope ${\dot{T}}_{\mathrm{dur}}$ . This is indicated with the zoom-in on the bottom right. The TDV slope depends on the mutual orbital inclinations. Examining the TDVs of an ensemble of planets thus allows us to place constraints on the distribution of mutual inclinations.
Download figure:
Standard image High-resolution image

TDV signals have been detected for ∼30 Kepler planets (Holczer et al. 2016; Kane et al. 2019; Shahaf et al. 2021). Notable examples include Kepler-108 (Mills & Fabrycky 2017), Kepler-693 (Masuda 2017), and Kepler-9 (Freudenthal et al. 2018; Borsato et al. 2019). The TDVs have led to mutual inclination constraints in these systems. In other examples (e.g., Sanchis-Ojeda et al. 2012), an absence of TDVs has been used as evidence for low mutual inclinations. Recently, Boley et al. (2020) investigated observed planets that are the best candidates for exhibiting detectable TDV signals in near-future observations. In addition, Dawson (2020) described how to use TDVs to infer systems' three-dimensional architectures. While this paper was in review, Shahaf et al. (2021) published a new comprehensive analysis of TDVs of Kepler planets, building off of Holczer et al. (2016). Overall, their results are complementary to this work; where relevant, we will discuss specific comparisons.

Here we use TDVs to constrain the mutual inclination distribution of Kepler systems. We accomplish this by comparing the TDV statistics of the observed planet population to expectations from simulated planet populations constructed by the forward models of H19 and H20. Both models reproduce many aspects of the Kepler survey statistics, but they have significantly different intrinsic mutual inclination distributions. In Section 2, we describe the relevant equations for orbital-precession-induced TDVs, including an analytic calculation of ${\dot{T}}_{\mathrm{dur}}$ . In Section 3, we describe our methods for calculating the TDVs of the simulated planets and comparing their properties to the observed Kepler planets with TDVs. We summarize our results in Section 4. Namely, we show that the TDV statistics of the observed planets strongly support the nondichotomous model of H20 and disfavor the dichotomous model of H19. In Section 5, we review the proposed theories for generating mutual inclinations among Kepler planets and discuss which are most consistent with our results.

2. Analytic Calculation of TDVs

We begin with an analytic calculation of a planet's TDV, which is related to the time derivative of its transit duration, T_dur. Our derivation bears resemblance to that of Miralda-Escudé (2002), but we allow for arbitrary eccentricities rather than assuming circular orbits. We will consider the orbital evolution to be secular, with the semimajor axis, a, approximately constant. In this case, the time evolution of T_dur is driven by orbital precession of the longitude of the ascending node, Ω, and the argument of periapse, ω, as well as changes in the orbital eccentricity, e. Our goal is to relate ${\dot{T}}_{\mathrm{dur}}$ to $\dot{{\rm{\Omega }}}$ , $\dot{\omega }$ , and $\dot{e}$ .

Within the following derivation, there are three distinct planes to keep in mind: (1) the orbital plane, which is perpendicular to the planet's orbital angular momentum vector; (2) the invariable plane, which is perpendicular to the system's total angular momentum vector; and (3) the sky plane, which is perpendicular to the line of sight. Specifically, we define the invariable plane to be xy, with the x-axis the intersection of the invariable plane and the sky plane. The planet's orbital plane has an inclination, i, with respect to the invariable plane, and the line of nodes of the orbit forms an angle Ω with the x-axis. Note that here we will assume that i is constant. This is strictly only valid for a two-planet system, but it is a good approximation when considering TDVs over a baseline as short as the ∼4 yr Kepler mission, which is much shorter than the secular timescales over which i varies. For simplicity, we will also assume R_p ≪ R_⋆, a good approximation for the sub-Neptune-sized planets we are focusing on here. It is straightforward to generalize the treatment without this assumption (Miralda-Escudé 2002).

For a ≫ R_⋆, a planet's transit duration is given by

$\begin{eqnarray}&&{T}_{\mathrm{dur}}=\displaystyle \frac{2{R}_{\star }\sqrt{1-{b}^{2}}}{{v}_{\mathrm{mid}}},\end{eqnarray} \tag{ 1 }$

where R_⋆ is the stellar radius, b is the dimensionless impact parameter, and v_mid is the sky-projected orbital velocity at midtransit, given by

$\begin{eqnarray}&&{v}_{\mathrm{mid}}=\displaystyle \frac{{na}(1+e\sin {\omega }_{\mathrm{sky}})}{\sqrt{1-{e}^{2}}}.\end{eqnarray} \tag{ 2 }$

Here, n = 2π/P is the mean motion, e is the orbital eccentricity, and ω_sky is the argument of periapse measured with respect to the sky plane. The argument ω_sky is related to the corresponding invariable plane angle ω via the expression (Judkovsky et al. 2020)

$\begin{eqnarray}&&\tan ({\omega }_{\mathrm{sky}}-\omega )=\displaystyle \frac{\sin {\rm{\Omega }}\cos \beta }{\cos {\rm{\Omega }}\cos i\cos \beta +\sin i\sin \beta }.\end{eqnarray} \tag{ 3 }$

Here, β is the fixed angle between the invariable plane and the line of sight. Put another way, β = π/2 − i_inv, where i_inv is the angle between the invariable plane and the sky plane.

Taking the time derivative of Equation (1) yields

$\begin{eqnarray}&&{\dot{T}}_{\mathrm{dur}}=-{T}_{\mathrm{dur}}\left[\dot{b}\displaystyle \frac{b}{(1-{b}^{2})}+\displaystyle \frac{{\dot{v}}_{\mathrm{mid}}}{{v}_{\mathrm{mid}}}\right].\end{eqnarray} \tag{ 4 }$

The time derivatives $\dot{b}$ and ${\dot{v}}_{\mathrm{mid}}$ are driven by orbital precession of the longitude of the ascending node and the argument of periapse, as well as changes in the orbital eccentricity. These are parameterized by $\dot{\omega }$ , $\dot{{\rm{\Omega }}}$ , and $\dot{e}$ , respectively. As for the first term in Equation (4), we can compute $\dot{b}$ by using the definition of the dimensionless impact parameter,

$\begin{eqnarray}&&b=\displaystyle \frac{{r}_{\mathrm{mid}}\sin \alpha }{{R}_{\star }},\end{eqnarray} \tag{ 5 }$

where r_mid is the star−planet separation at midtransit,

$\begin{eqnarray}&&{r}_{\mathrm{mid}}=\displaystyle \frac{a(1-{e}^{2})}{1+e\sin {\omega }_{\mathrm{sky}}}.\end{eqnarray} \tag{ 6 }$

The quantity α is the angle between the orbital plane of the planet and the line of sight. (Alternatively, α = π/2 − i_sky, where i_sky is the orbital inclination with respect to the sky plane.) This angle is related to the other angles in the problem via (Judkovsky et al. 2020)

$\begin{eqnarray}&&\sin \alpha =-\sin i\cos {\rm{\Omega }}\cos \beta +\cos i\sin \beta .\end{eqnarray} \tag{ 7 }$

Using these relationships, we obtain

$\begin{eqnarray}&&\dot{b}=\displaystyle \frac{1}{{R}_{\star }}\left[{\dot{r}}_{\mathrm{mid}}\sin \alpha +{r}_{\mathrm{mid}}\dot{{\rm{\Omega }}}\sin i\cos \beta \sin {\rm{\Omega }}\right],\end{eqnarray} \tag{ 8 }$

where ${\dot{r}}_{\mathrm{mid}}$ can be calculated from Equation (6),

$\begin{eqnarray}&&{\dot{r}}_{\mathrm{mid}}=\displaystyle \frac{a[-2e\dot{e}(1+e\sin {\omega }_{\mathrm{sky}})-(1-{e}^{2})\tfrac{{\rm{d}}}{{\rm{d}}t}(e\sin {\omega }_{\mathrm{sky}})]}{{\left(1+e\sin {\omega }_{\mathrm{sky}}\right)}^{2}}.\end{eqnarray} \tag{ 9 }$

As for the second term in the expression for ${\dot{T}}_{\mathrm{dur}}$ (Equation (4)), ${\dot{v}}_{\mathrm{mid}}$ can be calculated by differentiating Equation (2),

$\begin{eqnarray}&&{\dot{v}}_{\mathrm{mid}}={na}\displaystyle \frac{\dot{e}(e+\sin {\omega }_{\mathrm{sky}})+{\dot{\omega }}_{\mathrm{sky}}e\cos {\omega }_{\mathrm{sky}}(1-{e}^{2})}{{\left(1-{e}^{2}\right)}^{\tfrac{3}{2}}}.\end{eqnarray} \tag{ 10 }$

We can relate this expression to ω and $\dot{\omega }$ by using Equation (3), from which $\dot{\omega }$ can be calculated straightforwardly, recalling that both β and i are held fixed.

2.1. Lagrange's Planetary Equations

Equipped with the analytic formula for ${\dot{T}}_{\mathrm{dur}}$ in terms of $\dot{{\rm{\Omega }}}$ , $\dot{\omega }$ , and $\dot{e}$ , we now move on to the analytic calculations of the three latter quantities. Lagrange's planetary equations yield (Murray & Dermott 1999)

$\begin{eqnarray}\begin{array}{rcl}\dot{e} & = & -\displaystyle \frac{\sqrt{1-{e}^{2}}}{{{na}}^{2}e}\displaystyle \frac{\partial { \mathcal R }}{\partial \varpi }\\ \dot{{\rm{\Omega }}} & = & \displaystyle \frac{1}{{{na}}^{2}\sqrt{1-{e}^{2}}\sin i}\displaystyle \frac{\partial { \mathcal R }}{\partial i}\\ \dot{\varpi } & = & \dot{{\rm{\Omega }}}+\dot{\omega }=\displaystyle \frac{\sqrt{1-{e}^{2}}}{{{na}}^{2}e}\displaystyle \frac{\partial { \mathcal R }}{\partial e}+\displaystyle \frac{\tan \tfrac{1}{2}i}{{{na}}^{2}\sqrt{1-{e}^{2}}}\displaystyle \frac{\partial { \mathcal R }}{\partial i}.\end{array}\end{eqnarray} \tag{ 11 }$

Here, ${ \mathcal R }$ is the disturbing function, or the non-Keplerian perturbing gravitational potential experienced by the planets owing to their mutual interactions. We adopt a secular expansion of ${ \mathcal R }$ to fourth order in e and $s\equiv \sin (i/2)$ using the Appendix tables from Murray & Dermott (1999). The resulting disturbing function expansion is provided in Appendix A.

2.2. Comparison of Analytic TDV to N-body

In Appendix B, we assess the accuracy of the analytic calculation of ${\dot{T}}_{\mathrm{dur}}$ by showing how it compares to a direct N-body computation. We demonstrate that the analytic calculation performs well when inclinations are low, i ≲ 10°. In this regime, the analytic ${\dot{T}}_{\mathrm{dur}}$ is within 50% of the N-body ${\dot{T}}_{\mathrm{dur}}$ in roughly three-quarters of cases. However, when inclinations are larger than ∼10°, there can be orders-of-magnitude discrepancies between the analytic calculation and the N-body result. This occurs because the disturbing function expansion breaks down for large e and i, leading to perturbations in the equations for $\dot{e}$ , $\dot{{\rm{\Omega }}}$ , and $\dot{\varpi }$ . As we will describe in Section 3.3.2, this accuracy problem ultimately limits the applicability of the analytic ${\dot{T}}_{\mathrm{dur}}$ calculation.

3. Methods

Our goal is to constrain the underlying distribution of Kepler planet mutual inclinations by comparing the TDVs of the observed Kepler planets to those of various simulated planet populations with different mutual inclination distributions. In this section, we describe our methods. We define our observed and simulated planet populations in Sections 3.1 and 3.2, respectively. We then explain the procedure for identifying detectable TDVs of the simulated planets in Section 3.3.

3.1. Observations: TDVs of Kepler Planets

We first identify a set of Kepler planets with significant TDVs. These planets will later be compared to the simulated planets. As part of their comprehensive catalog of TTV measurements for 2599 Kepler Objects of Interest (KOIs) across the full 17 quarters of the Kepler mission, Holczer et al. (2016, hereafter H16) also measured transit depths and durations. The duration and depth measurements were limited to cases with T_dur > 1.5 hr and S/N > 10, where S/N is the signal-to-noise ratio per transit. As a result, a total of 779 KOIs have duration and depth measurements. These were then analyzed for any TDVs or transit depth variations by identifying potentially significant periodicities or long-term trends. The long-term trends, quantified by the TDV slope ${\dot{T}}_{\mathrm{dur}}$ , are the most relevant for our analysis, as these are the expected signal of duration drifts induced by orbital precession.

H16 summarized their results for TDVs in the TDV_statistics.txt file at ftp://wise-ftp.tau.ac.il/pub/tauttv/TTV/ver_112. We will refer to this as the "H16 TDV catalog." In particular, the quantities that are relevant to our analysis are the slope of the TDV measurements and the estimated error of the slope. We apply several cuts to the H16 TDV catalog to establish a list of observed planets with detected TDVs.

1.
Significant TDV slope: We consider a detectable transit duration drift to be one with $| {\dot{T}}_{\mathrm{dur}}| \gt 3\ {\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ . A total of 31 KOIs pass this criterion.¹¹
2.
Cuts on planet properties: The simulated planet population that we will introduce in Section 3.2 is restricted to planets with 3 days < P < 300 days and 0.5 R_⊕ < R_p < 10 R_⊕. Applying these cuts leaves 21 remaining KOIs. Most of the removed planets are cut for having R_p > 10 R_⊕. Three have P < 3 days.
3.
Cuts on stellar properties: The simulated planet population is built using a curated sample of FGK dwarf stars, developed using a series of quality cuts on the Kepler DR25 target list in conjunction with target information from Gaia DR2 (H19; H20). We require the host stars of the observed planets to be in this cleaned stellar input catalog. After this requirement, we are left with 16 remaining KOIs.

The final list of 16 KOIs with significant TDV slopes is shown in Table 1. The 16 KOIs are part of 15 systems. Of these, 7 are single-transiting systems, 4 are double-transiting systems, 3 are triple-transiting systems, and 1 is a quintuple-transiting system. Examples of TDV drift signals observed for Kepler-9 b and c are shown in Figure 2.

**Figure 2.** Examples of the TDV drift signal, shown here for Kepler-9 b and c (KOI 377.01 and 377.02). The data for the TDVs, given by $\mathrm{TDV}\,=({T}_{\mathrm{dur}}-{\bar{T}}_{\mathrm{dur}})/{\bar{T}}_{\mathrm{dur}}$ , are taken from H16. The green curves represents the best-fit lines, with slopes equal to those from the H16 TDV catalog.
Download figure:
Standard image High-resolution image

**Figure 2.** Examples of the TDV drift signal, shown here for Kepler-9 b and c (KOI 377.01 and 377.02). The data for the TDVs, given by $\mathrm{TDV}\,=({T}_{\mathrm{dur}}-{\bar{T}}_{\mathrm{dur}})/{\bar{T}}_{\mathrm{dur}}$ , are taken from H16. The green curves represents the best-fit lines, with slopes equal to those from the H16 TDV catalog.
Download figure:
Standard image High-resolution image

Table 1. KOIs with Significant TDV Slopes

KOI	Kepler Name	N_obs	P (days)	R_p (R_⊕)	${\dot{T}}_{\mathrm{dur}}$ (minutes yr⁻¹)
103.01	⋯	1	14.911	${2.55}_{-0.05}^{+0.06}$	5.0 ± 1.1
137.02	Kepler-18 d	3	14.859	${5.16}_{-0.11}^{+0.09}$	1.73 ± 0.41
139.01	Kepler-111 c	2	224.779	${6.91}_{-0.17}^{+0.16}$	8.1 ± 2.0
142.01	Kepler-88 b	1	10.916	${3.84}_{-0.31}^{+0.08}$	1.98 ± 0.58
209.02	Kepler-117 b	2	18.796	${8.36}_{-0.52}^{+0.33}$	2.98 ± 0.82
377.01	Kepler-9 b	3	19.271	${7.9}_{-0.16}^{+0.16}$	1.64 ± 0.3
377.02	Kepler-9 c	3	38.908	${8.14}_{-0.18}^{+0.18}$	−4.58 ± 0.56
460.01	Kepler-559 b	2	17.588	${3.41}_{-0.09}^{+0.11}$	6.0 ± 1.5
806.01	Kepler-30 d	3	143.206	${8.79}_{-0.31}^{+0.49}$	6.3 ± 1.9
841.02	Kepler-27 c	5	31.33	${6.24}_{-0.28}^{+0.24}$	3.6 ± 1.1
872.01	Kepler-46 b	2	33.601	${7.45}_{-0.26}^{+0.2}$	5.78 ± 0.83
1320.01	Kepler-816 b	1	10.507	${9.44}_{-0.36}^{+0.46}$	−1.88 ± 0.44
1423.01	Kepler-841 b	1	124.42	${5.02}_{-0.19}^{+0.25}$	3.7 ± 1.1
1856.01	⋯	1	46.299	${2.24}_{-0.05}^{+0.08}$	−11.8 ± 3.1
2698.01	Kepler-1316 b	1	87.972	${3.4}_{-0.08}^{+0.11}$	18.1 ± 4.5
2770.01	⋯	1	205.386	${2.26}_{-0.08}^{+0.13}$	19.4 ± 6.4

Note. List of KOIs from the H16 TDV catalog with $| {\dot{T}}_{\mathrm{dur}}| \gt 3\ {\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ that pass our additional cuts on planet and stellar properties. N_obs is the observed transit multiplicity of the system according to the Kepler DR25 catalog (Thompson et al. 2018). Orbital periods are drawn from DR25, planetary radii are from Berger et al. (2020), and the ${\dot{T}}_{\mathrm{dur}}$ measurements are from the H16 TDV catalog.

Download table as: ASCII Typeset image

3.2. Simulations: SysSim Planet Populations

We consider populations of simulated planets from the SysSim (short for "Planetary Systems Simulator") forward modeling framework. SysSim is an empirical model that generates simulated planetary systems according to flexibly parameterized statistical descriptions. It was first introduced by Hsu et al. (2018) and has undergone continuous development by Hsu et al. (2019), H19, H20, and He et al. (2021). The model has been calibrated using a range of summary statistics of the observed Kepler planet population, including (but not limited to) the observed distributions of multiplicities, orbital periods and period ratios, transit depths and depth ratios, and transit durations. SysSim is implemented in the Julia programming language (Bezanson et al. 2014) within the ExoplanetsSysSim.jl package. Both the core SysSim code (https://github.com/ExoJulia/ExoplanetsSysSim.jl) and the specific forward models explored in this study (https://github.com/ExoJulia/SysSimExClusters) are publicly available. More information can be found in the key papers mentioned above. We provide a brief description of the model below.

The process of the SysSim forward modeling framework is to first generate an underlying population of planetary systems ("physical catalog") according to a statistical model of the intrinsic distribution of systems. The physical catalog is a simulated realization designed to resemble the entire population of Kepler planetary systems, including planets that are not observed. The next step is to generate an observed population of planetary systems ("observed catalog") from the physical catalog by simulating the full Kepler detection pipeline and determining which planets would be detected by the pipeline and labeled as planet candidates during the automated vetting process. This simulated observed catalog is then compared with the true Kepler planet population using a set of summary statistics and a distance function. The preceding steps are repeated iteratively in order to optimize the distance function and identify the best-fit parameters of the statistical model. Finally, approximate Bayesian computing is applied to approximate the posterior distributions of model parameters.

In this work, we do not implement extensions of SysSim or fit new statistical models. Rather, we examine sets of simulated catalogs from two previously optimized models from H19 and H20. These models effectively describe many aspects of the Kepler population statistics, and they fit the data nearly equally well. By construction, the two models are similar in most ways, but they differ primarily in the distributions of eccentricities, inclinations, and number of planets per star.

3.2.1. Two-Rayleigh Model

H19 considered independent distributions of eccentricities and inclinations. The eccentricities were drawn from a Rayleigh distribution, e ∼ Rayleigh(σ_e). The mutual inclinations were drawn from one of two Rayleigh distributions, representing a low- and high-inclination population, with the fraction of systems belonging to the high-inclination population being ${f}_{{\sigma }_{i,\mathrm{high}}}$ .¹² That is,

$\begin{eqnarray}i\sim \left\{\begin{array}{ll}\mathrm{Rayleigh}({\sigma }_{i,\mathrm{high}}), & u\lt {f}_{{\sigma }_{i,\mathrm{high}}}\\ \mathrm{Rayleigh}({\sigma }_{i,\mathrm{low}}), & u\geqslant {f}_{{\sigma }_{i,\mathrm{high}}},\end{array}\right.\end{eqnarray} \tag{ 12 }$

where u ∼ Unif(0, 1). Because of the dichotomous nature of the inclination distribution, we will call this the "two-Rayleigh model."

After fitting to the Kepler data, H19 identified best-fit parameters equal to ${\sigma }_{i,\mathrm{low}}={1.40}_{-0.39}^{+0.54}$ deg, ${\sigma }_{i,\mathrm{high}}={48}_{-17}^{+17}$ deg, and ${f}_{{\sigma }_{i,\mathrm{high}}}={0.42}_{-0.07}^{+0.08}$ . The high-inclination component makes up nearly half of the total population of systems (and nearly a third of all planets), but the mutual inclination scale σ_i,high is poorly constrained, although clearly greater than ∼10°. We note that, in this work, we use the two-Rayleigh model of He et al. (2021), which is similar to H19 but with an additional dependence of the fraction of stars with planets on the host star color. This feature is also present in the model introduced in the next section, allowing for a better comparison.

3.2.2. Maximum AMD Model

The second model uses a joint distribution in which the eccentricity and inclination distributions are dependent and correlated with the number of planets in a system. This contrasts with the first model, where eccentricities and inclinations are independent of one another and of the number of planets in a system. The second approach is based on the argument that a system's long-term orbital stability is governed by its AMD. The AMD is defined as the difference between the total orbital angular momentum of the system and what it would be if all orbits were circular and coplanar (e.g., Laskar & Petit 2017). For a system of N planets, the AMD can be written as

$\begin{eqnarray}&&\mathrm{AMD}=\sum _{k=1}^{N}{{\rm{\Lambda }}}_{k}(1-\sqrt{1-{e}_{k}^{2}}\cos {i}_{k}),\end{eqnarray} \tag{ 13 }$

where ${{\rm{\Lambda }}}_{k}={M}_{p,k}\sqrt{{{GM}}_{\star }{a}_{k}}$ . This quantity is effectively a measure of the dynamical excitation of the system, and it is a conserved quantity when the orbits are evolving secularly. The AMD has been shown to be a reasonable predictor of long-term stability, at least when mean motion resonances (MMRs) are absent and when distant and dynamically detached planets are not included in the AMD calculation (Laskar & Petit 2017).

The assumption of H20 is that all systems have the critical (i.e., maximum) AMD for stability, calculated using the analytic criteria from Laskar & Petit (2017) and Petit et al. (2017). The physical justification for this assumption is that collisional events during planet formation decrease the total AMD, such that systems may generally evolve from outside the stability limit to just inside after a sequence of collisions. While this "maximum AMD model" assumption breaks down at a detailed level (i.e., not all systems are exactly at the critical AMD), it is a useful physical framework for assigning system orbital properties and has been shown to reproduce many aspects of the Kepler data (H20). In this work, we are primarily interested in the model's mutual inclination and eccentricity distributions.

Given a set of planet masses and orbital periods of a system, H20 assigned orbital properties by (1) calculating the system's critical AMD, (2) distributing the AMD among the planets per unit mass, and (3) further randomly dividing each planet's AMD into eccentricity and inclination components. This approach results in eccentricities and inclinations that are correlated with one another at the population level. Another key feature that emerges from the model is an inverse trend of inclination and eccentricity dispersion with intrinsic multiplicity. In particular, H20 found that the median mutual inclination, ${\tilde{\mu }}_{i,n}$ , of systems with n = 2, 3,..., 10 planets is well modeled by a power-law relationship of the form

$\begin{eqnarray}&&{\tilde{\mu }}_{i,n}={\tilde{\mu }}_{i,5}{\left(\displaystyle \frac{n}{5}\right)}^{\alpha },\end{eqnarray} \tag{ 14 }$

with ${\tilde{\mu }}_{i,5}={1.10}_{-0.11}^{+0.15}$ and $\alpha =-{1.73}_{-0.08}^{+0.09}$ . This power-law relationship is qualitatively similar but shallower than that inferred by Zhu et al. (2018). Equation (14) illustrates that the typical mutual inclinations of systems within the maximum AMD model are less than a few degrees, in sharp contrast with the two-Rayleigh model.

3.2.3. Model Comparison

The two-Rayleigh model and maximum AMD model both fit the Kepler data roughly equally well (see H20 Section 3.1), but they have very different underlying inclination and eccentricity distributions. To illustrate these differences, Figure 3 displays scatter plots of inclinations and eccentricities of intrinsic multiplanet systems in SysSim physical catalogs. We observe that eccentricities and inclinations are strongly correlated when distributed according to the maximum AMD model. Moreover, the maximum AMD model contains considerably fewer systems with very high inclination configurations than the two-Rayleigh model. For instance, 97% of the mutual inclinations are below 10° in the maximum AMD example plotted in Figure 3, whereas 68% and 63% of the inclinations are below 10° in the σ_i,high = 20° and σ_i,high = 60° two-Rayleigh model examples, respectively.

**Figure 3.** Scatter plots of orbital inclinations referenced to the invariable plane (i) vs. eccentricity (e) for intrinsic multiplanet systems in SysSim physical catalogs. The top panel considers a single physical catalog from the maximum AMD model. The middle and bottom panels consider physical catalogs from the two-Rayleigh model using σ_i,high = 20° and σ_i,high = 60°, respectively. To illustrate the scatter plot density, we display a set of contours calculated using a Gaussian kernel density estimation.
Download figure:
Standard image High-resolution image

As a caveat, we note that the simulated systems have not been evaluated for long-term orbital stability in either model. The systems in the maximum AMD model obey the AMD stability criterion by construction, but that does not guarantee that they are all long-term stable (Cranmer et al. 2021). They are, however, more likely to be stable than the systems in the two-Rayleigh model, which frequently have very high inclinations (sometimes even >90°). This difference means that a direct comparison of the two models is slightly misleading. A stability analysis would require significant computation but would be interesting to address in future SysSim models.

3.3. TDV Calculations of SysSim Planets

For the two-Rayleigh model and the maximum AMD model, we aim to compare the TDVs of the simulated planets with those of the observed Kepler planets. We thus need to simulate the TDV detection process for the SysSim planets in a similar way to the true observations, which are obtained from the H16 TDV catalog (Section 3.1). Clearly, we can only measure TDVs for planets in the observed catalog (i.e., that have been "detected" in transit), but we also need to consider whether they have parameters suitable for individual transit duration measurements and whether their TDV signal has a detectable slope. We will discuss this process in three steps: TDV measurability (Section 3.3.1), TDV calculation (Section 3.3.2), and TDV slope detectability (Section 3.3.3). We will then summarize the process end to end in Section 3.4.

3.3.1. TDV Measurability

We must first determine which planets in a given observed catalog have transits with sufficiently high S/N such that individual transit durations would be measurable and the planet would enter into the H16 TDV catalog. We summarize this as a planet having "TDV measurability." H16 derived individual transit durations only when T_dur > 1.5 hr and S/N > 10; we thus apply these same thresholds to the SysSim planets.

We calculate the individual transit S/N of planets in the SysSim observed catalogs as

$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{\delta }{{\mathrm{CDPP}}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 15 }$

Here δ is the transit depth and CDPP_eff is the effective combined differential photometric precision (CDPP), given by

$\begin{eqnarray}&&{\mathrm{CDPP}}_{\mathrm{eff}}={f}_{\mathrm{CDPP}}\times {\mathrm{CDPP}}_{1.5\ \mathrm{hr}}\sqrt{\displaystyle \frac{1.5\ \mathrm{hr}}{{T}_{\mathrm{dur}}}}.\end{eqnarray} \tag{ 16 }$

CDPP_{1.5 hr} is the mission average of the 1.5 hr duration CDPP for the target star, taken from the Kepler Q1–Q17 DR25 stellar catalog (Thompson et al. 2018). Equation (16) also contains a scaling factor, f_CDPP < 1, that accounts for a systematic offset in the calculated S/N compared to that of H16. The systematic offset arises because H16 derived the transit S/N using measurement uncertainties, rather than the CDPP. The CDPP is larger because it includes both photon noise and variability due to the star and the instrument. We identify the optimal scaling factor by calculating the S/N for the planets in the H16 TDV catalog and minimizing the difference with respect to H16's reported S/N values. The resulting value is f_CDPP = 0.8707. Figure 4 shows the calculated S/N versus the H16 S/N, before and after the scaling factor correction.

**Figure 4.** Calibration of the S/N calculation. We show the calculated S/N (using Equation (15)) vs. the H16 S/N for the planets in the H16 TDV catalog. The top panel is before the scaling correction (f_CDPP = 1), and the bottom panel is after the correction (f_CDPP = 0.8707). The black line is one-to-one.
Download figure:
Standard image High-resolution image

3.3.2. TDV Calculation

After a planet is deemed to have sufficiently high transit S/N such that the individual transit durations would be measurable, we obtain an estimate of its TDV slope, ${\dot{T}}_{\mathrm{dur}}$ . All planets in the system (transiting or not) are accounted for in this calculation. For systems with two or more planets, we calculate ${\dot{T}}_{\mathrm{dur}}$ in two ways. The first is the analytic calculation outlined in Section 2. The second is a numerical approach aided by the REBOUND N-body gravitational dynamics software (Rein & Liu 2012). We use the REBOUND Wisdom–Holman integrator (Rein & Tamayo 2015) to calculate a short orbital evolution over a 4 yr baseline (approximately the length of the Kepler prime mission). The integration uses a time step equal to 0.1P₁. We do not account for general relativity (or any other additional forces), since the perturbations on ${\dot{T}}_{\mathrm{dur}}$ from general relativistic apsidal precession are negligible relative to the influences of nodal precession. Following the integration, we calculate T_dur as a function of time using Equation (1) (along with the other relevant relationships for b and v_mid). Finally, we perform a least-squares linear fit to T_dur versus time to calculate ${\dot{T}}_{\mathrm{dur}}$ .

As discussed in Section 2.2 and Appendix B, the analytic calculation of ${\dot{T}}_{\mathrm{dur}}$ is adequate (relative to N-body) when inclinations are low (i.e., i ≲ 10°, as in the maximum AMD model) but can have orders-of-magnitude discrepancies at high inclinations (i.e., in the two-Rayleigh model). In order to avoid systematic errors in our study of the two models, we opt to use the N-body ${\dot{T}}_{\mathrm{dur}}$ calculation for all analyses going forward. The analytic derivation, however, is still useful for physical insight and for quick calculations in low-inclination systems.

3.3.3. TDV Slope Detectability

The next step is to determine whether the calculated TDV slope, ${\dot{T}}_{\mathrm{dur}}$ , would be detectable according to the $| {\dot{T}}_{\mathrm{dur}}| \gt 3\ {\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ threshold identified above for the observed systems. We thus require an estimate of the TDV slope uncertainty, ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ , to use in conjunction with the calculated ${\dot{T}}_{\mathrm{dur}}$ . Using generalized least squares, ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ can be related to the individual transit duration uncertainty, ${({\sigma }_{{T}_{\mathrm{dur}}})}_{\mathrm{ind}}$ , via the expression

$\begin{eqnarray}&&{\sigma }_{{\dot{T}}_{\mathrm{dur}}}^{2}\approx \displaystyle \frac{{({\sigma }_{{T}_{\mathrm{dur}}}^{2})}_{\mathrm{ind}}}{{P}^{2}{f}_{0}\sum _{j=-M}^{M}{j}^{2}}.\end{eqnarray} \tag{ 17 }$

Here f₀ is the duty cycle (the fraction of data cadences with valid data), and the transits are assumed to run from − M to M in the absence of data gaps. (In principle, one could estimate the increased uncertainty due to some transits not being observed by not including the corresponding j² terms in the denominator of Equation (17).) We approximate M ≈ t_obs/(2P), where t_obs is the total observation time span. Next, we can simplify Equation (17) by expressing ${({\sigma }_{{T}_{\mathrm{dur}}})}_{\mathrm{ind}}$ in terms of the composite transit duration uncertainty, ${\sigma }_{{T}_{\mathrm{dur}}}$ , and the number of transits, ${N}_{\mathrm{tr}}$ , yielding ${({\sigma }_{{T}_{\mathrm{dur}}}^{2})}_{\mathrm{ind}}={\sigma }_{{T}_{\mathrm{dur}}}^{2}{N}_{\mathrm{tr}}={\sigma }_{{T}_{\mathrm{dur}}}^{2}{t}_{\mathrm{obs}}{f}_{0}/P$ . Substituting this into Equation (17) and adding a scaling factor, ${f}_{{\sigma }_{{\dot{T}}_{\mathrm{dur}}}}$ , yields

$\begin{eqnarray}&&{\sigma }_{{\dot{T}}_{\mathrm{dur}}}^{2}=\displaystyle \frac{{\sigma }_{{T}_{\mathrm{dur}}}^{2}{t}_{\mathrm{obs}}{f}_{{\sigma }_{{\dot{T}}_{\mathrm{dur}}}}^{2}}{{P}^{3}\sum _{j=-M}^{M}{j}^{2}}.\end{eqnarray} \tag{ 18 }$

Similar to the calculation of CDPP_eff in Equation (16), the scaling factor is required to calibrate ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ to the H16 TDV catalog. Using Equation (18) without a scaling factor leads to a systematic overestimation of the calculated ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ compared to the H16 values because of differences in the ${\sigma }_{{T}_{\mathrm{dur}}}$ calculation. It appears that H16's use of measurement uncertainties rather than the CDPP led them to underestimate ${\sigma }_{{T}_{\mathrm{dur}}}$ and that this effect was greater than the impact of Kepler's duty cycle being less than 100%. We solve for the optimal value of ${f}_{{\sigma }_{{\dot{T}}_{\mathrm{dur}}}}$ by calculating ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ for the planets in the H16 TDV catalog and minimizing the difference with respect to H16's reported values. The resulting value is ${f}_{{\sigma }_{{\dot{T}}_{\mathrm{dur}}}}=0.7378$ .

We use Equation (18) to estimate ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ for the SysSim planets with calculated TDVs. Finally, we use the $| {\dot{T}}_{\mathrm{dur}}| \gt 3\ {\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ criterion to determine whether a planet qualifies as having a detectable TDV signal.

3.4. Summary of the Full Process

We synthesize the preceding series of steps and calculate TDVs for simulated planet populations in both the two-Rayleigh model and the maximum AMD model. We consider a set of 100 physical/observed catalog pairs for each model with the parameters of the statistical models sampled according to their posterior distributions (see H20; He et al. 2021). The physical catalogs contain orbital elements (inclinations, longitudes of ascending node, etc.) referenced to the sky plane. We transform these orbital elements to be referenced to the invariable plane using Equations (3) and (7).

The observed catalogs are first used to specify which planets to perform the TDV calculation for based on whether their TDVs are measurable (Section 3.3.1). The physical catalogs determine the details of a given planet's TDV calculation (Section 3.3.2). Finally, transit properties from the observed catalog then help quantify whether the resulting TDV signal is detectable (Section 3.3.3). Each physical/observed catalog pair yields a subset of the planet population with detectable TDV slopes.

4. Results

4.1. Number of Planets with Detected TDVs

The most direct way of comparatively evaluating the two-Rayleigh model and maximum AMD model is to examine each model's total number of simulated planets with detected TDV signals. This quantity can then be compared to the number of observed planets with detected TDV signals, thereby determining which model is a better fit in terms of TDV statistics. In this section, we consider this simple tabulation; in the next section, we look deeper into the properties of the simulated and observed planets with detected TDVs.

Figure 5 presents the tabulation of TDV detections. We show histograms of the number of planets with detected TDVs for 100 physical and observed catalog pairs of each model. (That is, each physical/observed catalog pair corresponds to a single number of planets with detected TDVs.) The medians of the distributions and confidence intervals representing the 16th and 84th percentiles are shown in the figure; these values are ${43}_{-13}^{+18}$ planets with detected TDVs for the two-Rayleigh model and ${22}_{-6}^{+10}$ for the maximum AMD model. The observed number of planets with detected TDVs according to Kepler observations (16; see Section 3.1) is also shown. The maximum AMD model yields a quantity of planets with detected TDVs that is in agreement with the observations. In contrast, the two-Rayleigh model yields too many planets with detected TDVs to be compatible with the observations. Thus, the TDV statistics support the maximum AMD model (and its associated mutual inclination distribution) over the two-Rayleigh model. We will revisit this conclusion in the Discussion (Section 5).

Although the majority of outcomes for the two-Rayleigh model lead to more simulated systems with detectable TDVs than are observed, the low end of the distribution (∼20 planets with detected TDVs) is close to the observed number. This suggests that it may be difficult to rule out the two-Rayleigh model entirely. It is helpful to better understand these cases within the context of the broader distribution. To do this, we can study the relationships between the number of planets with detected TDVs and the parameters of the low- and high-inclination population components. These include the fraction of systems belonging to the high-inclination population, ${f}_{{\sigma }_{i,\mathrm{high}}}$ , and the Rayleigh distribution scale parameters of the low and high components, σ_i,low and σ_i,high.

Figure 6 shows scatter plots of these relationships. We plot the number of planets with detected TDVs versus σ_i,low, σ_i,high, and ${f}_{{\sigma }_{i,\mathrm{high}}}$ . The number of planets with detected TDVs is positively (albeit weakly) correlated with σ_i,low and ${f}_{{\sigma }_{i,\mathrm{high}}}$ , while it is negatively correlated with σ_i,high. These relationships can be understood by noting that $| \dot{{\rm{\Omega }}}| \sim \cos i$ , such that $| {\dot{T}}_{\mathrm{dur}}| \sim \sin i\cos i\sim \sin 2i$ (Equations (4), (8), and (11); see also Figure 7). When i is small, $| {\dot{T}}_{\mathrm{dur}}| \sim i$ , creating the positive correlation between the number of planets with detected TDVs and σ_i,low (left panel of Figure 6). Since $| {\dot{T}}_{\mathrm{dur}}| \sim \sin 2i$ peaks at 45°, the high-inclination component of the two-Rayleigh model has larger average TDV signals, which results in the positive trend between the number of TDV detections and ${f}_{{\sigma }_{i,\mathrm{high}}}$ (right panel of Figure 6). However, we might also expect that the number of TDV detections would peak at σ_i,high ∼ 45° rather than show a negative correlation (middle panel of Figure 6). The negative correlation arises because ${f}_{{\sigma }_{i,\mathrm{high}}}$ is inversely correlated with σ_i,high, such that a low σ_i,high is associated with many more systems in the high-inclination population, and these systems show more detectable TDVs on average.

**Figure 6.** Scatter plots of the correlations between the number of planets with detected TDVs in the two-Rayleigh model and the parameters of the model. The x-axes of the left and middle panels are the Rayleigh distribution scale parameters of the low- and high-inclination components, σ_i,low and σ_i,high. The x-axis of the right panel is the fraction of systems belonging to the high-inclination population, ${f}_{{\sigma }_{i,\mathrm{high}}}$ . To guide the eye, we include linear regression lines for each scatter plot. The horizontal dashed line corresponds to the number of planets with detected TDVs in the Kepler observations.
Download figure:
Standard image High-resolution image

**Figure 7.** Absolute value of the TDV slope, $| {\dot{T}}_{\mathrm{dur}}|$ , as a function of the inclination of the planet's orbit relative to the invariable plane. The gray points indicate all TDV calculations (without regard to TDV detection) within a single observed catalog for the maximum AMD model (top panel) and two-Rayleigh model (bottom panel). The yellow points represent the mean and standard deviation of the gray points within log-uniform bins of i. The dashed black lines represent $| {\dot{T}}_{\mathrm{dur}}| \sim \sin (2i)$ with a normalization such that the curve approximately passes through the data. While $| {\dot{T}}_{\mathrm{dur}}| \sim \sin (2i)$ is a good rough approximation, there is significant scatter due to the other factors within the ${\dot{T}}_{\mathrm{dur}}$ calculation (planet masses, semimajor axes, etc.).
Download figure:
Standard image High-resolution image

The middle panel of Figure 6 shows that in most cases where the number of planets with detected TDVs is on the low end of the distribution (≲20), the scale parameter of the high-inclination population is near its maximum, σ_i,high ≳ 50°. However, systems with extreme mutual inclinations of the level σ_i,high ≳ 50° are disfavored from a physical standpoint, since the orbits are sometimes retrograde and are more likely to be unstable owing to secular planet–planet interactions, as noted by H20. The trend shown in the middle panel of Figure 6 thus further disfavors the two-Rayleigh model; although the trend itself is weak, it is clear that obtaining a plausible number of TDV detections generally requires less plausible values of σ_i,high for stability.

4.2. Properties of Planets with Detected TDVs

In addition to the simple tabulation of the number of planets with detected TDVs, it is also valuable to compare the specific properties of the simulated planets to the corresponding observed planets with detected TDVs. Relevant properties include the magnitude of the TDV slopes, as well the planets' orbital periods and transit depths. Figures 8 and 9 (as well as Table 2, discussed later) show these properties.

**Figure 8.** Distributions of the absolute value of the TDV slope, $| {\dot{T}}_{\mathrm{dur}}|$ , for planets with detected TDVs. The top and bottom panels show the results for the two-Rayleigh model and the maximum AMD model, respectively. The solid thick line represents the median of the individual histograms for each of the 100 sets of physical and observed catalog pairs. The shaded regions represent the intervals between the 16th and 84th percentiles. The histogram of $| {\dot{T}}_{\mathrm{dur}}|$ values for the observations is shown with the dashed thick gray line.
Download figure:
Standard image High-resolution image

Table 2. Average Properties of Planets with Detected TDVs

	Average Values
	Observations	Maximum AMD Model	Two-Rayleigh Model
Planet properties:
R_p (R_⊕)	5.7 ± 0.6	${4.7}_{-0.4}^{+0.4}$	${4.8}_{-0.4}^{+0.3}$
M_p (M_⊕)	⋯	${13.1}_{-3.4}^{+10.3}$	${18.1}_{-6.1}^{+7.1}$
depth (ppm)	4327 ± 929	${2975}_{-531}^{+672}$	${3096}_{-437}^{+420}$
P (days)	65.2 ± 17.8	${21.8}_{-5.9}^{+7.7}$	${22.5}_{-4.6}^{+3.7}$
T_dur (hr)	6.1 ± 0.8	${3.0}_{-0.2}^{+0.3}$	${3.2}_{-0.3}^{+0.2}$
$\| {\dot{T}}_{\mathrm{dur}}\|$ (minutes yr⁻¹)	6.4 ± 1.4	${3.6}_{-0.9}^{+0.8}$	${4.2}_{-0.8}^{+0.9}$

System properties:
N_obs	2.0 ± 0.3	${1.8}_{-0.2}^{+0.2}$	${1.6}_{-0.2}^{+0.3}$
N_phys	⋯	${3.8}_{-0.4}^{+0.5}$	${5.6}_{-0.4}^{+0.4}$
σ_i (deg)	⋯	${1.3}_{-0.2}^{+0.3}$	${16.7}_{-6.2}^{+6.5}$

Catalog properties:
Single-to-multi ratio	0.78	${0.91}_{-0.3}^{+0.5}$	${1.77}_{-0.6}^{+0.8}$

Note. Means of the distributions of several planetary and system properties of planets with detected TDVs. The left column corresponds to the observed KOIs listed in Table 1. We report the mean and standard error of the mean. The middle and right columns correspond to the SysSim simulated planets in the maximum AMD model and two-Rayleigh model, respectively. We calculate the mean values associated with each of the 100 catalogs and report the medians and 16th and 84th percentiles of the distributions of means. Here N_obs is the observed transit multiplicity of the system, and N_phys is the intrinsic multiplicity of the system. The quantity σ_i is the standard deviation of the system's orbital inclinations with respect to the invariable plane, including nondetected planets. N_phys and σ_i are unknown for the observed systems. The final row is the ratio of the number of TDV detections that are in single-transiting systems compared to multiple-transiting systems.

Download table as: ASCII Typeset image

Figure 8 presents histograms of the absolute value of the TDV slope, $| {\dot{T}}_{\mathrm{dur}}|$ , for planets with detected TDVs. The range and typical values of $| {\dot{T}}_{\mathrm{dur}}|$ for both the two-Rayleigh model and the maximum AMD model agree well with the observed planets. In all three distributions (observations and two models), the majority of values of $| {\dot{T}}_{\mathrm{dur}}|$ fall in the range of 1 − 10 minutes yr⁻¹. The simulations have a slightly greater proportion of detections with $| {\dot{T}}_{\mathrm{dur}}| \lesssim 1\,\mathrm{minutes}\,{\mathrm{yr}}^{-1}$ . The two-Rayleigh model is inconsistent in terms of the observed count, but this is just a restatement of the finding from Figure 5 that the two-Rayleigh model produces too many detected TDVs.

Figure 9 shows additional dimensions of this comparison between simulated and observed planets. We plot $| {\dot{T}}_{\mathrm{dur}}|$ versus P and ${\mathrm{log}}_{10}(\mathrm{depth})$ for planets with detected TDVs in the maximum AMD model, the two-Rayleigh model, and the Kepler observations. The left column illustrates that $| {\dot{T}}_{\mathrm{dur}}|$ is positively correlated with P for both the simulated and observed planets.

It is illuminating to discuss the physical origins of this positive correlation. In the limit of circular orbits, the expression for ${\dot{T}}_{\mathrm{dur}}$ (Equation (4)) becomes

$\begin{eqnarray}&&{\dot{T}}_{\mathrm{dur}}\approx -P\left(\displaystyle \frac{b}{\pi \sqrt{1-{b}^{2}}}\right)\dot{{\rm{\Omega }}}\sin i\cos \beta \sin {\rm{\Omega }}.\end{eqnarray} \tag{ 19 }$

If the typical period ratio between adjacent planets, P_i+1/P_i, is independent of P, then $| \dot{{\rm{\Omega }}}| \propto {P}^{-1}$ , and the dependence of $| {\dot{T}}_{\mathrm{dur}}|$ on P vanishes.¹³ Indeed, $| {\dot{T}}_{\mathrm{dur}}|$ is completely uncorrelated with P in the full distribution (that includes the TDV slopes that are too small to be detected). However, ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ is positively correlated with P (due to the decreasing number of transits with increasing P), and since TDV detection requires $| {\dot{T}}_{\mathrm{dur}}| \gt 3{\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ (Section 3.1), this conspires to yield a positive correlation between $| {\dot{T}}_{\mathrm{dur}}|$ and P. That is, for large P, only the steepest values of $| {\dot{T}}_{\mathrm{dur}}|$ are detectable.

Shahaf et al. (2021) also identified a positive slope between $| {\dot{T}}_{\mathrm{dur}}|$ and P within the Kepler-detected TDVs, and they additionally pointed to a paucity of large $| {\dot{T}}_{\mathrm{dur}}|$ detections at small P, the latter of which might not be explained by the ${\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ bias. The simulated planets at small P in Figure 9 show a slight deficiency of large $| {\dot{T}}_{\mathrm{dur}}|$ detections (∼10 minutes yr⁻¹) compared to more moderate slopes (∼1–3 minutes yr⁻¹), which is due to the fact that moderate values are more common in the underlying distribution. However, it is not clear that the deficiency in the simulations matches that of the data. Altogether, the simulations reproduce the correlation, which is the dominant feature of the data, but they do not clearly reproduce the paucity of large $| {\dot{T}}_{\mathrm{dur}}|$ at small P.

The $| {\dot{T}}_{\mathrm{dur}}|$ versus P distribution also shows that the simulated planets are more concentrated at smaller orbital periods in the range P = 3–10 days than the observed planets, which are mostly at P > 10 days. In particular, the observations have more planets with detected TDVs with P = 100–300 days than represented in the simulations. This may be small number statistics of the data, or it could be because SysSim is less well calibrated in the range of P = 100–300 days owing to fewer Kepler planet detections there. In addition, the observed systems may contain perturbing companion planets with P > 300 days, but the simulated systems do not. We will return to this in Section 4.3.

The right column of Figure 9 shows that $| {\dot{T}}_{\mathrm{dur}}|$ versus ${\mathrm{log}}_{10}(\mathrm{depth})$ exhibits a weak negative correlation, which is due to the fact that a larger transit S/N can allow for the detection of a smaller TDV slope. The typical transit depths of the simulated planets are in general agreement with the observed planets, although the observed planets appear to have a larger spread in ${\mathrm{log}}_{10}(\mathrm{depth})$ than the simulated planets. Altogether, these results indicate that the agreement between the simulated and observed planets with detected TDVs is satisfactory.

Table 2 shows a quantitative comparison of the average properties of the planets with detected TDVs (and their systems) between the SysSim simulated planets and the observed KOIs. The table summarizes several of the key takeaways from Figures 8 and 9, including the general agreement of $| {\dot{T}}_{\mathrm{dur}}|$ between simulations and observations and the larger average P for the observed planets. In addition, we note that the observed transit multiplicity of systems with detected TDVs shows good agreement, with ${N}_{\mathrm{obs}}={1.8}_{-0.2}^{+0.2}$ for the simulated systems in the maximum AMD model versus 2.0 ± 0.3 for the observed systems. The average intrinsic multiplicity of the maximum AMD model simulated systems hosting planets with detected TDVs is ${N}_{\mathrm{phys}}={3.8}_{-0.4}^{+0.5}$ (but clearly unknown for the observed systems). This is marginally smaller than the average of the full population of simulated systems with observed planets, which is N_phys ∼ 4.5.¹⁴ The maximum AMD model simulated systems hosting planets with detected TDVs have an average standard deviation of the inclinations (measured with respect to the invariable plane) equal to ${\sigma }_{i}={1.3}_{-0.2}^{+0.3}$ deg, marginally larger than the average of the full population, σ_i ∼ 1 deg.

4.3. Sensitivity of Results to Model Assumptions

In this section, we investigate the robustness of our results with respect to variations in the models, specifically the planet radius and period ranges, and H20's assumption that systems are at the critical AMD. To begin, we note that SysSim considers planet radii in the range 0.5 R_⊕ < R_p < 10 R_⊕, but it is best constrained for R_p < 4 R_⊕ owing to the large number of detections of sub-Neptune-sized planets. We can repeat our analysis by changing the R_p upper limit from 10 R_⊕ to 4 R_⊕. This reduces the number of observed Kepler TDV detections from 16 (Section 3.1) to 6. Figure 10 (top panel) displays the distributions of the number of simulated planets with detected TDVs, as in Figure 5. We observe a similar or perhaps even better agreement between the maximum AMD model and the observed number of planets with detected TDVs. Meanwhile, the two-Rayleigh model again produces too many detected TDVs. Our results are therefore robust with respect to the choice of the upper limit on R_p.

**Figure 10.** Same as Figure 5, except with the planet radius upper limit equal to 4 R_⊕ rather than 10 R_⊕ (both panels) and with the orbital period upper limit equal to 100 days rather than 300 days (bottom panel).
Download figure:
Standard image High-resolution image

In addition to the finite planet radius range, SysSim also uses a finite range in orbital periods, 3 days < P < 300 days. This leads to an important limitation of our analysis, in that our population of simulated TDV detections is missing cases that would have arisen from planetary perturbers with P < 3 days or P > 300 days. This may be responsible for the underabundance of simulated TDV detections with P ∼ 100–300 days relative to the observations (Figure 9). Accounting for this underprediction of planets with detected TDVs in the simulated population would tend to further disfavor the two-Rayleigh model, although it may also lead to a tension with the maximum AMD model if the difference is significant. The best way to address this would be to generalize the H20 model to include planets with orbital periods beyond the current cutoffs and then use SysSim to generate new simulated catalogs for comparison. Generalizing the H20 model and performing parameter estimation on the new model's parameters is beyond the scope of this study. However, we can repeat the analysis on a restricted sample of TDV detections with P < 100 days (bottom panel of Figure 10), which is less affected by the 300-day cutoff. We find that the overall conclusions are unchanged.

Another simulation feature to consider is the assumption of the maximum AMD model that all multiplanet systems are at the AMD stability limit, with AMD = AMD_crit. H20 tested variations of their model in which AMD = f_crit × AMD_crit, where f_crit is some factor in the range [0, 2]. They found that values between 0.4 and 2 are all acceptable and that the fit did not significantly improve with the extra f_crit parameter. However, it is possible that this factor could affect the TDV distribution more strongly than the other Kepler observables. We find this to be unlikely. As shown by Figure 7 in H20, when allowing for f_crit ≠ 1, the mutual inclination distribution maintains the same shape and shifts by less than a factor of two. Specifically, the median mutual inclination shifts to ∼0.82° (1.64°) for f_crit = 0.5 (2). Meanwhile, $| {\dot{T}}_{\mathrm{dur}}|$ varies by up to three orders of magnitude among the planets with detected TDVs (with $| {\dot{T}}_{\mathrm{dur}}|$ ranging from ∼0.1 to 100 minutes yr⁻¹; see Figure 9). Including the nondetected TDVs, there is an even larger range, with $| {\dot{T}}_{\mathrm{dur}}|$ as low as ∼10⁻⁴ minutes yr⁻¹. Given the dynamic range in $| {\dot{T}}_{\mathrm{dur}}|$ , the small changes in inclinations cannot significantly change the distribution of TDVs.

5. Discussion

5.1. Small Mutual Inclinations for the Kepler Planets

The TDV statistics of the Kepler population agree very well with the simulated planet population constructed by the maximum AMD model, whereas the two-Rayleigh model produces too many detected TDVs (Figures 5 and 10). Given that the primary difference between these two models is the mutual inclination distribution, which TDVs are sensitive to (e.g., Figures 6 and 7), we conclude that the TDV statistics support the mutual inclination distribution of the maximum AMD model. We emphasize that the TDV statistics are not capable of assessing H20's assumption that systems are at the critical AMD. However, they appear consistent with the implications of that assumption.

As discussed in Section 3.2.2 (and originally in Section 3.3 of H20), the maximum AMD model yields median mutual inclinations, ${\tilde{\mu }}_{i,n}$ , of systems with n = 2, 3,..., 10 planets that are well modeled by a power law of the form ${\tilde{\mu }}_{i,n}=1.10^\circ {(n/5)}^{-1.73}$ . This equation evaluates to 5 fdg 4, 2 fdg 7, and 1 fdg 6 for two-, three- and four-planet systems, respectively. These inclinations are quite modest relative to the high mutual inclinations that have previously been invoked to explain the Kepler dichotomy (e.g., Mulders et al. 2018, H19). Of course, systems with higher mutual inclinations do exist (e.g., Kepler-108; Mills & Fabrycky 2017), but according to our results, they do not make up a large fraction of the population.

A final point to emphasize about H20's mutual inclination distribution is that it is nondichotomous; it does not bifurcate the population into "dynamically cool" and "dynamically hot" systems. The fact that the TDV statistics support this model is therefore evidence for a continuum of architectures rather than a dichotomy. As a caveat, we note that there are many possible combinations in the true multiplicity/inclination distribution that are not encapsulated in the two models we tested. Reality is likely more complicated than "there is" or "there is not" a dichotomy. However, we can say that (1) the continuous maximum AMD model is consistent with the data and (2) replacing any of the model's low-inclination systems with high-inclination configurations would tend to increase the number of TDV detections, thus leading to eventual inconsistency with the small number of observed detections. The two models tested here serve as a basis of comparison for future, more complex models. For now, we focus on the implications of the favored maximum AMD model.

The mutual inclinations of the maximum AMD model are consistent with previous works that constrained the distribution by supplementing the Kepler transit statistics with additional observations. The most similar result is that of Zhu et al. (2018; see also Yang et al. 2020), who used Kepler TTV statistics as their additional constraint. The Zhu et al. (2018) model assumed a power law of the mutual inclination dispersion with the intrinsic multiplicity, σ_i,n = σ_i,5(n/5)^α, and identified best-fit parameters equal to α = −3.5 and σ_i,5 = 0 fdg 8. This is steeper than that found by H20 but a similar result indicating relatively small inclinations that depend on the intrinsic multiplicity.

Rather than TTVs, Tremaine & Dong (2012) and Figueira et al. (2012) used RV survey data as their extra observational constraint. By combining statistical information from Kepler and RV surveys, they found evidence that suggested that the mean mutual inclinations are ≲5°. Detailed analyses of RVs of individual multiplanet systems with strong gravitational interactions also provide evidence for small mutual inclinations (Laughlin et al. 2002; Nelson et al. 2014, 2016; Millholland et al. 2018).

Finally, the inference of small mutual inclinations of the Kepler planets is also consistent with population-level constraints on Kepler stellar obliquities. Studies using photometric variability observations (Mazeh et al. 2015) and $v\sin i$ measurements (Winn et al. 2017; Louden et al. 2021) showed that stellar spin–orbit misalignments are small for Kepler planets orbiting cool stars (T_eff ≲ 6250 K). All of these results are broadly consistent with H20's maximum AMD model. Thus, we are seeing agreement between Kepler transit statistics, TTVs, RVs, stellar obliquities, and TDVs, indicating a collection of robust evidence that the vast majority of inner planetary systems around Sun-like stars have small mutual inclinations.

5.2. Physical Origins of the Mutual Inclinations

The predominantly small (but nonzero) inclinations, i ≲ 5°–10°, must have been excited through one or more physical processes. While several dynamical mechanisms have been proposed, it is still unclear which dominate. Here we will review the relevant theories and discuss which are favored by the observational evidence for small inclinations.

5.2.1. Excitation during Late-stage Planet Assembly

To begin, we note that the requisite mutual inclinations can probably be acquired primordially during the planet formation epoch. Many studies have used N-body simulations to study the formation of close-in, compact systems in their final assembly stage, during which planetary embryos undergo mutual scatterings and collisional growth. Some of these models have considered in situ formation in a gas-poor or gas-empty environment, typically starting with a dense planetesimal disk with various disk masses and radial profiles (e.g., Hansen & Murray 2013; Dawson et al. 2016; Moriarty & Ballard 2016; Matsumoto & Kokubo 2017; Poon et al. 2020). In contrast to in situ formation models, migration models consider formation of protoplanets beyond ∼1 au. Disk–planet interactions lead to inward migration and the formation of long chains of short-period protoplanets in MMRs. Once the gas disk disperses, these compact resonant chains can become dynamically unstable and collide, promoting further planetary growth (e.g., Terquem & Papaloizou 2007; Coleman & Nelson 2014; Cossou et al. 2014; Izidoro et al. 2017, 2019; Carrera et al. 2018, 2019).

In both the in situ and migration scenarios, gravitational scattering of planetary embryos leads to dynamical excitation of the mutual inclinations and eccentricities of the final planetary system, potentially allowing for the production of systems that agree with the Kepler multiplicity distribution.¹⁵ These models are generally capable of producing inclinations in the range of ∼1°–10° (e.g., Hansen & Murray 2013; Moriarty & Ballard 2016; Izidoro et al. 2019; Carrera et al. 2019; Poon et al. 2020), in agreement with the mutual inclinations of the maximum AMD model from H20. For instance, the migration simulations of Izidoro et al. (2019) find an excess of single-transiting systems that matches the Kepler multiplicity distribution and arises primarily from systems of two to three planets with inclinations between ∼4° and 10°. In contrast to Izidoro et al. (2019), some other variations of late-stage planet formation simulations do not reproduce the large number of single-transiting systems, at least not with a uniform underlying model (e.g., Hansen & Murray 2013; Poon et al. 2020). However, this appears equally if not more influenced by the intrinsic multiplicities of the simulated systems being too high (average of ∼5 planets per system rather than ∼3 as found by H20) than by the inclinations being too low, and the multiplicities are strongly influenced by the initial conditions of the simulations. Finally, in addition to the proper range of inclinations, these models also naturally produce an inverse correlation between inclination (and eccentricity) dispersion and intrinsic multiplicity, as shown in the analysis of the simulations of Carrera et al. (2018) by H20.

5.2.2. External Perturbations: Stellar Oblateness and Distant Giant Planets

While dynamical instabilities during the late-stage planet assembly process appear capable of delivering the required level of mutual inclination excitation, this is not guaranteed, and the primordial formation conditions are poorly constrained. This has led some to postulate the importance of external gravitational perturbations from the star or other planets in the system as sources of dynamical excitation. These external perturbations additionally offer solutions to systems with very high (≳20°) mutual inclinations (e.g., Mills & Fabrycky 2017) and/or misaligned stellar obliquities (e.g., Kunovac Hodžić et al. 2021), which are observed in a handful of systems and cannot be attributed to primordial collisional excitation among the inner Kepler planets.

One source of external perturbation is the quadrupolar gravitational potential of a tilted host star, which is particularly strong when the star is young and rapidly rotating (Spalding & Batygin 2016). Shortly after disk dispersal, a tilted star can drive differential nodal precession of the inner planetary orbits that leads to mutual inclinations comparable to the magnitude of the stellar obliquity. The excitation also depends on the coupling of the Kepler planets compared to the stellar forcing. To match the ≲10° mutual inclinations, spin–orbit misalignments of this magnitude must be widespread during the disk-hosting stage. There are numerous pathways for generating star–disk tilts (e.g., Batygin 2012; Lai 2014; Spalding & Batygin 2014; Spalding et al. 2014; Zanazzi & Lai 2018). Moreover, the 6° obliquity of the Sun is of the requisite magnitude, and recent data suggest that spin–orbit misalignments in systems of small Kepler planets are more common than previously thought, particularly among hot stars (Louden et al. 2021). The oblateness-driven inclination excitation theory is also consistent with the relatively large (up to ∼10°) observed inclinations of ultra-short-period planets (Dai et al. 2018; Li et al. 2020), provided that their inward migration occurs early enough (Millholland & Spalding 2020). The biggest unknown with the theory is that it requires rapid disk dispersal timescales, such that a primordial star–disk misalignment can naturally become a spin–orbit misalignment (Spalding & Millholland 2020).

In addition to the host star, another gravitational source that can drive differential precession is a distant (≳1 au) giant planet (or "cold Jupiter") on an inclined orbit (Becker & Adams 2017; Lai & Pu 2017; Hansen 2017; Pu & Lai 2018). The giant can similarly generate mutual inclinations among the inner planets up to roughly the magnitude of the giant's orbital tilt, depending on the coupling of the inner planets to one another relative to the perturbations from the giant. The inclination excitation can, however, become much larger than the giant's inclination when there exists a secular resonance between a nodal precession frequency and a system eigenfrequency (Granados Contreras & Boley 2018; Pu & Lai 2018). In addition to dynamically quiet secular interactions, an outer system containing multiple giant planets can undergo a violent epoch of planet–planet scattering and collisions that reduce the multiplicity and/or dynamically heat the inner system (e.g., Gratia & Fabrycky 2017; Mustill et al. 2017; Poon et al. 2020). Recently, Masuda et al. (2020) showed that cold Jupiters have ∼10° mutual inclinations relative to inner transiting systems, with lower inclinations relative to the multis than singles. This is consistent with a scenario in which inclined giants perturb inner systems. However, it is equally consistent with other mechanisms that generate inclinations among the inner planets and produce an inclination with respect to the cold Jupiter as a by-product.

Cumulatively, the data suggest that distant giants likely play a role in dynamically heating some subset of inner planetary systems, with prime examples including HAT-P-11 (Yee et al. 2018), π Mensae (e.g., Xuan & Wyatt 2020; Kunovac Hodžić et al. 2021), and WASP-107 (e.g., Piaulet et al. 2021), all of which have short-period Neptunes (or sub-Neptunes) with large stellar obliquities and eccentric cold Jupiters. However, as far as generating the Kepler dichotomy through dynamical excitation of inner systems, distant giant planets are unlikely to be the dominant solution. There are several reasons for this. First, roughly three-quarters of Kepler planetary systems of super-Earths/sub-Neptunes have just a single transiting planet, while roughly one-third of these systems contain distant giant planets (Zhu & Wu 2018; Bryan et al. 2019). Provided that the distant giant fraction is similar around transit singles and transit multis (Masuda et al. 2020), the statistical constraints make it difficult to generate a sufficient number of transit singles through giant planet perturbations. Moreover, it is unlikely that the cold Jupiter occurrence is significantly higher among the transit singles compared to the transit multis, given that the stellar metallicity distributions are indistinguishable (Munoz Romero & Kempton 2018).

A second consideration is that, even when one or more cold Jupiters are present, they do not always (in the absence of secular resonances) lead to enhanced mutual inclinations in inner systems, since the inner planet orbits are often tightly coupled (e.g., Gratia & Fabrycky 2017; Mustill et al. 2017). However, when mutual inclinations do result, they are generally larger than the i ≲ 5°–10° scale required by H20, particularly when there are scattering interactions in the outer system (Mustill et al. 2017). Finally, the gravitational perturbations from distant giants are often overpowered by stellar oblateness soon after disk dispersal (Spalding & Millholland 2020), such that a planetary system initially dominated by the star can evade giant-induced inclination excitation by adiabatically realigning to the giant's orbital plane during stellar spin-down.

6. Conclusion

The distribution of inclinations in multiplanet systems encodes fundamental clues about planet formation. Unfortunately, it is inherently difficult to infer from RV and transit observations. In particular, the observed transiting multiplicity distribution is strongly affected by the underlying mutual inclination distribution, but it also depends on the intrinsic multiplicity distribution, yielding a near degeneracy that makes inferences of inclinations difficult. In this work, we used TDVs of the Kepler planet population to break the near degeneracy and constrain the mutual inclination distribution of close-in, multiplanet systems. TDVs often arise when a transiting planet's transit chord drifts owing to orbital precession induced by torques from perturbing planets (Figure 1). The signal is sensitive to the mutual orbital inclinations and can yield detectable drifts on the order of ${\dot{T}}_{\mathrm{dur}}\sim 1\mbox{--}10\,\mathrm{minutes}\,{\mathrm{yr}}^{-1}$ . Several dozen Kepler planets exhibit TDV drift signals (Table 1; Figure 2 for an example). Our work is the first to exploit these detections statistically to characterize the mutual inclination distribution.

We compared the observed TDV detections of Kepler planets to expectations from simulated planet populations subject to different assumptions about the mutual inclination distribution. These simulated planetary systems were drawn from two population models built using the "SysSim" empirically calibrated forward modeling framework: (1) the "two-Rayleigh model" (He et al. 2019), which assumes a dichotomous mutual inclination distribution with low-inclination (σ_i,low ∼ 1°–2°) and high-inclination (σ_i,high ∼ 30°–65°) components, and (2) the "maximum AMD model" (He et al. 2020), in which the mutual inclination distribution is broad, continuous, and multiplicity dependent, with small inclinations on the scale of a few degrees. (See Figure 3 for a comparison of the two models.) To analyze the simulated planet population, we considered both analytic and N-body calculations of ${\dot{T}}_{\mathrm{dur}}$ , along with a simulated TDV detection pipeline to identify which planets would have TDVs that are both measurable (requiring transits with sufficiently large S/N) and detectable (requiring a significant TDV slope, $| {\dot{T}}_{\mathrm{dur}}| \gt 3\ {\sigma }_{{\dot{T}}_{\mathrm{dur}}}$ ).

Our main result is shown in Figure 5, where the key diagnostic is the number of planets with detected TDVs. The maximum AMD model yields a quantity ( ${22}_{-6}^{+10}$ ) that is in good agreement with the observed number from Kepler (16 after cuts have been applied). The two-Rayleigh model, by contrast, consistently overpredicts the number of planets with detected TDVs ( ${43}_{-13}^{+18}$ ). This is because it has too many high-inclination systems, which generally produce larger TDV signals (Figure 7). These results are robust with respect to model assumptions, such as the upper limit of R_p. When restricting to planets with R_p < 4 R_⊕ (rather than R_p < 10 R_⊕), there are six observed Kepler planets with detected TDVs, compared to ${8}_{-3}^{+4}$ with the maximum AMD model and ${16}_{-5}^{+7}$ with the two-Rayleigh model (Figure 10).

Given these results, our key takeaway is that the TDV statistics support a continuous distribution of relatively low mutual inclinations (i ≲ 5°–10°) rather than a dichotomous distribution with many high-inclination systems. These results are consistent with Zhu et al. (2018), who considered TTV statistics rather than TDVs. Moreover, small inclinations are also supported by studies that combined Kepler and RV survey statistics (Figueira et al. 2012; Tremaine & Dong 2012). Cumulatively, this work is further evidence that the apparent excess of single-transiting systems relative to expectations from the multis, an observation known as the "Kepler dichotomy," does not actually provide evidence for a dichotomy in the underlying architectures (He et al. 2020). Rather, the observations are naturally explained by and more consistent with a "Kepler continuum" of intrinsic multiplicities and low mutual inclinations.

Even with predominantly small (few-degree) mutual inclinations, there must have been one or more physical processes that disrupted these close-in systems from perfect coplanarity. Primordial excitation during the planet assembly phase is perhaps the most favored mechanism for producing a nondichotomous, multiplicity-dependent distribution of low mutual inclinations, although the optimal disk surface densities, gas damping timescales, and planetary embryo properties are still uncertain. Stellar oblateness can also drive small mutual inclinations if ∼5°–10° spin–orbit misalignments (like the 6° solar obliquity) are widespread and disk dispersal timescales are rapid. Both of these are observable properties that will become better understood over time. Distant giant planets are almost certainly responsible for dynamical heating of a subset of inner planetary systems, but they are probably not the dominant origin of the transit-reducing mutual inclinations. Regardless, the results in this paper can be exploited in future work to constrain the prevalence of these various mechanisms.

Statistical characterization of TDV signals will improve in the future as observational time baselines lengthen. In particular, the PLATO mission (Rauer et al. 2014) will enable the detection of many more Kepler planet TDVs when it revisits the Kepler field. Moreover, some planets may be seen to transit into or out of view (Freudenthal et al. 2018), as has already been observed in a handful of cases (Hamann et al. 2019; Judkovsky et al. 2020). Modeling this expanded set of TDV observations will allow an even more detailed characterization of the mutual inclination distribution.

In summary, TDVs provide a window into the three-dimensional properties of planetary systems, which are otherwise difficult to probe. This is another axis with which to study our solar system in the context of exoplanetary systems, and, in this case, it appears that the near coplanarity of the solar system is indeed the rule.

We thank Tsevi Mazeh and Dimitri Veras for comments and questions that improved the paper. We also thank the anonymous referee for their insightful and helpful report. S.C.M. was supported by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51465 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. M.Y.H. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference No. PGSD3-516712-2018. This work was supported by a grant from the Simons Foundation/SFARI (675601, E.B.F.). E.B.F. acknowledges the support of the Ambrose Monell Foundation and the Institute for Advanced Study. M.Y.H. and E.B.F. acknowledge support from the Penn State Eberly College of Science and Department of Astronomy & Astrophysics, the Center for Exoplanets and Habitable Worlds, and the Center for Astrostatistics.

Appendix A: Fourth-order Expansion of the Secular Disturbing Function

Here we provide the expression for the secular disturbing function expansion of two point-mass planets with masses m and $m^{\prime}$ . The semimajor axes are a and $a^{\prime}$ , with $a\lt a^{\prime}$ . The disturbing functions for the inner planet and outer planet, respectively, are given by

$\begin{eqnarray}\begin{array}{rcl}{ \mathcal R } & = & \displaystyle \frac{\mu ^{\prime} }{a^{\prime} }{{ \mathcal R }}_{{\rm{D}}}+\displaystyle \frac{\mu ^{\prime} }{a^{\prime} }\alpha {{ \mathcal R }}_{{\rm{E}}}\\ { \mathcal R }^{\prime} & = & \displaystyle \frac{\mu }{a^{\prime} }{{ \mathcal R }}_{{\rm{D}}}+\displaystyle \frac{\mu }{a^{\prime} }\displaystyle \frac{1}{{\alpha }^{2}}{{ \mathcal R }}_{{\rm{I}}},\end{array}\end{eqnarray} \tag{ A1 }$

where $\mu \equiv { \mathcal G }m$ and $\alpha \equiv a/a^{\prime}$ . ${{ \mathcal R }}_{{\rm{D}}}$ is the direct part of the disturbing function, and ${{ \mathcal R }}_{{\rm{E}}}$ and ${{ \mathcal R }}_{{\rm{I}}}$ are the indirect parts. In our case, the indirect parts are zero because there are no secular terms. As for ${{ \mathcal R }}_{{\rm{D}}}$ , we retain secular terms up to fourth order in e and $s\equiv \sin (i/2)$ . These terms can be found in the Appendix tables of Murray & Dermott (1999), yielding

$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal R }}_{D}={f}_{1}+{f}_{2}({e}^{2}+e{{\prime} }^{2})+{f}_{3}({s}^{2}+s{{\prime} }^{2})+{f}_{4}{e}^{4}+{f}_{5}{e}^{2}e{{\prime} }^{2}\\ \quad +\,{f}_{6}e{{\prime} }^{4}+{f}_{7}({e}^{2}{s}^{2}+e{{\prime} }^{2}{s}^{2}+{e}^{2}s{{\prime} }^{2}+e{{\prime} }^{2}s{{\prime} }^{2})\\ \quad +\,{f}_{8}({s}^{4}+s{{\prime} }^{4})+{f}_{9}{s}^{2}s{{\prime} }^{2}+[{f}_{10}{ee}^{\prime} +{f}_{11}{e}^{3}e^{\prime} \\ \quad +\,{f}_{12}{ee}{{\prime} }^{3}+{f}_{13}{ee}^{\prime} ({s}^{2}+s{{\prime} }^{2})]\cos (\varpi ^{\prime} -\varpi )\\ \quad +\,[{f}_{14}{ss}^{\prime} +{f}_{15}{ss}^{\prime} ({e}^{2}+e{{\prime} }^{2})+{f}_{16}{ss}^{\prime} ({s}^{2}+s{{\prime} }^{2})]\\ \quad \cos ({\rm{\Omega }}^{\prime} -{\rm{\Omega }})+{f}_{17}{e}^{2}e{{\prime} }^{2}\cos (2\varpi ^{\prime} -2\varpi )\\ \quad +\,{f}_{18}{e}^{2}{s}^{2}\cos (2\varpi -2{\rm{\Omega }})+{f}_{19}{ee}^{\prime} {s}^{2}\\ \quad \cos (\varpi ^{\prime} +\varpi -2{\rm{\Omega }})+{f}_{20}e{{\prime} }^{2}{s}^{2}\cos (2\varpi ^{\prime} -2{\rm{\Omega }})\\ \quad +\,{f}_{21}{e}^{2}{ss}^{\prime} \cos (2\varpi -{\rm{\Omega }}^{\prime} -{\rm{\Omega }})\\ \quad +\,{f}_{22}{ee}^{\prime} {ss}^{\prime} \cos (\varpi ^{\prime} -\varpi -{\rm{\Omega }}^{\prime} +{\rm{\Omega }})\\ \quad +\,{f}_{23}{ee}^{\prime} {ss}^{\prime} \cos (\varpi ^{\prime} -\varpi +{\rm{\Omega }}^{\prime} -{\rm{\Omega }})\\ \quad +\,{f}_{24}{ee}^{\prime} {ss}^{\prime} \cos (\varpi ^{\prime} +\varpi -{\rm{\Omega }}^{\prime} -{\rm{\Omega }}).\end{array}\end{eqnarray} \tag{ A2 }$

The f_i coefficients in this expression are functions of α that depend on the Laplace coefficients and their derivatives.

The disturbing function expansion in Equation (A2) is written explicitly for two planets. It is straightforward to generalize this treatment to N-planet systems. For each planet in the system, the disturbing functions associated with each pairwise planet–planet interaction are summed, using ${ \mathcal R }$ when the perturbing planet is external and ${ \mathcal R }^{\prime}$ when the perturbing planet is internal.

Appendix B: Comparison of Analytic and N-Body TDV Calculation

To assess the performance of the analytic calculation of ${\dot{T}}_{\mathrm{dur}}$ , here we show how it compares to a computation using a direct N-body integration. The analytic approach was described in Section 2. The N-body calculation uses the REBOUND gravitational dynamics software (Rein & Liu 2012) and was described in Section 3.3.2. Rather than build a set of systems on which to test the calculations, we simply use the SysSim systems from both the maximum AMD model and the two-Rayleigh model.

Figure 11 shows the analytic ${\dot{T}}_{\mathrm{dur}}$ versus the N-body ${\dot{T}}_{\mathrm{dur}}$ for a set of systems generated by the maximum AMD model, and Figure 12 shows a histogram of the fractional difference. Overall, we observe a fairly strong agreement. The best half of the distribution has analytic ${\dot{T}}_{\mathrm{dur}}$ values within 25% of the N-body values; the best three-quarters of the distribution agree within 50%. Figure 13 shows the fractional difference between the analytic and N-body values as a function of the period ratio between the planet for which ${\dot{T}}_{\mathrm{dur}}$ is calculated and its nearest neighbor. The analytic calculation is less accurate for smaller period ratios and for planets near MMRs. This is expected because the disturbing function expansion leaves out resonant terms. Figure 14 shows the comparisons for the two-Rayleigh model. There are orders-of-magnitude discrepancies between the analytic and N-body calculations at large inclinations, where the accuracy of the secular expansion breaks down.