SPIN–ORBIT MISALIGNMENT AS A DRIVER OF THE KEPLER DICHOTOMY

In our solar system, the orbits of all eight confirmed planets are confined to the same plane with an rms inclination of ∼1°–2°, inspiring the notion that planets arise from protoplanetary disks (Kant 1755; Laplace 1796). By inference, one would expect extrasolar planetary systems to form with a similarly coplanar architecture. However, it is unknown whether such low mutual inclinations typically persist over billion-year timescales. Planetary systems are subject to many mechanisms capable of perturbing coplanar orbits out of alignment, including secular chaos (Laskar 1996; Lithkwick & Wu 2012), planet–planet scattering (Ford & Rasio 2008; Beaugé & Nesvorný 2012), and Kozai interactions (Naoz et al. 2011).

Despite numerous attempts, mutual inclinations between planets are notoriously difficult to measure directly (Winn & Fabrycky 2015). In light of this, investigations have turned to indirect methods. For example, by comparing the transit durations of co-transiting planets, Fabrycky et al. (2014) inferred generally low mutual inclinations (∼10–22) within closely packed Kepler systems. Additionally, within a subset of systems (e.g., 47 Uma and 55 Cnc) stability arguments have been used to limit mutual inclinations to ≲40° (Laughlin et al. 2002; Veras & Armitage 2004; Nelson et al. 2014). On the other hand, Dawson & Chiang (2014) have presented indirect evidence of unseen, inclined companions based upon peculiar apsidal alignments within known planetary orbits. Obtaining a better handle on the distribution of planetary orbital inclinations would lend vital clues to planet formation and evolution.

Recent work has attempted to place better constraints upon planet–planet inclinations at a population level, by comparing the number of single- to multi-transiting systems within the Kepler data set (Johansen et al. 2012; Ballard & Johnson 2016). Owing to the nature of the transit technique, an intrinsically multiple planet system might be observed as a single if the planetary orbits are mutually inclined. An emerging picture is that although a distribution of small ∼5° mutual inclinations can explain the relative numbers of double- and triple-transiting systems, a striking feature of the planetary census is a significant over-abundance of single-transiting systems. Furthermore, the singles generally possess larger radii (more with R_p ≳ 4 Earth radii), drawing further contrast.

The problem outlined above has been dubbed the "Kepler Dichotomy," and is interpreted as representing at least two separate populations; one with low mutual inclinations and another with large mutual inclinations that are observed as singles. The physical origin of this dichotomy remains unresolved (Morton & Winn 2014; Becker & Adams 2016). To this end, Johansen et al. (2012) proposed the explanation that planetary systems with higher masses undergo dynamical instability, leaving a separate population of larger, mutually inclined planets, detected as single transits. While qualitatively attractive, this model has two primary shortcomings. First, it cannot explain the excess of smaller single-transiting planets. Second, unreasonably high-mass planets are needed to induce instability within the required approximately gigayear timescales. Accordingly, the dichotomy's full explanation requires a mechanism applicable to a more general planetary mass range. In this paper, we propose such a mechanism—the torque arising from the quadrupole moment of a young, inclined star.

The past decade has seen a flurry of measurements of the obliquities, or spin–orbit misalignments, of planet-hosting stars (Winn et al. 2010; Albrecht et al. 2012; Huber et al. 2013; Morton & Winn 2014; Mazeh et al. 2015; Li & Winn 2016). A trend has emerged whereby hot stars (T_eff ≳ 6200 K) hosting hot Jupiters possess obliquities ranging from 0° to 180°, as opposed to their more modestly inclined, cooler (lower-mass) counterparts. Further investigation has revealed a similar trend among stars hosting lower-mass and multiple-transiting planets (Huber et al. 2013; Mazeh et al. 2015). Most relevant to the Kepler Dichotomy, Morton & Winn (2014) concluded at 95% confidence that single-transiting systems possess enhanced spin–orbit misalignment compared to multi-transiting systems.

Precisely when these spin–orbit misalignments arose in each system's evolution is still debated (Albrecht et al. 2012; Lai 2012; Storch et al. 2014; Spalding & Batygin 2016). However, the presence of stellar obliquities within currently coplanar, multi-planet systems hints at an origin during the disk-hosting stage (Huber et al. 2013; Mazeh et al. 2015). Indeed, many studies have demonstrated viable mechanisms for the production of disk–star misalignments, including turbulence within the protostellar core (Bate et al. 2010; Spalding et al. 2014; Fielding et al. 2015) and torques arising from stellar companions (Batygin 2012; Batygin & Adams 2013; Lai 2014; Spalding & Batygin 2014, 2015). Furthermore, Spalding & Batygin (2016) proposed that differences in magnetospheric topology between high and low-mass T Tauri stars (Gregory et al. 2012) may naturally account for the dependence of obliquities upon stellar (main sequence) T_eff. Crucially, if the star is inclined relative to its planetary system while young, fast-rotating, and expanded (Shu et al. 1987; Bouvier 2013), its quadrupole moment can be large enough to perturb a coplanar system of planets into a mutually inclined configuration after disk dissipation.

In what follows, we analyze this process quantitatively. First, we calculate the mutual inclination induced between two planets as a function of stellar oblateness (J₂), demonstrating a proof-of-concept that stellar obliquity suffices as a mechanism for over-producing single-transiting systems. Following this, we use N-body simulations to subject the famed, six-transiting system Kepler-11 to the quadrupole moment of a tilted, oblate star. We show that not only are the planetary orbits mutually inclined, but for nominal parameters the system itself can undergo a dynamical instability, losing three to five of its planets, with larger mass planets preferentially retained. In this way, we naturally account for the slightly larger observed size of singles (Johansen et al. 2012).

In order to motivate the following discussion, consider two planets orbiting in a shared plane around an inclined, oblate (high J₂) star. The effect of the stellar potential is to force a precession of each planetary orbit about the stellar spin pole, with the precession rate higher for the inner planet. If planet–planet coupling is negligible, the subsequent evolution would excite a mutual inclination between the planets of twice the stellar inclination (assuming fixed stellar orientation and negligible eccentricities). Alternatively, if planet–planet coupling is very strong, they will retain approximate coplanarity. Below, we analytically compute the system's evolution between these two extreme regimes (i.e., for general J₂).

2.1. Assumptions

We restrict our analytic calculation to small mutual inclinations between the planets and utilize Laplace–Lagrange secular theory (Murray & Dermott 1999). This framework assumes the planets to be far from mean motion resonances, allowing one to average over the orbital motion. Consequently, each planetary orbit becomes dynamically equivalent to a massive wire, a concept that is due to Gauss (Murray & Dermott 1999; Morbidelli 2002). Furthermore, we set all eccentricities to zero.¹

The star's orientation will be held fixed. The validity of this assumption can be demonstrated by considering the ratio of stellar spin to planetary orbital angular momenta:

$$\begin{eqnarray}&&\displaystyle \frac{{J}_{\star }}{{{\rm{\Lambda }}}_{{\rm{p}}}}\equiv \displaystyle \frac{{I}_{\star }{M}_{\star }{R}_{\star }^{2}{{\rm{\Omega }}}_{\star }}{{m}_{p}\sqrt{{{GM}}_{\star }{a}_{p}}}\end{eqnarray} \tag{ 1 }$$

where I_⋆ ≈ 0.21 is the dimensionless moment of inertia appropriate for a fully convective, polytropic star (Chandrasekar 1939), and the stellar rotation rate is Ω_⋆ = 2π/P_⋆. Consider a young, Sun-like star, possessing a rotation period of P_⋆ = 10 days (on the slower end of observations; Bouvier 2013) and a radius of roughly 2R_⊙ (Shu et al. 1987). A 10 Earth-mass object would need to orbit at over ∼100 au in order to possess the angular momentum of the star. Thus, provided we deal with compact, relatively low-mass systems, the stellar orientation can be safely fixed to zero.

A further assumption is that the dynamical influence of stellar oblateness may be approximated using only the leading order quadrupole terms, neglecting those of order ${ \mathcal O }({J}_{2}^{2})$. Therefore, the disturbing part of the stellar potential (with e = 0) may be written as (Danby 1992)

$$\begin{eqnarray}{ \mathcal R } & = & \displaystyle \frac{{{Gm}}_{{\rm{p}}}{M}_{\star }}{2{a}_{p}}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{p}}\right)}^{2}{J}_{2}\left(\displaystyle \frac{3}{2}{\sin }^{2}{i}_{p}-1\right)\\ & \approx & \displaystyle \frac{{{Gm}}_{{\rm{p}}}{M}_{\star }}{2{a}_{p}}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{p}}\right)}^{2}(6{s}_{p}^{2}-1),\end{eqnarray} \tag{ 2 }$$

where the second step has made the assumption of small planetary inclination i_p and defined a new variable ${s}_{p}\equiv \sin ({i}_{p}/2)$ (this definition is introduced to maintain coherence with traditional notation in celestial mechanics, e.g., Murray & Dermott 1999).

Finally, it is essential to define an initial configuration for the planetary system. Both numerical and analytic modeling of planet–disk interactions suggest that embedded protoplanets have their inclinations and eccentricities damped to small values within the disk-hosting stage (Tanaka & Ward 2004; Cresswell et al. 2007; Kley & Nelson 2012). Furthermore, any warping of the disk in response to a stellar companion is expected to be small (Fragner & Nelson 2010; Batygin et al. 2011). Therefore, throughout this work, we assume that the planets emerge from the disk with circular, coplanar orbits that are inclined by some angle β_⋆ relative to the star. Note that we will always fix the stellar spin direction to be aligned with the z axis, so β_⋆, the stellar obliquity, constitutes the initial inclination of the planetary orbits in our chosen frame.

2.2. Two-planet System

Incorporating the above assumptions, we may now write down the Hamiltonian (${ \mathcal H }$) that governs the dynamical evolution of the planetary orbits. To second order in inclinations (and dropping constant terms), we have (Murray & Dermott 1999)

$$\begin{eqnarray}&&{ \mathcal H }=\mathop{\underbrace{\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{a}_{2}}[({s}_{1}^{2}+{s}_{2}^{2}){f}_{3}+{s}_{1}{s}_{2}{f}_{14}\cos ({{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2})]}}\limits_{{\rm{Planet\mbox{--}planet}}\,{\rm{interaction}}}\\ &&\quad \,\,\mathop{\underbrace{-\displaystyle \frac{3{{Gm}}_{1}{M}_{\star }}{{a}_{1}}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{1}}\right)}^{2}{s}_{1}^{2}-\displaystyle \frac{3{{Gm}}_{2}{M}_{\star }}{{a}_{2}}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{2}}\right)}^{2}{s}_{2}^{2}}}\limits_{{\rm{Planet\mbox{--}stellar}}\,{\rm{quadrupole}}\,{\rm{interaction}}},\end{eqnarray} \tag{ 3 }$$

where the prefactors are

$$\begin{eqnarray}&&{f}_{3}=-\displaystyle \frac{1}{2}{f}_{14}=-\displaystyle \frac{1}{2}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right){b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right),\end{eqnarray} \tag{ 4 }$$

and ${b}_{3/2}^{(1)}$ is the Laplace coefficient

$$\begin{eqnarray}&&{b}_{3/2}^{(1)}(\alpha )\equiv \displaystyle \frac{1}{\pi }{\int }_{0}^{2\pi }\left[\displaystyle \frac{\cos \psi }{{(1+{\alpha }^{2}-2\alpha \cos \psi )}^{3/2}}\right]d\psi .\end{eqnarray} \tag{ 5 }$$

Because we are using Hamiltonian mechanics, the dynamics must be described in terms of canonical variables. Traditional Keplerian orbital elements do not constitute a canonical set, so we transform to Poincaré (or, modified Delauney; Murray & Dermott 1999) variables, defined as

$$\begin{eqnarray}&&{Z}_{{\rm{p}}}\equiv {m}_{{\rm{p}}}\sqrt{{{GM}}_{\star }{a}_{{\rm{p}}}}(1-\cos ({i}_{{\rm{p}}}))\quad {z}_{{\rm{p}}}\equiv -{{\rm{\Omega }}}_{{\rm{p}}}.\end{eqnarray} \tag{ 6 }$$

Physically, Z_p is the angular momentum of a circular orbit after subtracting its component in the z-direction. Notice that in the small angle limit,

$$\begin{eqnarray}&&{Z}_{{\rm{p}}}\approx \displaystyle \frac{1}{2}{m}_{{\rm{p}}}\sqrt{{{GM}}_{\star }{a}_{{\rm{p}}}}{i}_{{\rm{p}}}^{2}\equiv \displaystyle \frac{1}{2}{{\rm{\Lambda }}}_{{\rm{p}}}{i}_{{\rm{p}}}^{2}.\end{eqnarray} \tag{ 7 }$$

After substituting, we arrive at the governing Hamiltonian

$$\begin{eqnarray}{ \mathcal H } & = & -\,{ \mathcal E }\left[\displaystyle \frac{{Z}_{1}}{{{\rm{\Lambda }}}_{1}}+\displaystyle \frac{{Z}_{2}}{{{\rm{\Lambda }}}_{2}}-2\sqrt{\displaystyle \frac{{Z}_{1}{Z}_{2}}{{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}}}\cos ({z}_{1}-{z}_{2})\right]\\ & & -\displaystyle \frac{3}{2}{n}_{1}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{1}}\right)}^{2}{Z}_{1}-\displaystyle \frac{3}{2}{n}_{2}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{2}}\right)}^{2}{Z}_{2}\end{eqnarray} \tag{ 8 }$$

where, for compactness, we define

$$\begin{eqnarray}&&{ \mathcal E }\equiv \displaystyle \frac{{{Gm}}_{1}{m}_{2}}{4{a}_{2}}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right){b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right).\end{eqnarray} \tag{ 9 }$$

In order to complete the calculation, we define a complex variable that represents the inclination of each planet

$$\begin{eqnarray}{\eta }_{{\rm{p}}} & \equiv & \sqrt{\displaystyle \frac{{Z}_{{\rm{p}}}}{{{\rm{\Lambda }}}_{{\rm{p}}}}}(\cos ({z}_{{\rm{p}}})+{\imath }\sin ({z}_{{\rm{p}}}))\\ & \approx & \displaystyle \frac{1}{\sqrt{2}}{i}_{{\rm{p}}}(\cos ({{\rm{\Omega }}}_{{\rm{p}}})-{\imath }\sin ({{\rm{\Omega }}}_{{\rm{p}}})),\end{eqnarray} \tag{ 10 }$$

where ${\imath }=\sqrt{-1}$. The purpose is to cast Hamilton's equations into an eigenvector/eigenvalue problem. Specifically, in terms of these new variables, we must solve

$$\begin{eqnarray}&&{\dot{\eta }}_{{\rm{p}}}={\imath }\displaystyle \frac{\partial { \mathcal H }}{\partial {\eta }_{{\rm{p}}}^{* }}\displaystyle \frac{1}{{{\rm{\Lambda }}}_{{\rm{p}}}},\end{eqnarray} \tag{ 11 }$$

in which "*" denotes complex conjugation, yielding the matrix equation

$$\begin{eqnarray*}&&\displaystyle \frac{d}{{dt}}\left(\begin{array}{c}{\eta }_{1}\\ {\eta }_{2}\end{array}\right)=-{\imath }\left(\begin{array}{cc}{B}_{1}+{\nu }_{1} & -{B}_{1}\\ -{B}_{2} & {B}_{2}+{\nu }_{2}\end{array}\right)\left(\begin{array}{c}{\eta }_{1}\\ {\eta }_{2}\end{array}\right),\end{eqnarray*}$$

where we have defined four frequencies as

$$\begin{eqnarray}{B}_{1} & \equiv & \displaystyle \frac{1}{4}{n}_{1}{\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)}^{2}{b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\displaystyle \frac{{m}_{2}}{{M}_{\star }}\\ {B}_{2} & \equiv & \displaystyle \frac{1}{4}{n}_{2}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right){b}_{3/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\displaystyle \frac{{m}_{1}}{{M}_{\star }}\\ {\nu }_{{\rm{p}}} & = & \displaystyle \frac{3}{2}{n}_{{\rm{p}}}{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{{\rm{p}}}}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

The equation above may be solved using standard methods, whereby the solution is written as a sum of eigenmodes

$$\begin{eqnarray}&&{\eta }_{{\rm{p}}}=\displaystyle \sum _{j=1}^{2}{\eta }_{{\rm{p}},j}\exp ({\imath }{\lambda }_{j}t).\end{eqnarray} \tag{ 13 }$$

Indeed, the problem may be easily extended to N planets, though writing down all eigenvectors η_p,j and eigenmodes λ_j rapidly becomes cumbersome.

2.2.1. Initial Conditions and Solution

A choice must be made for the initial conditions of the problem. As already mentioned above, here we choose the condition that both orbits are initially coplanar, having recently emerged from their natal disk, with the star inclined by some angle β_⋆ relative to them. Accordingly, all four boundary conditions may be satisfied by requiring that

$$\begin{eqnarray}&&{\eta }_{{\rm{p}}}{| }_{t=0}=\displaystyle \frac{{\beta }_{\star }}{\sqrt{2}}.\end{eqnarray} \tag{ 14 }$$

What we seek is the mutual planet–planet inclination (denoted β_rel). In the small angle approximation, we can compute this quantity using the relation

$$\begin{eqnarray}(1-\cos ({\beta }_{{\rm{rel}}})) & \approx & {\eta }_{1}{\eta }_{1}^{* }+{\eta }_{2}{\eta }_{2}^{* }-({\eta }_{1}{\eta }_{2}^{* }+{\eta }_{2}{\eta }_{1}^{* })\\ & \approx & \displaystyle \frac{1}{2}{\beta }_{{\rm{rel}}}^{2}.\end{eqnarray} \tag{ 15 }$$

After solving for eigenvalues, eigenvectors, and matching the initial conditions, we arrive at the solution for the mutual inclination of the two planets as a function of time, which takes the rather compact form

$$\begin{eqnarray}&&\boxed{{\beta }_{{\rm{rel}}}(t)=2{\beta }_{\star }{ \mathcal G }({J}_{2})\sin ({\omega }_{0}t/2)},\end{eqnarray} \tag{ 16 }$$

where we define the (semi-)amplitude of the oscillations between planets

$$\begin{eqnarray}&&{ \mathcal G }={ \mathcal L }{\left[1+{{ \mathcal L }}^{2}+2\left(\displaystyle \frac{{{\rm{\Lambda }}}_{2}-{{\rm{\Lambda }}}_{1}}{{{\rm{\Lambda }}}_{2}+{{\rm{\Lambda }}}_{1}}\right){ \mathcal L }\right]}^{-1/2}\end{eqnarray} \tag{ 17 }$$

in terms of the ratio of frequencies

$$\begin{eqnarray}{ \mathcal L } & \equiv & \displaystyle \frac{{\nu }_{1}-{\nu }_{2}}{{B}_{1}+{B}_{2}}\\ & = & 6{J}_{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{1}}\right)}^{2}{\left(\displaystyle \frac{{a}_{2}}{{a}_{1}}\right)}^{2}\displaystyle \frac{1}{{b}_{3/2}^{(1)}(\alpha )}\displaystyle \frac{{M}_{\star }}{{m}_{2}}\displaystyle \frac{1-{\alpha }^{7/2}}{1+{{\rm{\Lambda }}}_{1}/{{\rm{\Lambda }}}_{2}}.\end{eqnarray} \tag{ 18 }$$

A convenient consequence of the aligned initial conditions is that the oscillations are purely sinusoidal, evolving with half the frequency

$$\begin{eqnarray}{\omega }_{0} & \equiv & [{({B}_{1}+{B}_{2})}^{2}+({\nu }_{1}-{\nu }_{2})\\ & & {\times ({\nu }_{1}-{\nu }_{2}+2({B}_{1}-{B}_{2}))]}^{1/2}.\end{eqnarray} \tag{ 19 }$$

One aspect to notice is that the amplitude is maximized when the equality

$$\begin{eqnarray}&&{\nu }_{1}+{{ \mathcal K }}^{2}{B}_{1}={\nu }_{2}+{{ \mathcal K }}^{2}{B}_{2}\end{eqnarray} \tag{ 20 }$$

is satisfied, where ${ \mathcal K }\equiv ({{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2})/({{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2})$. The maximum amplitude, scaled by β_⋆, can then be written as

$$\begin{eqnarray}&&2{{ \mathcal G }}_{{\rm{\max }}}=\displaystyle \frac{{{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2}}{\sqrt{{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}}}.\end{eqnarray} \tag{ 21 }$$

The significance of this result is best seen upon considering the outer planet to be a test particle, such that B₁ = 0 and ${{\rm{\Lambda }}}_{2}\to 0$. In such a scenario, the maximum of the amplitude ${{ \mathcal G }}_{{\rm{\max }}}\to \infty$. Such an unphysical result occurs as a consequence of a secular resonance (Murray & Dermott 1999; Morbidelli 2002; Spalding & Batygin 2014; Batygin et al. 2016), whereby the inner and outer bodies precess at similar rates. In reality, our earlier approximation that mutual inclinations are small breaks down in this regime and the inclusion of higher order terms is required.

In principle, one may also set ${{\rm{\Lambda }}}_{1}\to 0$ and conclude that the above resonance persists when the inner planet is a test particle. However, the resonant criterion in terms of stellar oblateness reads

$$\begin{eqnarray}{J}_{2}{| }_{{\rm{res}}} & \approx & \displaystyle \frac{1}{6}{\left(\displaystyle \frac{{a}_{1}}{{R}_{\star }}\right)}^{2}{\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)}^{2}\left(\displaystyle \frac{{m}_{2}}{{M}_{\star }}\right)\\ & & \times \left(\displaystyle \frac{{b}_{3/2}^{(1)}(\alpha )}{{{\rm{\Lambda }}}_{2}}\right)\displaystyle \frac{{({{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2})}^{2}}{({{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2})(1-{\alpha }^{7/2})},\end{eqnarray} \tag{ 22 }$$

which is negative when Λ₁ < Λ₂. Accordingly, the condition for secular resonance can only be satisfied when Λ₁ > Λ₂, i.e., when the inner planet possess more orbital angular momentum than the outer planet.

As an illustration, we plot the semi-amplitude ${ \mathcal G }$ in Figure 1 appropriate for a configuration where the two planetary orbits are situated at 0.05 and 0.1 au, both orbiting a solar-mass star with radius 0.01 au (about twice the Sun's radius). Three cases are shown: the red line depicts a 1 Earth-mass planet interior to a 10 Earth-mass planet while the blue line has the planets interchanged. The former configuration possesses a positive ${J}_{2}{| }_{{\rm{res}}}$, appearing as a maximum in the amplitude. The third situation (the cyan line) represents a 100 Earth-mass planet interior to a 1 Earth-mass body, illustrating that higher-mass inner planets may more easily misalign their outer companion (Equation (21)), but ${J}_{2}{| }_{{\rm{res}}}$ is correspondingly higher.

Zoom In Zoom Out Reset image size

Figure 1 demonstrates that misalignments of the order twice the stellar obliquity can be readily excited for reasonable values of J₂. By geometric arguments, the potential for such misalignments to take one of the planets out of transit depends upon the ratio R_⋆/a. However, for the cases considered above, only ∼4° of stellar obliquity are required to remove the two planets from a co-transiting configuration (less than the ∼7° present in the solar system; Lissauer et al. 2011). Conversely, planets may remain co-transiting if the innermost planet is sufficiently distant, the planets are very massive and/or tightly packed, or the stellar quadrupole moment is particularly low.

Several crucial aspects of real systems were neglected in order to obtain the analytic solution (16) above. Principally, we included only two planets whose orbits were assumed to be circular and only slightly mutually inclined. Additionally, in averaging over short-term motion, our adopted secular approach is unable to describe the full dynamics. A more subtle aspect was that we considered a constant J₂, when in reality, stars are expected to spin down and shrink over time until J₂ is essentially negligible (Irwin et al. 2008; Bouvier 2013; McQuillan et al. 2013).

In order to test our hypothesis within a more general framework, we now turn to N-body simulations. To carry out the calculations, we employed the well-tested Mercury6 symplectic integration software package (Chambers 1999).² In addition to standard planet–planet interactions, we modified the code to include the gravitational potential arising from a tilted star of given J₂, along with a term to produce general relativistic precession (following Nobili & Roxburgh 1986).

For the sake of definiteness, the parameters of our modeled system were based off of Kepler-11, a star around which six transiting planets have been discovered (Lissauer et al. 2011). Detailed follow-up studies, using Transit Timing Variations, have constrained the masses of the innermost five planets and placed upper limits upon the mass of Kepler-11g, the outermost planet,³ making this system ideal for dynamical investigation. Though choosing one system is not exhaustive, our goal is to demonstrate the influence of a tilted star upon a general coplanar system of planets. We follow Lissauer et al. (2013) and use their best-fit mass of 8 Earth Masses for Kepler-11g, with the stellar mass given by 0.961M_⊙ (see Table 1).

3.1. N-body Simulation

For our numerical runs, we choose 10 values of stellar obliquity and 11 of initial stellar J₂ = J_2,0 (i.e., the oblateness immediately as the disk dissipates). Once again, we fix the stellar orientation aligned with the z axis, but choose the initial planet–star misalignments:

$$\begin{eqnarray}&&{\beta }_{\star }\in \{5,10,20,30,40,50,60,70,80,85\}.\end{eqnarray} \tag{ 23 }$$

The value of J₂ for a star deformed by its own rotation may be related to its spin rate ω_⋆ and Love number (twice the apsidal motion constant) k₂ as follows (Sterne 1939):

$$\begin{eqnarray}&&{J}_{2}=\displaystyle \frac{1}{3}{\left(\displaystyle \frac{{\rm{\Omega }}}{{{\rm{\Omega }}}_{b}}\right)}^{2}\,{k}_{2},\end{eqnarray} \tag{ 24 }$$

where Ω_b is the break-up spin frequency at the relevant epoch. The Love number can be estimated by modeling the star as a polytope with index χ = 3/2 (i.e., fully convective; Chandrasekar 1939), which yields k₂ ≈ 0.28. Observations constrain the spin-periods of T-Tauri stars to lie within the range of ∼1–10 days (Bouvier 2013), while the break-up period is given by

$$\begin{eqnarray}&&{T}_{{\rm{b}}}=\displaystyle \frac{2\pi }{{{\rm{\Omega }}}_{b}}\approx \displaystyle \frac{1}{3}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{\star }}{2{R}_{\odot }}\right)}^{3/2}{\rm{days}}.\end{eqnarray} \tag{ 25 }$$

In our simulations below, we use the current mass of Kepler-11 for the star, but suppose its radius to be somewhat inflated relative to its current radius (R_⋆ = 2R_⊙), reflecting the T Tauri stage (Shu et al. 1987). With these parameters, we arrive at a reasonable range of J_2,0 of

$$\begin{eqnarray}&&{10}^{-4}\lesssim {J}_{\mathrm{2,0}}\lesssim {10}^{-2},\end{eqnarray} \tag{ 26 }$$

within which we choose 11 values uniformly separated in log-space:

$$\begin{eqnarray}&&{J}_{\mathrm{2,0}}\in \{{10}^{-4},{10}^{-3.8}\,\ldots \,{10}^{-2}\}.\end{eqnarray} \tag{ 27 }$$

In all runs, rather than allowing both R_⋆ and J₂ to vary, we simply left R_⋆ as a constant, letting J₂ decay exponentially over a timescale of τ = 1 Myr

$$\begin{eqnarray}&&{J}_{2}(t)={J}_{\mathrm{2,0}}\,{{\rm{e}}}^{(-t/\tau )}.\end{eqnarray} \tag{ 28 }$$

The choice for τ is essentially arbitrary, provided J₂ decays over many precessional timescales, owing to the adiabatic nature of the dynamics (Lichtenberg & Lieberman 1992; Morbidelli 2002). Our choice of 1 Myr roughly coincides with a Kelvin–Helmholtz timescale (Batygin & Adams 2013) but is chosen also to save computational time.

For each case, our integrations span 22 million years, beginning with the initial condition of a coplanar system possessing the current semi-major axes of the Kepler-11 system, but with eccentricities set to zero (Lissauer et al. 2013). In order to analyze the results, we sample the system six times between 19 and 22 Myr, and at each step calculate the maximum number of transiting planets that could be observed from a single direction. The results at all six times were then averaged.

The determination of the maximum number of transits was accomplished as follows. We begin by checking whether all possible pairs within the six planets mutually transit, where the criterion for concluding a pair of planets to be non-transiting is

$$\begin{eqnarray}&&| \sin ({\beta }_{{\rm{rel}}})| \gt \sin ({\beta }_{{\rm{crit}}})\approx \displaystyle \frac{{R}_{\star }}{{a}_{1}}+\displaystyle \frac{{R}_{\star }}{{a}_{2}},\end{eqnarray} \tag{ 29 }$$

where in the above formula we used the current radius of Kepler-11 (R_⋆ = 1.065R_⊙, as opposed to the inflated value relevant to the T Tauri stage). If any single pair of planets was non-transiting, we proceeded to choose each possible combination of five out of the six and performed a similar pairwise test to identify potentially observable five-transiting systems. If no set of five passed the test, we chose all sets of four, etc., until potentially finding that only one planet could be seen transiting. We note that the criterion above neglects the possibility of fortunate orbital configurations allowing two mutually inclined orbits to intersect along the line of sight. This complication, however, does not affect our qualitative picture.

The numerical results are presented in Figure 2, where the x-axis depicts the initial stellar J₂ = J_2,0. The y-axis refers to the misalignment β_⋆ between the stellar spin axis and the initial plane of the six planets. Each run has been given its own rectangular box, within which the color represents the maximum number of co-transiting planets that may be observed around the star (as discussed above). The numerics verify our analytic result, in that the observable multiplicity may be significantly reduced solely as a consequence of stellar obliquity.

Zoom In Zoom Out Reset image size

As expected, higher values of J_2,0 result in fewer transiting planets, provided the star is tilted relative to the planetary orbits. As with our analytic results, even small stellar obliquities are sufficient to reduce the transit count, with 5° of obliquity reducing the transit number to as little as three (Figure 2). However, planet–planet mutual inclination was not the only source of the reduction in transit number. A crucial finding was that for large enough J_2,0 and β_⋆, the stellar quadrupole potential caused the system to go unstable, casting three to five planets out of the system or into the central body, with planet–planet collisions existing as an additional possibility not captured in our simulations (Boley et al. 2016).

The region of instability (i.e., where at least one planet was lost) is outlined by a dotted line in Figure 2. Interestingly, the areas of instability map closely onto the regions where only a single transit remains. In other words, almost every single-transiting system coming out of the integration had lost planets through dynamical instability. Furthermore, each time instability occurred, the two lowest mass members were lost: Kepler-11b and f, with the next lowest mass body, Kepler-11c often joining them. Such a preference for retaining more massive planets is indeed reflected in the data as a slightly larger typical radius for single-transiting systems (Johansen et al. 2012). However, more testing is required to determine whether this is a generic feature of our model.

We showed that for any given two-planet system, there exists a resonant J₂ if the inner planet has more angular momentum than the outer. However, the picture becomes much more complicated in a multi-planet system, where each planet introduces two additional secular modes (one for eccentricity and one for inclination; Murray & Dermott 1999), increasing the density of resonances in Fourier space. As the stellar J₂ decays, its influence sweeps across each resonance, providing ample opportunity to excite mutual inclinations. If the two planets, Kepler-11d and f, were alone, they could become resonant at ${J}_{2}{| }_{{\rm{res}}}\sim {10}^{-2.4}$, which coincides approximately with the onset of instability in the low-inclination runs (Figure 2) but not exactly, for the reasons mentioned above.

In recent years, the Kepler data set has grown sufficiently comprehensive to facilitate statistically robust investigations at a population level. Out of this data has emerged a so-called "Kepler Dichotomy"; the notion that single-transiting systems are too common to be explained as resulting from a simple distribution of mutual inclinations within systems of higher multiplicity (Johansen et al. 2012; Ballard & Johnson 2016).

In separate studies, significant misalignments have been detected between the stellar spin axes and the planetary orbits they host, particularly around stars with main sequence T_eff ≳ 6200 K (Winn et al. 2010; Albrecht et al. 2012; Mazeh et al. 2015; Li & Winn 2016). This trend initially became apparent within the hot-Jupiter data set and was consequently often interpreted as evidence for a post-disk, high-eccentricity migration pathway for hot-Jupiter formation.⁴ However, a similar trend has now emerged within Kepler systems, including the multi-transiting sub-population (Huber et al. 2013; Mazeh et al. 2015), with little evidence supporting a tidal origin (Li & Winn 2016). These observations cumulatively suggest that many of the misalignments originated from directly tilting the protoplanetary disk, thereby inclining all planets in the system at once (Batygin 2012; Lai 2014; Spalding & Batygin 2014, 2015; Fielding et al. 2015).

A consequence of primordially generated spin–orbit misalignments is that stellar obliquity would be present at the end of the natal disk's life, leaving the planetary orbital architecture subject to the quadrupole moment of their young, rapidly spinning, and expanded host star. This paper has demonstrated that such a configuration naturally misaligns close-in systems and, furthermore, provides a mechanism for dynamical instability that by-passes the problem encountered in earlier work that unreasonably large planets were required to induce instability (Johansen et al. 2012).

The observable multiplicity of transiting systems can be reduced either by inclining planetary orbits relative to each other, or by intrinsically reducing the number of planets. Here, we have shown that both can be at play, with modest J₂ and stellar obliquity causing misalignments, whereas sufficiently large values thereof lead to dynamical instability, shedding planets. The origin of the instability is likely secular in nature, and significant planet–planet inclinations have been shown to reduce the inherent stability of planetary systems in numerous previous works (Laughlin et al. 2002; Veras & Armitage 2004; Nelson et al. 2014). In support of this interpretation, our simulations resulted in planetary instability at much smaller J₂ when obliquity was high ≳40°. Accordingly, we would expect multiplicity (both transiting and intrinsic) to be lower around hot stars, which tend to possess higher obliquities (Winn et al. 2010; see below).

5.1. Predictions

5.2. Kepler Data

To obtain data on confirmed, Kepler systems, we downloaded the data from the "Confirmed Planets" list (as of 2016 July) on the NASA Exoplanet Archive website.⁵ The systems were filtered to include only those in the Kepler field, though the conclusions that follow do not change qualitatively if non-Kepler detections are included.

In Figure 3, we split the data into "hot" stars with T_eff > 6200 K (132 in total) and"cool" stars, with T_eff < 6200 K (1504 in total). For each sub-population, we illustrate the fraction of systems as a function of the number of planets observed in transit. The hot stars clearly demonstrate a larger fraction of singles and a smaller fraction of multiples for each value of multiplicity, in agreement with the predictions of our model. In order to quantify the significance of this agreement, we carry out a statistical test that quantitatively compares the two populations.

Zoom In Zoom Out Reset image size

Using Bayes' theorem, with a uniform prior, we generated a binomial distribution for hot stars and cool stars separately that illustrates the probability of the data, given an assumption about what fraction of systems are single. Such an argument is similar to determining the fairness of a coin flip, where heads equates to a single system and tails a multi system. Specifically, we plot

$$\begin{eqnarray}&&{ \mathcal P }(\{{\rm{data}}\}| S)={{AS}}^{{N}_{{\rm{s}}}}{(1-S)}^{{N}_{{\rm{t}}}-{N}_{{\rm{s}}}},\end{eqnarray} \tag{ 30 }$$

where N_s is the number of single systems within a population of N_t total stars and A is a normalization coefficient (Sivia 1996). The variable S is the single bias weighting; the probabilistic tendency for a population of planet-hosting stars to display single transits as opposed to multiples. The quantity ${ \mathcal P }(\{{\rm{data}}\}| S)$ gives the probability of reproducing the data if the underlying tendency is S. In the hot population, N_t = 132 and N_s = 110, whereas the cool population had N_t = 1504 and N_s = 1098.

As can be seen from Figure 4, the two distributions are visually distinct, with hotter stars possessing a stronger bias toward singles, with a significance of 2.9σ.⁶ More data are needed to tease out this relationship further and to isolate the influence of a tilted star versus other confounding factors, such as the dependence upon stellar mass of the occurrence rate of giant planets. For now, we conclude that the data supports our general prediction, that hotter, more oblique stars possess a relatively greater abundance of single-transiting planets.

Zoom In Zoom Out Reset image size

A separate prediction relates to the distance of the planets from the host star. Specifically, the quadrupole moment falls off as ${R}_{\star }^{2}/{a}^{2}$, whereas the coplanarity required for transit grows as a/R_⋆, and so the overall magnitude of our proposed mechanism should become negligible within more distant systems. In consequence, we predict that the "Kepler Dichotomy" signal will weaken for systems at larger orbital distances. As future missions, such as TESS collect more data, this unique aspect of our model will become amenable to observational tests.

A caveat to the above analysis is that the dichotomy appears to persist within the population of planets around M-dwarfs (Ballard & Johnson 2016). This is problematic as these stars are expected to exhibit lower inclinations, being cooler on the main sequence. We interpret this as stellar oblateness being effective across all stellar masses, but being relatively more important within the hotter, more inclined population. This is supported by our numerical simulations, where even small obliquities reduced the transit number if oblateness was large enough (Figure 2).

5.3. The Origin of Spin–Orbit Misalignments

This paper investigates the origin of the "Kepler Dichotomy," within the context of primordially generated spin–orbit misalignments. We have shown that the quadrupole moment of such misaligned, young, fast-rotating stars is typically capable of exciting significant mutual inclinations between the hosted planetary orbits. In turn, the number of planets available for observation through transit around such a star is reduced, either through dynamical instability or directly as a result of the mutual inclinations, leaving behind an abundance of single-transiting systems (Johansen et al. 2012). The outcome is an apparent reduction in the multiplicity of tilted, hot stars, with their observed singles being slightly larger, as a consequence of many having undergone dynamical instabilities, in accordance with observations.

Through the conclusions of this work, the origins of hot Jupiters, of compact Kepler systems, the Kepler Dichotomy, and spin–orbit misalignments, are all placed within a common context.

This research is based in part upon work supported by NSF grant AST 1517936 and the NESSF Graduate Fellowship in Earth and Planetary Sciences (C.S). We would like to thank Erik Petigura, Henry Ngo, and Peter Gao for helpful discussions, along with Fred Adams and Joe O'Rourke for useful suggestions. We are grateful to the anonymous reviewer for helpful comments that significantly improved the manuscript.