Stellar Oblateness versus Distant Giants in Exciting Kepler Planet Mutual Inclinations

In order to describe our mathematical approach, we begin with n_p = 2 inner planets, but we generalize the method in the next section. We assume that the planets' orbital periods are far from integer ratios, allowing the use of secular dynamics, whereby each orbit is effectively approximated as a massive wire (Morbidelli 2002). Furthermore, we restrict attention to small inclinations, such that Laplace–Lagrange secular theory is appropriate (Murray & Dermott 1999).

The interior planets are influenced by secular perturbations from the exterior giant, the host star's quadrupole, and the other close-in planets. The ratio of the stellar angular momentum to the orbital angular momentum of planet p is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{l}_{\star }{M}_{\star }{\omega }_{\star }{R}_{\star }^{2}}{{m}_{p}\sqrt{{{GM}}_{\star }{a}_{p}}} & \sim & 10{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{{P}_{\star }}{10\mathrm{day}}\right)}^{-1}\\ & & \times \ {\left(\displaystyle \frac{{a}_{p}}{0.3\mathrm{au}}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{{m}_{p}}{10{M}_{\oplus }}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 1 }$$

where we define the stellar mass, M_⋆, radius, R_⋆, spin rate, ω_⋆ = 2π/P_⋆, and moment of inertia factor, l_⋆ = 0.2 (Batygin & Adams 2013). This ratio is much greater than unity within the regime we consider. Similarly, the distant giant's orbital angular momentum greatly exceeds that of the inner planetary system. Thus, we assume that both the giant's orbit and the stellar orientation are fixed throughout the secular calculation (an assumption that is relaxed in Section 6).

The evolution of the system is conveniently described by Hamiltonian dynamics in terms of the canonical Poincaré variables (Morbidelli 2002):

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{i} & \equiv & {m}_{i}\sqrt{{{GM}}_{\star }{a}_{i}}\left[1-\cos ({I}_{i})\right]\approx \displaystyle \frac{{I}_{i}^{2}}{2}{m}_{i}\sqrt{{{GM}}_{\star }{a}_{i}}\\ {{\rm{\Lambda }}}_{i} & \equiv & {m}_{i}\sqrt{{{GM}}_{\star }{a}_{i}}\\ {z}_{i} & \equiv & -{{\rm{\Omega }}}_{i},\end{array}\end{eqnarray} \tag{ 2 }$$

where Ω_i are the longitudes of ascending node.

Without loss of generality, we set Ω_⋆ = 0 and Z_G = 0. In the case thus described, we may write the Laplace–Lagrange Hamiltonian to lowest order as (Murray & Dermott 1999)

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal H }=\displaystyle \frac{3}{2}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{1}}\right)}^{2}\displaystyle \frac{{{GMm}}_{1}}{{a}_{1}}\left[\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}-{\beta }_{\star }\sqrt{\displaystyle \frac{2{Z}_{1}}{{{\rm{\Lambda }}}_{1}}}\cos ({z}_{1})\right]\\ +\displaystyle \frac{3}{2}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{2}}\right)}^{2}\displaystyle \frac{{{GMm}}_{2}}{{a}_{2}}\left[\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}}-{\beta }_{\star }\sqrt{\displaystyle \frac{2{Z}_{2}}{{{\rm{\Lambda }}}_{2}}}\cos ({z}_{2})\right]\\ +\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{4{a}_{2}}{b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\left[\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}+\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}}\right.\\ -\left.2\sqrt{\displaystyle \frac{{Z}_{1}{Z}_{2}}{{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}}}\cos ({z}_{1}-{z}_{2})\right]\\ +\displaystyle \frac{3{{Gm}}_{G}{m}_{2}}{4{a}_{G}}{\left(\displaystyle \frac{{a}_{2}}{{a}_{G}}\right)}^{2}\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}+\displaystyle \frac{3{{Gm}}_{G}{m}_{2}}{4{a}_{G}}{\left(\displaystyle \frac{{a}_{2}}{{a}_{G}}\right)}^{2}\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}},\end{array}\end{eqnarray} \tag{ 3 }$$

where the first two terms represent the gravitational interaction between the planets and the host star's quadrupole. The third term represents the interaction between the inner two planets, and the final two terms arise from the distant giant's coupling to the two close-in planets. Note that we have made the assumption that a_G ≫ a_i.

It is convenient to introduce the following frequencies:

$$\begin{eqnarray}\begin{array}{rcl}{B}_{12} & \equiv & {n}_{1}\displaystyle \frac{{m}_{2}}{4{M}_{\star }}{b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right){\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)}^{2}\\ {B}_{21} & \equiv & {n}_{2}\displaystyle \frac{{m}_{1}}{4{M}_{\star }}{b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\\ {\nu }_{\star ,i} & \equiv & \displaystyle \frac{3}{2}{n}_{i}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{i}}\right)}^{2}\\ {\nu }_{G,i} & \equiv & {n}_{i}\displaystyle \frac{3{m}_{G}}{4{M}_{\star }}{\left(\displaystyle \frac{{a}_{i}}{{a}_{G}}\right)}^{3}.\end{array}\end{eqnarray} \tag{ 4 }$$

These represent, respectively, the precession rate induced on planet 1 by planet 2, the precession rate induced on planet 2 by planet 1, the star-induced precession rate on planet i, and the giant-induced rates. The equations of motion may be solved in their most concise form by a canonical transformation to the complex variables (Morbidelli 2002; Batygin & Adams 2013),

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{1} & \equiv & \sqrt{\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}}\cos ({z}_{1})+i\sqrt{\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}}\sin ({z}_{1})\\ {\eta }_{2} & \equiv & \sqrt{\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}}}\cos ({z}_{2})+i\sqrt{\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}}}\sin ({z}_{2}),\end{array}\end{eqnarray} \tag{ 5 }$$

from which we may derive the inclinations

$$\begin{eqnarray}&&{I}_{i}=\sqrt{\displaystyle \frac{2{Z}_{i}}{{{\rm{\Lambda }}}_{i}}}=\sqrt{2{\eta }_{i}{\eta }_{i}^{* }}.\end{eqnarray} \tag{ 6 }$$

In terms of these new variables, Hamilton's equations read as (see, e.g., Batygin & Adams 2013)

$$\begin{eqnarray}&&{\dot{\eta }}_{i}=i\displaystyle \frac{\partial { \mathcal H }}{\partial {\eta }_{i}^{* }}\displaystyle \frac{1}{{{\rm{\Lambda }}}_{i}}.\end{eqnarray} \tag{ 7 }$$

Given the above definitions, the dynamical evolution of the system is governed by the linear differential equation,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\begin{array}{c}{\eta }_{1}\\ {\eta }_{2}\end{array}\right)={iM}\left(\begin{array}{c}{\eta }_{1}\\ {\eta }_{2}\end{array}\right)-i\displaystyle \frac{{\beta }_{\star }}{\sqrt{2}}\left(\begin{array}{c}{\nu }_{\star ,1}\\ {\nu }_{\star ,2}\end{array}\right),\end{eqnarray} \tag{ 8 }$$

where the matrix M takes the form

$$\begin{eqnarray}M=\left(\begin{array}{cc}{\nu }_{1}+{B}_{12} & -{B}_{12}\\ -{B}_{21} & {\nu }_{2}+{B}_{21}\end{array}\right),\end{eqnarray} \tag{ 9 }$$

and we define the sum of the precession rates induced by the star and giant as

$$\begin{eqnarray}&&{\nu }_{i}\equiv {\nu }_{\star ,i}+{\nu }_{G,i}.\end{eqnarray} \tag{ 10 }$$

Additional perturbers that are coplanar with the giant may be included by simply adding them to ν_i, such as the disk's potential, which we describe in Section 5.1.

The dynamical equations derived in the previous section are straightforward to generalize to n_p inner planets. Specifically, the matrix M becomes n_p × n_p in dimensionality, with elements

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{ii}} & = & {\nu }_{i}+\displaystyle \sum _{j=1,j\ne i}^{N}{B}_{{ij}}\\ {M}_{{ij}} & = & -{B}_{{ij}}.\end{array}\end{eqnarray} \tag{ 11 }$$

The frequencies B_ij are defined similarly to those in Equations (4),

$$\begin{eqnarray}&&{B}_{{ij}}\equiv {n}_{i}\displaystyle \frac{{m}_{j}}{4{M}_{\star }}{b}_{3/2}^{(1)}({\alpha }_{{ij}}){\alpha }_{{ij}}{\bar{\alpha }}_{{ij}},\end{eqnarray} \tag{ 12 }$$

where ${\bar{\alpha }}_{{jk}}=1$ if a_j > a_k and a_j/a_k otherwise. Finally, we define the complex inclination vector with elements ${\xi }_{i}\approx \sqrt{2}{\eta }_{i}$, the magnitude of which is equal to the orbital inclination I_i, leading to the generalized evolution equation

$$\begin{eqnarray}&&\boxed{\displaystyle \frac{d{\boldsymbol{\xi }}}{{dt}}=i{\boldsymbol{M}}{\boldsymbol{\xi }}-i{\beta }_{\star }{{\boldsymbol{\nu }}}_{\star }}.\end{eqnarray} \tag{ 13 }$$

Here we have introduced the vector ${{\boldsymbol{\nu }}}_{\star }$ with elements ${\nu }_{\star ,i}$ given by Equation (4). Physically, Equation (13) describes a system of n_p planetary-mass rings, perturbed by an exterior ring with zero inclination (the giant) and an inner ring of inclination β_⋆ (the stellar equatorial bulge).

The combined perturbations from the stellar quadrupole and the distant giant cause each of the close-in planets to possess an equilibrium tilt, the steady-state configuration of orbital inclinations. In the absence of a distant giant this equilibrium is star aligned, whereas in the absence of a J₂ the equilibrium is giant aligned. Solving Equation (13) for $\dot{{\boldsymbol{\xi }}}=0$ yields the equilibrium inclinations of each planet, given by

$$\begin{eqnarray}&&{{\boldsymbol{\xi }}}_{s}={\beta }_{\star }{{\boldsymbol{M}}}^{-1}{{\boldsymbol{\nu }}}_{\star }.\end{eqnarray} \tag{ 14 }$$

Note that neither the components of ${\boldsymbol{M}}$ nor ${{\boldsymbol{\nu }}}_{\star }$ depend on the stellar obliquity β_⋆ (Equation (4)). Accordingly, within the small-angle approximation, the equilibrium inclinations scale in direct proportion to the stellar obliquity.

The simplest possible system consists of one interior planet at a semimajor axis of a₁. (Later in Section 3.4, we will generalize to multiple interior planets.) In the single-planet case, the steady-state inclination coincides with the normal to the "Laplace Surface" (Tremaine et al. 2009). Solving the one-planet case of Equation (14) yields the equilibrium inclination

$$\begin{eqnarray}&&{\xi }_{s,1}=\displaystyle \frac{{\beta }_{\star }}{1+X},\end{eqnarray} \tag{ 15 }$$

where X is the ratio between the star-induced and giant-induced precession rates (as defined in Equation (4)),

$$\begin{eqnarray}&&X\equiv \displaystyle \frac{{\nu }_{G,1}}{{\nu }_{\star ,1}}\approx \displaystyle \frac{{m}_{G}{a}_{1}^{5}}{2{J}_{2}{R}_{\star }^{2}{a}_{G}^{3}{M}_{\star }}.\end{eqnarray} \tag{ 16 }$$

The Laplace surface may be visualized as the plane about which the inner planet's orbit precesses. Specifically, if the giant dominates, X ≫ 1 and ${\xi }_{s,1}\to 0$. In this case, the inner planet's orbital plane remains unaltered, assuming that it begins in a disk-aligned state (i.e., with I₁ = 0). On the other hand, if the star dominates, X ≪  1 and ${\xi }_{s,1}\to {\beta }_{\star }$. In this case, the inner planets that inherit the disk's original plane will have their orbital planes altered as they precess about the star's (tilted) spin axis.

Crucially, even if the giant and star exert comparable influences, i.e., X ∼ 1, the equilibrium inclination still differs from the natal disk's plane by ${\xi }_{s,1}\sim {\beta }_{\star }/2$. Accordingly, the role of the stellar quadrupole is to displace the equilibrium orientations of the inner planets away from disk aligned. In order to deduce typical values of ${\xi }_{s,1}/{\beta }_{\star }$, or equivalently X, we first estimate values of J₂ for young stars at the point of disk dispersal. In the next section, we compare the inferred quadrupole moments to observed examples of distant giant planets.

The magnitude of J₂ has not yet been measured for young stars; however, the precession rates of planets orbiting the rapidly rotating main-sequence stars WASP-33 (Watanabe et al. 2019), Kepler-13A, and HAT-P-7 (Masuda 2015) reveal associated values of J₂ ≈ 10⁻⁴. With a much slower rotation rate, the Sun's J₂ is lower, at ∼10⁻⁷ (Pireaux & Rozelot 2003).

Stellar models serve as the only available estimates of J₂ immediately after planet formation, when rapid rotation rates drive J₂ to larger values. Specifically, polytopic stellar structure models may be utilized to relate J₂ to the stellar spin rate and the tidal Love number k₂ (twice the apsidal motion constant⁶ ), through the following approximate relationship (Sterne 1939; Ward et al. 1976; Spalding & Batygin 2016):

$$\begin{eqnarray}\begin{array}{rcl}{J}_{2} & \approx & \displaystyle \frac{1}{3}{k}_{2}\displaystyle \frac{{\omega }_{\star }^{2}}{{{GM}}_{\star }/{R}_{\star }^{3}}\\ & \sim & {10}^{-3}\left(\displaystyle \frac{{k}_{2}}{0.2}\right){\left(\displaystyle \frac{{P}_{\star }}{\mathrm{day}}\right)}^{-2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{3}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 17 }$$

where the quantity $\sqrt{{{GM}}_{\star }/R{\star }^{3}}$ is the breakup rotational velocity. Inserting nominal parameters for WASP-33 yields J₂ ≈ 4 × 10⁻⁴ (Iorio 2011), slightly larger than that inferred empirically by Watanabe et al. (2019) but in agreement to an order of magnitude.

Pre–main-sequence stars are initially fully convective, corresponding to k₂ ≈ 0.2 (Batygin & Adams 2013). Eventually, those with M_⋆ ≳  0.3M_⊙ develop a radiative core, approaching k₂ ≈ 0.02 in the fully radiative limit. A Sun-like star likely remains fully convective throughout the disk-hosting stage, but more massive stars make the transition earlier on. We illustrate the equilibrium inclination of a single planet, scaled by the stellar obliquity in Figure 2 for three different semimajor axes and a range of J₂ values. We fix the distant giant to reside on a circular orbit of a_G = 2 au and possess a mass ${m}_{G}=2{M}_{J}$, which represent typical values (Bryan et al. 2019; Masuda et al. 2020).

Zoom In Zoom Out Reset image size

The stellar quadrupole is comparable to the giant's influence for a wide range of parameters. Specifically, a planet with a_p ≲ 0.05 au is forced to precess approximately around the stellar spin axis, even for the relatively small J₂ ∼ 10⁻⁵. Planets orbiting at 0.2 au feel a somewhat larger influence from the distant giant but nonetheless are driven toward inclinations of roughly half of the stellar obliquity at J₂ ∼ 10⁻³, values attainable during the first few million years of stellar evolution (Bouvier et al. 2014; Spalding & Batygin 2016).

Having discussed the typical properties of young stars, we now compare the stellar quadrupole to the distant giant orbital parameters required to provide a comparable secular influence. Specifically, Bryan et al. (2019) indicated that up to ∼30% of systems of close-in transiting planets are likely to possess an exterior giant with mass 0.5 < m_G/M_J < 20 and orbital distance $1\lt {a}_{G}/\mathrm{au}\lt 20$. Using nominal parameters of a_G = 2 au and m_G = 2 M_J, we find that

$$\begin{eqnarray}&&X\approx \frac{1}{20}{\left(\frac{{a}_{1}}{0.1\mathrm{au}}\right)}^{5}{\left(\frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-5}\left(\frac{{m}_{G}}{2{M}_{J}}\right){\left(\frac{{a}_{G}}{2\mathrm{au}}\right)}^{-3}{\left(\frac{{P}_{\star }}{\mathrm{day}}\right)}^{2}.\end{eqnarray} \tag{ 18 }$$

Thus, as illustrated in Figure 2, the stellar quadrupole dominates at early times for fiducial parameters. However, the strong dependence of X on ${({a}_{1}/{R}_{\star })}^{5}$ indicates that wider orbits may be more affected by the exterior giant (as indeed is the case in our solar system). Moreover, the tidal potential of the distant giant scales as m_G/a_G³, which ranges over many orders of magnitude (Bryan et al. 2016; Masuda et al. 2020).

In order to obtain a population-level idea of the relative influence of distant giants as compared to the stellar quadrupole, we turn to the NASA Exoplanet Archive (Akeson et al. 2013). We consider all confirmed, single-transiting planets possessing radii ${R}_{p}\lt 4{R}_{\oplus }$ and obtain their semimajor axes and host star radii, when available.

This procedure yielded 1211 systems. For each system, we compute the value of ${m}_{G}/{a}_{G}^{3}$ that would be required to generate X = 1 (Equation (16)). We consider two extreme cases for the stellar quadrupole. First, we considered a "strong quadrupole" case, where we assumed that the young star's period was P_⋆ = 1 day and its radius was twice the present-day value. The second case is the "weak quadrupole" regime, where we chose P_⋆ = 10 days and used the modern stellar radius.

In the left panel of Figure 3, we present histograms of the values of ${m}_{G}/{a}_{G}^{3}$ required to enforce X = 1 in the weak (blue) and strong (red) quadrupole cases. Superimposed, we illustrate the values of 10 giant planets, with well-constrained orbits, known to reside exterior to close-in transiting planets, as listed in Table 3 of Bryan et al. (2019). We included a red vertical line denoting Jupiter. Note that Jupiter's tidal potential is weaker than all 10 from the Bryan et al. (2019) sample. This is likely a result of biases, intrinsic to RV surveys, which favor giants that are closer-in and more massive. Accordingly, the typical influence of distant giants is likely lower than that inferred from the current observed sample.

Zoom In Zoom Out Reset image size

From inspection, the tidal potentials of the observed distant giants roughly coincide with the peak of the distribution corresponding to a weak stellar quadrupole (blue), but they fall far below the histogram illustrating a strong stellar quadrupole (red). These observations suggest that the secular influence of exterior giants is at best comparable to, but is often weaker than, the stellar quadrupolar influence.

Looking more specifically at the systems outlined by Bryan et al. (2019), we considered seven distant giant planets that possess inner super-Earths with well-characterized orbits.⁷ Of these close-in super-Earths, we considered the innermost member and computed X using the measured value of a₁, once again performing the calculation separately for a weak and strong stellar quadrupole. The right panel of Figure 3 illustrates the computed X values. As with the population-level histograms, X is typically smaller than unity unless the weaker limit for the stellar quadrupole is adopted, though possible exceptions are GJ 832 and Kepler-454.

In general, the observed distribution of distant giant parameters, coupled with the semimajor axes of close-in single-transiting planets, suggests that the stellar quadrupole exerts at least a comparable, but at most dominant, secular influence on close-in planetary systems. As mentioned above, this strong quadrupolar influence tends to force inner planetary orbital planes to precess about an axis that is displaced from disk aligned by a magnitude comparable to the stellar obliquity β_⋆. Close-in planets precess faster, leading to the excitation of mutual inclinations.

We now generalize the steady-state treatment from Section 3.1 to multiplanet systems by solving Equation (14) for n_p = 2 and 3. As in Figure 2, we choose three different values for ${a}_{1}=\{0.05,0.1,0.2\}$ au. For simplicity, we assume that the n_p close-in orbits are uniformly spaced with a distance of 20 mutual Hill radii, where we define a mutual Hill radius as

$$\begin{eqnarray}&&{R}_{H,\mathrm{mutual}}\equiv {\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{3{M}_{\star }}\right)}^{1/3}\displaystyle \frac{{a}_{1}+{a}_{2}}{2}.\end{eqnarray} \tag{ 19 }$$

This spacing is typical of Kepler compact multiplanet systems (Weiss et al. 2018). These planets also exhibit a high level of intrasystem uniformity in planetary radii (Weiss et al. 2018) and masses (Millholland et al. 2017); thus, we set all of the inner masses m_i equal. By assuming a constant separation in Hill radii, our calculations are relatively insensitive to the precise planetary mass; choosing a higher m_i would simply increase the separation between adjacent planets. Accordingly, we simply choose m_i = 5 M_⊕ throughout.

We illustrate the steady-state inclinations of n_p = 2 and 3 planet systems in Figure 4. By visual comparison to the case of n_p = 1 in Figure 2, the equilibrium inclinations of all n_p planets are similar to one another. Accordingly, in general, all planets within an inner system with a₁ ≲  0.2 au are initially forced to precess around a plane that differs significantly from that of the distant giant.

Zoom In Zoom Out Reset image size

Over time, the stellar quadrupole decays as a result of contraction, spin-down, and the formation of a radiative core. Indeed, the Sun possesses a value of J₂ ∼ 10⁻⁷ (Park et al. 2017), which is orders of magnitude lower than expected for pre–main-sequence stars. Figures 2 and 4 indicate that such small J₂ leads the inner planets to precess around the distant giant's plane at late times, regardless of the initial stellar J₂. Below, we show that stellar spin-down occurs over long-enough timescales that if the planets began in a star-aligned configuration, they adiabatically evolve to giant aligned in response to the loss of J₂.

We first suppose that the initial orientations of each planet are star aligned, given in terms of the complex inclination vector as

$$\begin{eqnarray}&&{\xi }_{i}{| }_{t=0}={\beta }_{\star },\end{eqnarray} \tag{ 22 }$$

i.e., each orbit normal is parallel to the stellar spin axis. For the sake of definiteness, we set the stellar obliquity to β_⋆ = 10° (similar to the Sun's 7°). We set τ_⋆ = 1 Myr and integrate Equation (8) for 10 Myr, or 10 e-folding times of the stellar quadrupole (which approximately amounts to a drop from J₂ = 10⁻² to J₂ < 10⁻⁶).

We integrate systems of n_p = 2 and 3 inner planets, choosing the innermost ${a}_{1}=0.1\,\mathrm{au}$, separating the planets by 20 mutual Hill spacings, as above. We present the time evolution of the planetary inclinations and their mutual inclinations in the left panels of Figure 5 (n_p = 2) and Figure 6 (n_p = 3). In the top left panels, colored solid lines depict the real part of the complex inclinations ξ_i, whereas colored dashed lines illustrate the equilibrium inclinations computed in the previous section. As the stellar quadrupole diminishes, the planets adiabatically reorient from star aligned to giant aligned. In the bottom left panels, we plot the mutual inclinations between all pairs of planets, which remain below 1°. Accordingly, in the star-aligned initial state, the inner planets remain coplanar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now suppose that the inner planets remain aligned with the disk's plane subsequent to disk dispersal, such that

$$\begin{eqnarray}&&{\xi }_{i}{| }_{t=0}=0.\end{eqnarray} \tag{ 23 }$$

As before, we simulate systems with n_p = 2 and 3 for 10 Myr and present the results in the right panels of Figures 5 and 6. In contrast to the left panels, a disk-aligned initial condition drives significant mutual inclinations between the inner planets (bottom right panels of Figures 5 and 6). Critically, throughout the subsequent evolution, these mutual inclinations remain high, exceeding the magnitude required to reduce the observed transit number (the horizontal dashed lines; Equation (21)). It is important to emphasize that the distant giant did not cause the mutual inclinations but still ended up inclined with respect to the inner planets.

Notice that as the stellar quadrupole decays, the planetary mutual inclinations decrease somewhat. This feature implies that while the initial, large J₂ drives mutual inclinations, later stellar spin-down plays a partial role in dynamically cooling close-in planetary systems. Similar behavior was highlighted in an analogous scenario by Anderson & Lai (2018), whereby stellar spin-down reduced the mutual inclination between a warm Jupiter and an exterior perturbing giant companion. Accordingly, in general, the mutual inclinations observed within modern-day, close-in planetary systems were likely once larger. The same may be said of satellite systems around giant planets that, like stars, lose their quadrupole moments with time (Batygin 2018). We defer a detailed discussion of this process to future work.

The primary findings of this section are summarized schematically in Figure 7. Specifically, an initial configuration of n_p planets inclined with respect to the stellar spin axis (but aligned with a distant giant) acquires substantial mutual inclinations. The converse scenario of initial alignment with the star leads to negligible mutual inclinations among close-in planets and alignment with the distant giant as J₂ decays. Accordingly, the picture described here is fully consistent with the inference by Masuda et al. (2020)—that giants orbiting exterior to single-transiting systems are inclined, whereas those exterior to multiple-transiting close-in systems are well aligned. The question is then, what determines whether star-aligned or giant-aligned initial conditions prevail subsequent to planet formation? In the next section, we will show that the disk dispersal timescale is the primary governing factor.

Zoom In Zoom Out Reset image size

In order to simulate inclination excitation during disk dispersal, we model the secular potential of the disk by way of a precession frequency ${\nu }_{d,i}$, given by (Hahn 2003)

$$\begin{eqnarray}&&{\nu }_{d,i}={n}_{i}\displaystyle \frac{\pi {\rm{\Sigma }}{a}_{i}^{2}}{{M}_{\star }\beta },\end{eqnarray} \tag{ 24 }$$

which is added to ν_i in Equation (10). We define the aspect ratio of the disk, β ≈ 0.05, and the surface density, Σ. This frequency is then used to compute the diagonal elements of the updated matrix ${\boldsymbol{M}}$.

The disk's surface density as a function of distance and time, ${\rm{\Sigma }}(a,t)$, is chosen to follow an infinite Mestel disk (Mestel 1963; Binney & Tremaine 2011; Schulz 2012):

$$\begin{eqnarray}&&{\rm{\Sigma }}(a,t)={{\rm{\Sigma }}}_{0}(t)\left(\displaystyle \frac{{a}_{0}}{a}\right).\end{eqnarray} \tag{ 25 }$$

The scaling factor Σ₀(t) is allowed to decay exponentially over a timescale τ_d following

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}(t)={{\rm{\Sigma }}}_{\mathrm{0,0}}\exp \left(-\displaystyle \frac{t}{{\tau }_{d}}\right).\end{eqnarray} \tag{ 26 }$$

We choose a₀ = 0.2 au and Σ_0,0 = 1300 g cm⁻², approximately in keeping with the minimum-mass extrasolar nebula (Chiang & Laughlin 2013). This is sufficient to enforce initial coplanarity with the disk plane.

Heuristically, disk dispersal should be adiabatic if τ_d is shorter than the average precession timescale of the inner n_p planets. This average rate roughly corresponds to the minimum eigenvalue λ_min of matrix M. In Figure 8, we illustrate typical values of ${\lambda }_{\min }^{-1}$ using systems of n_p = 3 inner planets with a₁ = 0.05 and 0.1 au. As can be seen, the timescale ${\lambda }_{\min }^{-1}$ can be as long as ∼10⁵ yr or as short as ∼10³ yr, depending on the giant's semimajor axis and the stellar J₂. Given the large J₂ ≳ 10⁻³ expected of young stars, we expect the system to evolve adiabatically if the disk's mass in the inner regions is removed over a timescale of τ_d ≳  10^3–4 yr. This is roughly in keeping with the viscous timescale of the inner 0.2 au of the disk gas, which we return to in the discussion (see Equation (31)).

Zoom In Zoom Out Reset image size

We test the dependence of planet–star inclinations on τ_d by performing a set of secular integrations that include the disk's secular potential. Specifically, we choose a system with a₁ = 0.1 au, n_p = 3 inner planets and an initial ${J}_{\mathrm{2,0}}={10}^{-3}$ (indicative of P_⋆ = 1 day; Equation (17)). These parameters correspond to an initial ${\lambda }_{\min }^{-1}\approx 3400$ yr. We simulate three different disk dispersal timescales: a rapid case where τ_d = 100 yr, a marginally adiabatic case with τ_d = 1000 yr, and a fully adiabatic case where τ_d = 10⁴ yr. The stellar obliquity is fixed at β_⋆ = 10° for each run, but note that in the Laplace–Lagrange approximation used here, all angles scale roughly linearly with β_⋆.

The evolution of the three inner planetary orbits during disk dispersal is presented in Figure 9. Orbital inclinations are presented in the top panels, with dashed lines tracking the equilibrium inclinations (Equation (14)). As the disk disperses, the equilibrium inclinations transition from disk aligned to a value of around 8°, close to the stellar obliquity of 10°. Mutual inclinations between all inner planet pairs are plotted in the bottom panels. The adiabatic case of a slowly dispersing disk is presented in the rightmost panel. Here, all three inner planets closely follow their forced equilibria, maintaining coplanarity. Thus, slow disk dispersal yields coplanar systems that are aligned with the distant giant.

Zoom In Zoom Out Reset image size

For shorter disk dispersal timescales (middle and left panels) the inner planets acquire their inclinations more impulsively. When the disk's dispersal is marginally adiabatic (τ_d = 1000 yr; middle panel), the planets track their equilibrium inclinations less exactly as compared to τ_d = 10⁴ yr, but they nevertheless maintain small mutual inclinations with one another. Only planets 1 and 3 develop mutual inclinations exceeding ∼2°, the typical rms spread of mutual inclinations in Kepler systems (Fabrycky et al. 2014). Dispersal must be strongly nonadiabatic in order to generate mutual inclinations exceeding half of the stellar obliquity of 10°, as illustrated by the case with τ_d = 100 yr.

When the stellar obliquity is set to β_⋆ = 10°, the most rapid disk dispersal timescale of 100 yr generates mutual inclinations that are only barely large enough to exceed the co-transiting criterion (horizontal dotted lines in Figure 9). We repeat the secular simulation above, this time with a larger stellar obliquity of 50°. The results are presented in Figure 10, where we illustrate the rapid disk dispersal time of τ_d = 100 yr. Moreover, we continue the simulation throughout the subsequent decay of the stellar quadrupole (with τ_⋆ = 0.25 Myr in this case). At the larger tilt of β_⋆ = 50°, mutual inclinations more often exceed the coplanarity criterion (alternatively, a larger ${J}_{\mathrm{2,0}}={10}^{-2}$ coupled with a smaller stellar obliquity suffices for generating nontransiting configurations; Figure 5).

Zoom In Zoom Out Reset image size

In general, our secular calculations show that stellar obliquities exceeding ${\beta }_{\star }\sim 10^\circ$, coupled with rapid disk dispersal, are able to reduce the transit multiplicities of planetary systems. Coplanarity is retained if the disk disperses slowly. Subsequently, the planets are finally drawn back to their original plane—aligned with the giant—as the stellar quadrupole decays.

Though the secular approach is critical for developing insight, it suffers several drawbacks. First, Laplace–Lagrange secular theory loses accuracy at high inclinations (Murray & Dermott 1999), such that our simulations with β_⋆ = 50° are unreliable at a quantitative level. Second, the secular frequencies of planetary systems vary with time as the disk disperses and J₂ changes (Nagasawa et al. 2005). This process can raise eccentricities and inclinations, or even drive instabilities (Ward 1981; Spalding et al. 2018). Lastly, we assumed that the giant's orbit and stellar spin axes are fixed in time, which is not correct in detail. Given these shortcomings, in the next section we check the validity of our analysis by turning to N-body simulations.

Stellar obliquities spanning the full range from 0° to 180° have been detected in a diverse collection of planetary systems (Winn et al. 2010, 2017; Albrecht et al. 2012; Huber et al. 2013; Dai & Winn 2017; Dalal et al. 2019). Large stellar obliquities were first widely observed among hot Jupiter hosts (Winn et al. 2010) but have since been detected in multiple-transiting systems (Huber et al. 2013) and in systems with planets down to Earth's size (Kamiaka et al. 2019). Therefore, nonzero stellar obliquities occur in many planetary systems at some point in their evolution; what is less certain is when these obliquities are excited.

Here, we have focused on scenarios whereby stellar obliquities are excited while the disk is still present, including gravitational perturbations from an exterior stellar companion (Batygin 2012; Lai 2014; Spalding & Batygin 2014; Zanazzi & Lai 2018), magnetic torques between the disk and star (Lai et al. 2011; Spalding & Batygin 2015), or simply side effects of star formation within a turbulent environment (Bate et al. 2010; Spalding et al. 2014; Fielding et al. 2015). However, it is currently unclear whether nonzero stellar obliquities typically emerge during these early times or through later, dynamical processes (Chatterjee et al. 2008; Naoz 2016; Ngo et al. 2016; Anderson & Lai 2018).

Unfortunately, most of the techniques leveraged to measure spin–orbit misalignments, such as the Rossiter-McLaughlin effect (Winn et al. 2005), gravity darkening (Barnes 2009), and asteroseismology (Huber et al. 2013), are rarely applicable to disk-hosting stars. An alternative approach is to photometrically constrain the stellar line broadening via $v\,\sin {i}_{\star }$ and divide by an estimated rotational velocity v to constrain the inclination i_⋆ (Winn et al. 2017). Among debris disks, this approach has suggested that stellar spin axes are usually aligned with the disk's plane to within a few tens of degrees (Watson et al. 2011; Greaves et al. 2013; Matthews et al. 2013).

More recently, Davies (2019) extended similar investigations to protoplanetary disks, while noting that large debris disks may be less common among systems possessing star–disk misalignments. Of 15 disks studied, 5 exhibited misalignments between the stellar spin axis and the disk's plane exceeding ∼10°. These measurements hint that primordial stellar obliquities are often small, but values of several tens of degrees are relatively abundant. Moreover, systems that are truly misaligned can appear more aligned depending on their sky-projected inclinations, potentially leading to an underestimation of their obliquities (Davies 2019).

Reliable measurements of star–disk misalignments remain sparse, leaving the relationship between stellar spin axes and their disks poorly constrained. Nevertheless, our work here generates numerous predictions that may be tested using observed planetary systems. Specifically, if no exterior giant is present, we suggest that in the case of a slowly dispersing disk, the inner planets adiabatically align with the stellar equator. This is consistent with the tendency thus far for multiplanet systems to exhibit low stellar obliquities (Sanchis-Ojeda et al. 2012; Winn et al. 2017; Wang et al. 2018; Hirano et al. 2020).

If, on the other hand, a distant giant is present during slow disk dispersal, the inner planets first adiabatically align with the stellar equator, before returning to the giant's plane during subsequent stellar spin-down. The outcome is coplanarity among the inner planets and distant giant, but a potential obliquity for the host star. Recent observations have shown that distant giants are typically aligned with close-in multiple-transiting systems (Masuda et al. 2020) but misaligned with transit singles.

The trends uncovered by Masuda et al. (2020) are consistent with both a giant-driven and stellar-oblateness-driven origin to mutual inclinations. Thus, stellar obliquity emerges as a key observational feature that may partially disentangle these two mechanisms. Specifically, consider a system of close-in planets, coplanar with a distant giant, but misaligned with the stellar spin axis. Such a configuration can only exist if either (1) the host star was never oblate enough to mutually incline the inner planets or (2) the disk dispersed slowly. Stellar evolution models may in principle be used to rule out the former option, leaving the modern-day stellar obliquity as an empirical signpost of slow disk dispersal in specific cases.

Unfortunately, measuring the stellar obliquity of close-in multiplanet systems that also possess a transiting exterior giant remains a significant observational challenge. The post–main-sequence star Kepler-56 possesses a giant companion of ${m}_{G}\,\sin \,{\rm{I}}=5.6{{\rm{M}}}_{J}$ and ${a}_{G}=2.2\,\mathrm{au}$ (corresponding to ${m}_{G}/{a}_{G}^{3}\approx 0.5{M}_{J}/{\mathrm{au}}^{3};$ Huber et al. 2013; Otor et al. 2016; see also Figure 3), orbiting exterior to two transiting planets that are inclined with respect to the stellar spin axis by over 40°. The exterior giant's orientation is unknown, but if shown to transit, it may support the mechanism described here.

In short, spin–orbit misalignments during the disk-hosting stage are expected on theoretical grounds (Batygin 2012; Lai 2014; Fielding et al. 2015; Spalding & Batygin 2015). Observational confirmation of this expectation has arisen for a limited number of cases (Davies 2019), but a substantially larger set of measurements is required in order to deduce a statistical distribution of star–disk misalignments.

In our analysis, we restricted attention to the case of a star with solar mass and radius. However, pre–main-sequence stars contract from large stellar radii of $\gt 2\,{R}_{\odot }$ to a final value that depends on the stellar mass. Moreover, this contraction occurs over a timescale comparable to the disk lifetime. Upon leaving their Hayashi tracks, stars with ${M}_{\star }\gtrsim 0.3\,{M}_{\odot }$ develop a radiative core, affecting k₂ (Gregory et al. 2016). Consequently, the stellar quadrupole weakens substantially during pre–main-sequence evolution.

In turn, disk lifetimes vary widely. Shorter-lived disks will disperse while the stellar quadrupole is stronger, as compared to longer-lived disks. More massive stars appear to lose their disks sooner (Ribas et al. 2015) and arrive onto the main sequence with larger radii (Gregory et al. 2016). Consequently, we expect a modest trend whereby higher-mass stars drive greater mutual inclinations among their close-in planets. A factor confounding this expectation is that a stronger stellar quadrupole augments the secular frequencies, such that faster disk dispersal timescales would be required to excite mutual inclinations. As discussed below, more information regarding disk dispersal timescales is required to delineate this expectation.

In contrast to radii, stellar spin rates are remarkably uniform across time and stellar type during the disk-hosting stage (Bouvier et al. 2014). However, over longer timescales, stars with M_⋆  ≳  1.2 M_⊙ tend to retain rapid rotations throughout their entire main-sequence lifetimes (Kraft 1967; Skumanich 1972). Combined with larger radii, these more massive stars exert a strong quadrupolar influence on close-in planets long after disk dispersal, preventing adiabatic alignment with an exterior giant. Accordingly, we expect stars with M_⋆  ≳  1.2 M_⊙ hosting multiplanet systems to exhibit smaller stellar obliquities, even if a distant giant is known to exist.

In previous investigations of star-driven misalignment, disk dispersal was assumed to be instantaneous (Spalding & Batygin 2016; Spalding et al. 2018; Li et al. 2020). Here, we have shown that this assumption is invalid unless the disk's gravity vanishes over a timescale comparable to or shorter than the system's slowest secular eigenmodes, typically between 10² and 10⁴ yr, with timescales below 10³ required for fiducial parameters (Figure 8). If disks typically disperse more slowly than this limit, then in general mutual inclinations among close-in planets are more likely to arise through later dynamical mechanisms that incline exterior giants (Gratia & Fabrycky 2017; Hansen 2017; Lai & Pu 2017; Pu & Lai 2018).

It is important to emphasize that the "disk dispersal timescale" mentioned here is distinct from the "disk lifetime." Disk lifetimes are relatively well constrained to lie between ∼1 and 10 Myr (Haisch et al. 2001; Mamajek 2009; Armitage 2011; Alexander et al. 2014; Ribas et al. 2015), with occasional longer-lived outliers (Silverberg et al. 2020). On the other hand, despite their significance, disk dispersal timescales are poorly constrained. Disks undoubtedly vanish over a timescale that is much shorter than their multiple million year lifetimes–a conclusion apparent from a lack of examples of disks that are in the process of dispersing (Cieza et al. 2008; Koepferl et al. 2013; Alexander et al. 2014). Moreover, the current inclinations of asteroids require that the inner several au of the solar nebula must have dispersed more rapidly than ∼10⁴–10⁵ yr (Ward 1981). Beyond these crude upper limits, little empirical data exist to constrain typical dispersal timescales.

Theoretical estimates of disk dispersal times are hampered by uncertainties regarding the dominant mechanism(s) driving their dispersal. A leading hypothesis is that UV photoevaporation of disk material at a few au drives a wind of material from the disk's surface (Alexander et al. 2014). As the disk ages and its accretion slows, the photoevaporative wind eventually exceeds the rate of inward disk-driven accretion. At this point, the inner few au of the disk is cut off from resupply and viscously drains onto the star. Modeling the photoevaporation process predicts that the inner 0.2 au can drain on timescales of ∼10 kyr (Alexander et al. 2006; Owen et al. 2010; Gorti et al. 2015).

A simple expression for the viscous time in terms of the dimensionless parameter α ∼ 10⁻³ may be written as (Shakura & Sunyaev 1973)

$$\begin{eqnarray}&&{\tau }_{\nu }\sim \displaystyle \frac{{a}^{2}}{\alpha {\rm{\Omega }}{h}^{2}}\approx 5700{\left(\displaystyle \frac{a}{0.2\mathrm{au}}\right)}^{3/2}\,\mathrm{yr}.\end{eqnarray} \tag{ 31 }$$

Though only a crude approximation, this timescale suggests that once starved from replenishment, the inner disk may disperse on timescales comparable to millennia. This is more rapid than the adiabatic threshold for small J₂ values (Figure 8) but may prevent mutual inclinations from being excited in systems with stronger J₂ and larger eigenvalues. Thus, paradoxically, larger stellar quadrupoles may inhibit mutual inclinations by increasing the eigenvalues and allowing planetary orbits to adiabatically reorient during disk dispersal.

In our simulations, we simply assumed that the disk's surface density decayed homogeneously and the disk remained coplanar throughout. However, this assumption is only valid if the inner disk mass is larger than that of the planets. We may solve for the disk mass required to dominate the quadrupolar potential felt by the planets by setting ${\nu }_{\star ,i}={\nu }_{d}$. This requirement yields a disk scaling factor of

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0,\mathrm{crit}}=\displaystyle \frac{3}{2}{J}_{2}\left(\displaystyle \frac{{a}_{1}}{{a}_{0}}\right){\left(\displaystyle \frac{{R}_{\star }}{{a}_{1}}\right)}^{2}\displaystyle \frac{\beta {M}_{\star }}{\pi {a}_{1}^{2}}.\end{eqnarray} \tag{ 32 }$$

The mass interior to radius a_out at this surface density may be computed as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{int}} & = & {\int }_{0}^{{a}_{\mathrm{out}}}2\pi a{{\rm{\Sigma }}}_{0,\mathrm{crit}}\left(\frac{{a}_{0}}{a}\right){da}\\ & = & 3{J}_{2}{\left(\frac{{R}_{\star }}{{a}_{1}}\right)}^{2}\left(\frac{{a}_{\mathrm{out}}}{{a}_{1}}\right)\beta {M}_{\star }\\ & \approx & \left(\frac{{J}_{2}}{{10}^{-3}}\right){\left(\frac{{a}_{1}}{0.1\mathrm{au}}\right)}^{-3}\left(\frac{{M}_{\star }}{{M}_{\odot }}\right)\left(\frac{{a}_{\mathrm{out}}}{\mathrm{au}}\right)\,{M}_{\oplus }.\end{array}\end{eqnarray} \tag{ 33 }$$

This result suggests that the disk's influence is comparable to the star's right down to M_int ∼ M_⊕. This is contradictory, since the planets in our simulations are 5 times that mass, violating the assumption that the disk remains coplanar. Moreover, the importance of disk dispersal time extends beyond the super-Earth regime considered here and into the regime of hot Jupiter inclinations (Zanazzi & Lai 2018), a mass regime long considered to carve out gaps in the natal disk (Goldreich & Tremaine 1980; Zhu et al. 2011). Accordingly, close-in planets undoubtedly contribute to the final disruption of the natal disk and may even hasten disk dispersal. The importance of planet–disk interactions for disk dispersal timescales is therefore an important outstanding problem that merits further study.