Formation of Ultra-short-period Planets by Obliquity-driven Tidal Runaway

The spin vectors of close-in planets are subject to a dissipative tidal torque that moves them toward equilibrium configurations. The tidal torque arises due to the gravitational deformation (or "bulge") raised on the planet from its host star; it is dissipative because it involves this bulge sweeping across the planet every orbit. If the orbit is static, the tidally relaxed equilibrium of the spin vector is a straightforward spin-synchronous and aligned state, where the spin rotation frequency, ω = 2π/P_rot, is equal to the orbital mean motion, n = 2π/P, and the obliquity, , is zero. However, most planetary orbits are not static; they undergo precession due to interactions with other planets, the oblate host star, and any other gravitational sources that cause deviations from a 1/r potential. The equilibrium configurations of the spin pole in a uniformly precessing orbit frame are known as "Cassini states." The dynamical origin and behavior of the Cassini states have been documented in many previous works (e.g., Colombo 1966; Peale 1969, 1974; Ward 1975; Ward & Hamilton 2004; Correia 2015; Su & Lai 2020). Here, we summarize the most relevant material.⁷

Cassini states are configurations in which the precession rate of the planet's spin axis exactly matches that of its orbital plane. More specifically, the planetary spin axis, $\hat{{\boldsymbol{s}}}$, and unit orbit normal vector, $\hat{{\boldsymbol{n}}}$, precess at the same rate about the axis of the total system angular momentum vector, $\hat{{\boldsymbol{k}}}$. In a dissipationless Cassini state, these three vectors are coplanar. Dissipation causes $\hat{{\boldsymbol{s}}}$ to shift out of the plane defined by $\hat{{\boldsymbol{n}}}$ and $\hat{{\boldsymbol{k}}}$. For a given orbital inclination, I, with respect to the invariable plane, the obliquities of Cassini states obey the relation (e.g., Ward 1975)

$$\begin{eqnarray}&&g\sin (\epsilon -I)+\alpha \cos \epsilon \sin \epsilon =0.\end{eqnarray} \tag{ 1 }$$

Here, $g=\dot{{\rm{\Omega }}}$ is the precession frequency of the longitude of the ascending node. This frequency is negative (corresponding to nodal recession) for the cases of interest here. The frequency α is the spin-axis precession constant, which sets the precession period, T_α  =  2π/(α cos ), of the spin axis due to the torque induced by the host star on the oblate planet. In the absence of satellites orbiting the planet, α is given by (Neron de Surgy & Laskar 1997)

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{1}{2}\displaystyle \frac{{M}_{\star }}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{3}\displaystyle \frac{{k}_{2}}{C}\displaystyle \frac{\omega }{{\left(1-{e}^{2}\right)}^{3/2}}.\end{eqnarray} \tag{ 2 }$$

Here, M_⋆ is the stellar mass, M_p the planet mass, R_p the planet radius, a the semimajor axis, e the eccentricity, k₂ the planetary Love number, and C the planet's moment of inertia normalized by ${M}_{p}{{R}_{p}}^{2}$.

It is important to note that Cassini states are strictly only defined for uniform orbital precession, that is, when $g=\dot{{\rm{\Omega }}}$ and I are constant. When the precession is nonuniform, the planet's orbital inclination/node solution is composed of a superposition of several modes with multiple frequencies {g_i} and amplitudes {I_i}. The resulting spin vector equilibria are "quasi-Cassini states," which behave approximately like Cassini states with g in Equation (1) equal to one of the g_i modes. For example, Saturn's proposed Cassini state is associated with the g₈ inclination/node fundamental frequency of the solar system, which is dominated by Neptune's nodal precession (Hamilton & Ward 2004; Ward & Hamilton 2004). This concept of multiple modes will be revisited several times in this work. For now, we will simply assume that the frequency g corresponds to one of the g_i secular modes of the system.

With g and α frequencies specified, there are either two or four well-defined Cassini states, depending on the frequency ratio ∣g∣/α in reference to the critical frequency ratio,

$$\begin{eqnarray}&&{\left(| g| /\alpha \right)}_{\mathrm{crit}}={\left({\sin }^{2/3}I+{\cos }^{2/3}I\right)}^{-3/2}.\end{eqnarray} \tag{ 3 }$$

This critical ratio decreases as a function of I for I < 45°. When $| g| /\alpha \lt {\left(| g| /\alpha \right)}_{\mathrm{crit}}$, Equation (1) has four roots, corresponding to Cassini states 1–4. When ∣g∣/α > (∣g∣/α)_crit, Equation (1) has two roots, corresponding to Cassini states 2 and 3. States 1 and 2 are stable; state 3 is linearly stable but unstable to tidal evolution (Fabrycky et al. 2007); and state 4 is unstable. Thus, Cassini states 1 and 2 will be our primary focus, as they are the only ones that are stable in the presence of tides.

Cassini state 1 corresponds to the configuration in which $\hat{{\boldsymbol{s}}}$ and $\hat{{\boldsymbol{n}}}$ are on the same side of $\hat{{\boldsymbol{k}}}$. In this case, the convention is for the obliquity to be defined as negative. In state 2, $\hat{{\boldsymbol{s}}}$ and $\hat{{\boldsymbol{n}}}$ are on opposite sides of $\hat{{\boldsymbol{k}}}$, and the obliquity is positive. In the limit ∣g∣/α ≪ (∣g∣/α)_crit, the state 1 and 2 equilibrium obliquities are

$$\begin{eqnarray}&&{\epsilon }_{1}\approx {\tan }^{-1}\left(\displaystyle \frac{\sin I}{1+\alpha /g}\right);\ \ {\epsilon }_{2}\approx {\cos }^{-1}\left(\displaystyle \frac{-g\cos I}{\alpha }\right).\end{eqnarray} \tag{ 4 }$$

Figure 1 shows the obliquities of the Cassini states as a function of ∣g∣/α, through solving Equation (1). The top panel shows all four states for I = 5°, and the bottom panel shows the evolution of states 1 and 2 and (∣g∣/α)_crit with respect to I. It is important to emphasize that both states 1 and 2 have nonzero obliquities whenever I > 0°, making a misaligned planetary spin axis unavoidable. However, state 1's obliquity is very close to zero for small I and for ∣g∣/α ≲ 0.1. Meanwhile, state 2 allows for very large obliquities. The absolute value of the obliquity of both states increases with I, and in the limit ∣g∣/α ≫ (∣g∣/α)_crit, the state 2 obliquity asymptotes at ₂ = I. Depending on the initial conditions, tides will carry the spin vector to either state 1 or state 2. When ∣g∣/α > (∣g∣/α)_crit, state 2 is required. Thus, a high-obliquity state is often the only option.

Zoom In Zoom Out Reset image size

Assume, for now, that a planet initially tidally relaxes into Cassini state 2. (Later we will show that this is often the case.) Once this happens, the same dissipative tidal torque that brought the planet into the Cassini state will continue to generate heat in the planetary interior due to the nonzero obliquity. This dissipation affects the orbit as well; it generates orbital decay because the thermal energy dissipated in the interior is balanced by the conversion of orbital energy. The orbital decay, in turn, leads to a decrease in the ratio ∣g∣/α (assuming that the orbital precession is arising from planet–planet interactions) and an increase in the equilibrium obliquity of Cassini state 2. This is depicted in the upper panel of Figure 1.

While tidal dissipation and orbital decay necessarily result from nonzero obliquities, quantifying the magnitude of energy dissipation and orbital decay is nontrivial and there are many available tidal models (e.g., Efroimsky & Williams 2009; Ferraz-Mello 2013; Correia et al. 2014). The simplest and most widely used approach is equilibrium tide theory (e.g., Darwin 1880; Goldreich & Soter 1966; Mignard 1979; Hut 1981; Eggleton et al. 1998), which we adopt here using the viscous approach (Levrard et al. 2007; Leconte et al. 2010).

The basic assumptions are that the planet's gravitational response to the tidal forces constitutes a hydrostatic deformation, or tidal bulge, and this bulge lags the star's position with a constant time lag, Δt. The constant time offset is often parameterized in terms of the annual tidal quality factor, Q, which is related to Δt through Q = (Δtn)⁻¹. Q quantifies the efficiency of tidal damping, and it is combined with k₂ into the "reduced tidal quality factor," $Q^{\prime} =3Q/2{k}_{2}$. Q is highly uncertain for individual planets but is known in an order-of-magnitude sense for different planetary archetypes. Rocky bodies in the solar system have Q ∼ 100 (Goldreich & Soter 1966; Murray & Dermott 1999), while extrasolar super-Earths and sub-Neptunes have been found with Q ∼ 10³–10⁵ (Morley et al. 2017; Puranam & Batygin 2018), similar to the estimated values for Uranus (Tittemore & Wisdom 1989) and Neptune (Zhang & Hamilton 2008). We assume a range of plausible planetary Q values in this work. As for k₂, constraints come from both the Solar System bodies (Lainey 2016) and from theoretical models (Kramm et al. 2011; Kellermann et al. 2018), which we will use to inform our fiducial values in this work.

Within the viscous equilibrium tide model, the tidal luminosity—or the rate at which orbital energy is converted into thermal energy—is given by the expression (Levrard et al. 2007; Leconte et al. 2010):

$$\begin{eqnarray}&&{L}_{\mathrm{tide}}(e,\epsilon )=2K\left[{N}_{a}(e)-\displaystyle \frac{{N}^{2}(e)}{{\rm{\Omega }}(e)}\displaystyle \frac{2{\cos }^{2}\epsilon }{1+{\cos }^{2}\epsilon }\right].\end{eqnarray} \tag{ 5 }$$

Here, N_a(e), N(e), Ω(e) are functions of eccentricity given by

$$\begin{eqnarray}&&{N}_{a}(e)=\displaystyle \frac{1+\tfrac{31}{2}{e}^{2}+\tfrac{255}{8}{e}^{4}+\tfrac{185}{16}{e}^{6}+\tfrac{25}{64}{e}^{8}}{{\left(1-{e}^{2}\right)}^{\tfrac{15}{2}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&N(e)=\displaystyle \frac{1+\tfrac{15}{2}{e}^{2}+\tfrac{45}{8}{e}^{4}+\tfrac{5}{16}{e}^{6}}{{\left(1-{e}^{2}\right)}^{6}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{\Omega }}(e)=\displaystyle \frac{1+3{e}^{2}+\tfrac{3}{8}{e}^{4}}{{\left(1-{e}^{2}\right)}^{\tfrac{9}{2}}}.\end{eqnarray} \tag{ 8 }$$

The magnitude of L_tide is dictated by K, the characteristic luminosity scale,

$$\begin{eqnarray}&&K=\displaystyle \frac{3n}{2}\displaystyle \frac{{k}_{2}}{Q}\left(\displaystyle \frac{G{{M}_{\star }}^{2}}{{R}_{p}}\right){\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{6}.\end{eqnarray} \tag{ 9 }$$

Equation (5) assumes that the planet's spin rotation frequency has reached its equilibrium rate, given by

$$\begin{eqnarray}&&{\omega }_{\mathrm{eq}}=n\displaystyle \frac{N(e)}{{\rm{\Omega }}(e)}\displaystyle \frac{2\cos \epsilon }{1+{\cos }^{2}\epsilon }.\end{eqnarray} \tag{ 10 }$$

The tidal luminosity is balanced by a decrease in the orbital energy, such that ${L}_{\mathrm{tide}}\,=\,-({{GM}}_{\star }{M}_{p}\dot{a})/(2{a}^{2})$. Using this, we may calculate the orbital decay timescale:

$$\begin{eqnarray}&&{\tau }_{a}=\displaystyle \frac{a}{\dot{a}}=-\displaystyle \frac{{{GM}}_{\star }{M}_{p}}{4{aK}}{\left[{N}_{a}(e)-\displaystyle \frac{{N}^{2}(e)}{{\rm{\Omega }}(e)}\displaystyle \frac{2{\cos }^{2}\epsilon }{1+{\cos }^{2}\epsilon }\right]}^{-1}.\end{eqnarray} \tag{ 11 }$$

The τ_a timescale can be used to delineate the regime in which planets undergo significant tidal migration for a given initial semimajor axis. Note that all planets in the system, not just the innermost one, can be migrating due to nonzero obliquities and/or eccentricities. However, it is generally only the innermost planet that can migrate fast enough to substantially separate itself from its neighbors within the age τ_age of the system, because the τ_a timescale depends strongly on a. Figure 2 shows ∣τ_a∣ as a function of P and for a fiducial set of rocky planet parameters. When ∣τ_a∣ ≲ τ_age ≈ 1–10 Gyr, corresponding to P ≲ 2–3 days for  ≈ 10°, the semimajor axis will decrease by order unity during the system lifetime. The instantaneous ∣τ_a∣ timescale is an overestimate of the total time to decay, however, because orbital migration is a runaway process in which ∣τ_a∣ decreases as the orbit shrinks and the obliquity (in Cassini state 2) increases. Accordingly, it is useful to calculate the time, τ_decay, for complete decay from some initial a = a_i to an ending position with a ≈ 0. Using the fact that $\dot{a}\propto {a}^{-11/2}$ for constant obliquity and eccentricity, one can show that

$$\begin{eqnarray}&&{\tau }_{\mathrm{decay}}\approx \displaystyle \frac{2}{13}| {\tau }_{a}(a={a}_{i})| .\end{eqnarray} \tag{ 12 }$$

Thus, for some initial a_i, the timescale for complete decay of order Δa ∼ a_i is nearly an order of magnitude smaller than the decay timescale at the initial separation, ∣τ_a(a = a_i)∣. In practice, a more accurate estimate of τ_decay than that provided by Equation (12) can be obtained through numerical integration of $\dot{a}$ in Equation (11) using an evolving e and . This will be our approach in Section 3.3.

Zoom In Zoom Out Reset image size

An orbit will only decay completely, however, if it maintains a nonzero eccentricity and/or obliquity. If these go to zero, the tidal dissipation and migration will stall. This stalling is expected due to the tidal breaking of Cassini state 2 at high obliquity.

Cassini state 2 cannot exist at an arbitrarily high obliquity in the presence of tides (Fabrycky et al. 2007; Peale 2008). As the obliquity increases due to the orbital decay, there becomes a point at which the dissipative tidal torque overwhelms the orbital precession torque; Cassini state 2 breaks and the obliquity damps down to Cassini state 1. Fabrycky et al. (2007) showed that the breaking obliquity is related to the limits of a specific phase shift that appears in Cassini state 2 in the presence of steady tidal dissipation. The spin axis, $\hat{{\boldsymbol{s}}}$, shifts out of the plane containing $\hat{{\boldsymbol{n}}}$ and $\hat{{\boldsymbol{k}}}$, and this phase shift provides a balance of the tidal torque, up until the breaking point.

In order to identify this breaking obliquity, we consider the secular equation of motion of the spin vector, ${\boldsymbol{\omega }}=\omega \hat{{\boldsymbol{s}}}$, using the framework of Eggleton & Kiseleva-Eggleton (2001) and Fabrycky et al. (2007) with zero orbital eccentricity. The equation may be written as the sum of two torques: a nondissipative torque due to the star's nonuniform gravitational force on the planet, and a dissipative tidal torque due to the lagged response of the tidal bulge raised in the planet. The nondissipative torque generates the spin axis precession, and the dissipative tidal torque drives the spin vector toward the Cassini states. Explicitly, we may write

$$\begin{eqnarray}&&\,\dot{{\boldsymbol{\omega }}}={\dot{{\boldsymbol{\omega }}}}_{\mathrm{prec}}+{\dot{{\boldsymbol{\omega }}}}_{\mathrm{tides}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{\dot{{\boldsymbol{\omega }}}}_{\mathrm{prec}}=\alpha \omega (\hat{{\boldsymbol{s}}}\cdot \hat{{\boldsymbol{n}}})(\hat{{\boldsymbol{s}}}\times \hat{{\boldsymbol{n}}})\\ {\dot{{\boldsymbol{\omega }}}}_{\mathrm{tides}}=\displaystyle \frac{n}{C}{\left(\displaystyle \frac{a}{{R}_{p}}\right)}^{2}\left[-\displaystyle \frac{{\boldsymbol{\omega }}}{2{{nt}}_{F}}+\left(1-\displaystyle \frac{{\boldsymbol{\omega }}\cdot \hat{{\boldsymbol{n}}}}{2n}\right)\displaystyle \frac{\hat{{\boldsymbol{n}}}}{{t}_{F}}\right].\end{array}\end{eqnarray} \tag{ 14 }$$

Here, t_F is the tidal friction timescale given by

$$\begin{eqnarray}&&{t}_{F}=\displaystyle \frac{4Q^{\prime} }{9}{\left(\displaystyle \frac{a}{{R}_{p}}\right)}^{5}\displaystyle \frac{{M}_{p}}{{M}_{\star }}\displaystyle \frac{1}{n}.\end{eqnarray} \tag{ 15 }$$

Inspecting the expressions for ${\dot{{\boldsymbol{\omega }}}}_{\mathrm{prec}}$ and ${\dot{{\boldsymbol{\omega }}}}_{\mathrm{tides}}$, we observe that both torques exhibit the same dependencies with most physical parameters of the problem, including M_⋆, M_p, R_p, a, k₂, and C. The two exceptions are Q, which enters into ${\dot{{\boldsymbol{\omega }}}}_{\mathrm{tides}}$ but not ${\dot{{\boldsymbol{\omega }}}}_{\mathrm{prec}}$, and I, which factors into the equations via $\hat{{\boldsymbol{n}}}$. These dependencies imply that the breaking obliquity depends only on Q and I. Accordingly, in order to identify the breaking obliquity, we can simply select arbitrary system parameters and numerically integrate Equation (13) in response to an evolving ∣g∣/α. Doing this for different values of Q and I and determining the maximum obliquity in each case will provide the full range of outcomes.

For the purposes of this calculation, we will assume the normal vector $\hat{{\boldsymbol{n}}}$ precesses uniformly about the total angular momentum vector $\hat{{\boldsymbol{k}}}$ with a constant inclination I between them. We will initialize the system with ∣g∣/α = 3 > (∣g∣/α)_crit and let ∣g∣/α exponentially decrease with a timescale equal to ten times the adiabatic limit given by Su & Lai (2020), i.e., firmly in the adiabatic regime. The planet starts in Cassini state 2 with an obliquity close to  ∼ I (see Figure 1). As ∣g∣/α decreases, the obliquity rises up the Cassini state 2 branch until the dissipative torque becomes too strong and the planet can no longer be maintained in the high-obliquity state. Figure 3 shows the numerically calculated breaking obliquity for a range of values of Q and I. We observe that the limit increases with both Q and I, indicating that planets with such properties can undergo more orbital decay before tidal breaking.

Zoom In Zoom Out Reset image size

After Cassini state 2 is destabilized, the obliquity damps down from its excited state and settles into Cassini state 1. From Equation (14), we see that this equilibration occurs on a fast timescale of roughly

$$\begin{eqnarray}&&\begin{array}{l}{\tau }_{\mathrm{equil}}\approx {t}_{f}C{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{2}\\ =\,135\ \mathrm{yr}\left(\displaystyle \frac{Q^{\prime} }{{10}^{3}}\right){\left(\displaystyle \frac{a}{0.03\mathrm{au}}\right)}^{\tfrac{9}{2}}\left(\displaystyle \frac{{\rho }_{p}}{{\rho }_{\oplus }}\right)\left(\displaystyle \frac{C}{0.35}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-\tfrac{3}{2}}.\end{array}\end{eqnarray} \tag{ 16 }$$

Once in Cassini state 1, the obliquity is small but nonzero (Equation (4) and Figure 1). For instance, when I = 5° and ∣g∣/α = 0.1, the obliquity of Cassini state 1 is equal to ₁ ≈ −0.5°. The planet may experience further orbital decay while in Cassini state 1, but it will generally be slow and stable because the Cassini state 1 obliquity decreases in magnitude as ∣g∣/α decreases.

We consider the innermost planet in a close-in, multiplanet system that is not perfectly coplanar. Mutual inclination constraints will be discussed in Section 5, but they do not matter in detail for now. In addition, we will adopt a three-planet system. This again does not strongly affect the overall picture. Although working with two-planet systems simplifies the dynamics, it is less generalizable, and three-planet systems are more representative of the observed multiplanet systems with USPs.

In the process of mapping out the parameter space, there are many system properties to consider, including M_p, R_p, a, period ratios P_j+1/P_j between neighboring planets, etc. We will reduce the exploration down to three essential parameters: M_p, a_1,i (the innermost planet's initial semimajor axis), and P_j+1/P_j. We assume that all planets in the system have uniform masses and orbital spacing, a simplifying assumption that we will later relax (Section 3.3.1) but which is justified on the basis of the observed intrasystem uniformity of Kepler multis (Weiss et al. 2018; Millholland et al. 2017). In addition, we assume that the inner planet has an Earth-like composition (approximately 30% Fe, 70% MgSiO₃), and we use this assumption to calculate R_p1 for a given M_p1 (where the subscript "1" refers to the innermost planet) using the mass–radius tables from Zeng et al. 2016. (In Section 6.2, we will discuss the implications of possible early mass loss from the inner planet, which would be expected if the planet formed with a primordial envelope of hydrogen/helium.)

Apart from the three essential parameters (M_p, a_1,i, and P_j+1/P_j), there are several additional parameters that we will hold fixed. We will take the inner planet's Love number and moment of inertia factor (which enter into α in Equation (2)) to represent fiducial values for terrestrial planets, as determined observationally for Solar System bodies (Murray & Dermott 1999; Lainey 2016) and theoretically for extrasolar super-Earths (Kramm et al. 2011; Kellermann et al. 2018). We will use Q = 10³, k₂ = 0.4, and C = 0.35, noting that changes in these parameters will only affect our results on a detailed level. We will also assume that the planet's spin rate is at equilibrium, ω = ω_eq, or approximately synchronous with the mean motion, ω = n, when eccentricities and obliquities are small.⁸ As for mutual inclinations, we will take the innermost planet to be misaligned by 10° with respect to the outer planets (approximately consistent with Dai et al. (2018)), with a 1° mutual inclination between the outer planets. Finally, we will consider two different stellar masses, M_⋆ = 0.6M_⊙ (K-dwarf star) and M_⋆ = 1.0M_⊙ (G-dwarf star).

Before proceeding, we note that evolution by way of obliquity tides must conserve angular momentum in spite of orbital decay. Accordingly, the inclination of the inner planet must change as its orbit shrinks (Fabrycky et al. 2007). Such inclination evolution does not affect this section substantially, since Cassini state 2 varies little with inclination at the high-obliquity end (Figure 1). However, the inclination evolution is important to incorporate into our theory because angular momentum constraints can limit the total extent of the migration. Thus, we redress this omission in Section 5 and Appendix, where we develop a secular model that self-consistently accounts for the tidal semimajor axis and inclination evolution.

The primary factor determining whether the innermost planet will become a USP is its Cassini state, and this depends most strongly on ∣g∣/α. Accordingly, our goal is to calculate this frequency ratio for the innermost planet across the parameter space we have just identified. While the spin-axis precession constant α has a straightforward analytic expression (Equation (2)), g is more complex. The planet's orbit nodal precession arises from gravitational interactions with its (potentially oblate) host star and neighboring planets. Assuming the orbits are nonresonant, the set of orbital eigenfrequencies, {g_i}, may be approximated using Laplace–Lagrange secular theory.⁹ To second order in the eccentricities and inclinations, this solution depends only on the masses and semimajor axes, and the eccentricity and inclination solutions are decoupled. As discussed in Section 2.1, the relevant nodal frequency g for the Cassini state will be one of the multiple eigenfrequencies identified in this solution.

We begin by constructing a planetary disturbing function for N planets orbiting an oblate host star (Murray & Dermott 1999). The disturbing function is the non-Keplerian perturbing gravitational potential experienced by the planets due to their mutual interactions. Keeping only terms associated with the inclinations to second order, the disturbing function takes the form

$$\begin{eqnarray}&&\left\langle {{ \mathcal R }}_{j}^{(\sec )}\right\rangle ={n}_{j}{a}_{j}^{2}\left[\displaystyle \frac{1}{2}{B}_{{jj}}{I}_{j}^{2}+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k=1}{k\ne j}}^{N}{B}_{{jk}}{I}_{j}{I}_{k}\cos ({{\rm{\Omega }}}_{j}-{{\rm{\Omega }}}_{k})\right],\end{eqnarray} \tag{ 17 }$$

where the subscript j is the planet number, n is the mean motion, I is the inclination, and Ω is the longitude of the ascending node. The quantities B_jj and B_jk represent the interaction coefficients within the matrix B and are given by

$$\begin{eqnarray}&&\begin{array}{l}{B}_{{jj}}=-{n}_{j}\left[\displaystyle \frac{3}{2}{J}_{2\star }{\left(\displaystyle \frac{{R}_{\star }}{{a}_{j}}\right)}^{2}-\displaystyle \frac{27}{8}{J}_{2\star }^{2}{\left(\displaystyle \frac{{R}_{\star }}{{a}_{j}}\right)}^{4}\right.\\ \,\left.+\,\displaystyle \frac{1}{4}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k=1}{k\ne j}}^{N}\displaystyle \frac{{M}_{{pk}}}{{M}_{\star }+{M}_{{pj}}}{\alpha }_{{jk}}{\bar{\alpha }}_{{jk}}{b}_{3/2}^{(1)}({\alpha }_{{jk}})\right]\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{B}_{{jk}}=\displaystyle \frac{1}{4}\displaystyle \frac{{M}_{{pk}}}{{M}_{\star }+{M}_{{pj}}}{n}_{j}{\alpha }_{{jk}}{\bar{\alpha }}_{{jk}}{b}_{3/2}^{(1)}({\alpha }_{{jk}})(j\ne k).\end{eqnarray} \tag{ 19 }$$

Here, M_pj is the mass of the jth planet. When a_j < a_k, ${\alpha }_{{jk}}={\bar{\alpha }}_{{jk}}={a}_{j}/{a}_{k}$, and when a_j > a_k, α_jk = a_k/a_j and ${\bar{\alpha }}_{{jk}}=1$. The quantity ${b}_{3/2}^{(1)}$ is a Laplace coefficient, defined by

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

In Equation (18), J_2⋆ is the star's second gravitational (quadrupole) moment. J_2⋆ can be expressed in terms of R_⋆, the stellar spin rate, ω_⋆ = 2π/P_⋆, and the tidal Love number k_2⋆ using the approximate relationship (Sterne 1939; Ward et al. 1976; Spalding & Batygin 2016)

$$\begin{eqnarray}&&\begin{array}{l}{J}_{2\star }\approx \displaystyle \frac{1}{3}{k}_{2\star }\displaystyle \frac{{\omega }_{\star }^{2}}{{{GM}}_{\star }/{R}_{\star }^{3}}\\ \sim {10}^{-3}\left(\displaystyle \frac{{k}_{2\star }}{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{ 21 }$$

where $\sqrt{{{GM}}_{\star }/{R}_{\star }^{3}}$ is the breakup rotational velocity. Here, we have used fiducial values appropriate to young, rapidly rotating stars (Batygin & Adams 2013; Spalding & Batygin 2016). As we will show, the inclusion of J_2⋆ is not required for the mechanism but does make it more efficient. Going forward, unless otherwise noted, we will use J_2⋆ = 10⁻⁴ to represent a typical star within the first several 100 Myrs of evolution.

Given the form of the disturbing function in Equation (17), it is convenient and customary to use a transformation to the inclination "vectors," defined by

$$\begin{eqnarray}\begin{array}{rcl}{p}_{j} & = & {I}_{j}\sin {{\rm{\Omega }}}_{j}\\ {q}_{j} & = & {I}_{j}\cos {{\rm{\Omega }}}_{j}.\end{array}\end{eqnarray} \tag{ 22 }$$

The solutions to the equations of motion are then

$$\begin{eqnarray}\begin{array}{rcl}{p}_{j}(t) & = & \displaystyle \sum _{i=1}^{N}{I}_{{ji}}\sin ({g}_{i}t+{\gamma }_{i})\\ {q}_{j}(t) & = & \displaystyle \sum _{i=1}^{N}{I}_{{ji}}\cos ({g}_{i}t+{\gamma }_{i}),\end{array}\end{eqnarray} \tag{ 23 }$$

where the {g_i} are the N eigenvalues of the matrix B and {I_ji} are the corresponding eigenvectors.¹⁰ Since the eigenvectors of B are only defined up to a scaling factor, one may use the initial conditions to determine the magnitudes of the eigenvectors and the phases γ_i. Finally, the time evolutions of the inclination and node are given by

$$\begin{eqnarray}\begin{array}{rcl}{I}_{j}(t) & = & {\left[{\left[{p}_{j}(t)\right]}^{2}+{\left[{q}_{j}(t)\right]}^{2}\right]}^{1/2}\\ {{\rm{\Omega }}}_{j}(t) & = & {\tan }^{-1}\left[\displaystyle \frac{{p}_{j}(t)}{{q}_{j}(t)}\right].\end{array}\end{eqnarray} \tag{ 24 }$$

Equations (23) and (24) indicate that the inclination/node solution is composed of a superposition of modes from the N secular eigenfrequencies, {g_i}. Any of these modes may be the g frequency that is dominant for a planet's Cassini state, such as in the case of Saturn, where the relevant g is that which is dominated by Neptune's nodal precession (Ward & Hamilton 2004; Hamilton & Ward 2004). Determining the important {g_i} mode is challenging, but for our regime of interest it will generally be the one that is closest to the spin-axis precession constant α. In the sections that follow, we will take the fastest frequency $| g{| }_{\max }$ as the dominant mode. We will show in Section 3.3.2 that this is a good approximation. A robust determination of which of the {g_i} modes dominates is the biggest area for follow-up of this work. This would include an investigation of when and how chaos arises from overlapping modes. We will return to this idea in the Discussion (Section 6).

With the M_p—a_1,i—P_j+1/P_j parameter space and the calculation of ∣g∣/α now specified, we can identify the region of this space that is susceptible to producing a USP via obliquity-driven tidal migration. To begin, we simply plot in Figure 4 ∣g∣/α as a function of M_p and a_1,i for a fixed P_j+1/P_j. There are several aspects to note. First, for a fixed M_p, the ratio ∣g∣/α increases with a_1,i, while both ∣g∣ and α independently decrease with increasing a_1,i. This highlights the fact that ∣g∣/α will shrink upon the inner planet's orbital decay, as depicted in Figure 1. (During the decay, however, the period ratio is also increasing, which leads to a more rapid decrease of ∣g∣/α than for a fixed P_j+1/P_j.)

Zoom In Zoom Out Reset image size

A second observation is that the ∣g∣/α = 1 contour crosses through the middle of this parameter space. For ∣g∣/α > (∣g∣/α)_crit ≈ 1, the inner planet is guaranteed to occupy Cassini state 2 (see Figure 1). This is not to say that the planet cannot initially be captured into Cassini state 2 if ∣g∣/α < (∣g∣/α)_crit. For ∣g∣/α < (∣g∣/α)_crit and when the obliquity is greater than the dashed separatrix curves in Figure 1, the obliquity tidally relaxes into state 2. In addition, the obliquity can also be resonantly excited into state 2 (e.g., Millholland & Laughlin 2019). However, in this work we focus on the ∣g∣/α > (∣g∣/α)_crit ≈ 1 regime, where capture into Cassini state 2 is inevitable.

Given that the ∣g∣/α = 1 contour delineates the region of guaranteed participation in Cassini state 2, we can isolate this contour and plot its variation with respect to the third parameter, P_j+1/P_j, which was held fixed in Figure 4. The result is shown in Figure 5, where the solid lines indicate the ∣g∣/α = 1 contours for a range of P_j+1/P_j given by the color bar. For increasing P_j+1/P_j, ∣g∣ decreases, so the ∣g∣/α = 1 contour moves to larger a_1,i.

Zoom In Zoom Out Reset image size

In addition to the ∣g∣/α ≳ 1 constraint, USP production also requires that a_1,i is small enough that τ_decay is sufficiently fast (i.e., ≲1 Gyr). To calculate τ_decay, we can integrate $\dot{a}$ from Equation (11) while asserting that the obliquity is in a Cassini state solved numerically using Equation (1). The obliquity begins in Cassini state 2,  = ₂. As discussed in Section 2.2, the orbital migration leads ₂ to increase, and the decay rate eventually reaches runaway. Accordingly, τ_decay is the time it takes until the tidal breaking of Cassini state 2 (recall Section 2.3). The obliquity subsequently settles into Cassini state 1,  = ₁, over a short timescale (∼100 yr, Equation (16)) that we approximate as instantaneous. Finally, in this low-obliquity state, the rapid orbital decay stalls.

For each point in the M_p—a_1,i—P_j+1/P_j grid, we integrate $\dot{a}$ from Equation (11) over 10 Gyr and calculate τ_decay as just described. Examples for three initial grid points are shown in Figure 6, where we plot the time evolution of P₁ and . All three examples undergo a period of runaway decay, although they reach different final orbital periods at different τ_decay times. For similar initial period ratios, these outcomes depend most strongly on the initial semimajor axes, a_1,i.

Zoom In Zoom Out Reset image size

As a result of performing these orbital decay evolutions for all grid points, the dashed lines in Figure 5 indicate the τ_decay = 1 Gyr contours for a range of P_j+1/P_j values. For a given contour, the area to the left corresponds to shorter τ_decay. Note that, for larger P_j+1/P_j at fixed a_1,i and M_p, the orbital precession frequency ∣g∣ is smaller, such that the initial ₂ is larger and τ_decay is smaller. This accounts for the observation that the τ_decay = 1 Gyr contours shift to the right for larger P_j+1/P_j. Similarly, comparing the top and bottom panels of Figure 5, we see that all contours are further to the right (larger a_1,i) for the M_⋆ = 1.0M_⊙ case compared to the M_⋆ = 0.6M_⊙ case. This is related to the fact that, with all other parameters held fixed, a larger M_⋆ yields a smaller ∣g∣/α and τ_decay.

Taken together, the regions with ∣g∣/α > 1 and τ_decay < 1 Gyr delineate the parts of parameter space (shaded in Figure 5) where the innermost planets are most susceptible to becoming a USP. (We note that these regions have been plotted for the illustrative set of period ratios and thus do not indicate the full parameter space.) It is intriguing to note that there is a larger susceptible parameter space (meaning more frequent USP production) for the case with M_⋆ = 0.6M_⊙ compared to that with M_⋆ = 1.0M_⊙. This trend is in the same direction as the observational occurrence rates derived by Sanchis-Ojeda et al. (2014), who found that USPs are more common around smaller-mass stars. While this comparison is suggestive, we would also need to know the planet occurrence rate as a function of M_p, a, and M_⋆, to robustly determine that the theory has made a correct prediction.

3.3.1. Unequal Planet Masses

In our analysis thus far, we have adopted the simplifying assumption of equal-mass planets in a compact, equally spaced system. This configuration is not always a good approximation, however, and it is useful to consider systems hosting more massive planets on wider exterior orbits. Still adopting a three-planet system, we keep the innermost planet's mass, M_p1, fixed and examine a range of masses and separations for the two exterior planets. We parameterize this using M_p,ext/M_p1 (where M_p,ext is the mass of each exterior planet) and P_j+1/P_j.

Figure 7 shows contours of ∣g∣/α = 1 in the M_p,ext/M_p1—P_j+1/P_j space for M_p1 = 6 M_⊕ and for different values of a_1,i. (Note that the ∣g∣/α ratio depends much more strongly on a_1,i than M_p1, so varying M_p1 does not change this picture much.) The contours illustrate that more massive exterior planets with wider separations can have the same effect on the orbital precession frequency as smaller and closer perturbers. Even for M_p,ext/M_p1 ∼ 100, however, ∣g∣/α > 1 requires that P_j+1/P_j < 10.

Zoom In Zoom Out Reset image size

3.3.2. Initial Obliquity Evolution

We have assumed throughout Section 3.3 that the fastest frequency $| g{| }_{\max }$ is the dominant of the {g_i} modes, but we have not yet justified this claim. Here, we conduct a numerical study of the innermost planet's initial obliquity evolution in order to examine which mode is dominant early on and which Cassini state the planet initially settles into. We take the secular orbital solution (e.g., Equation (24)) to represent the inclination and node evolution, and evolve the spin vector using the secular equations of motion (Equations (13) and (14)). We perform these evolutions for the same M_p—a_1,i—P_j+1/P_j grid of initial conditions presented earlier. We consider a 0° initial planet obliquity and a 0.5 day primordial planetary rotation period. We also assume Q₁ = 300.

Figure 8 shows three examples of the innermost planet's initial spin vector evolution. Here, we use M_p = 6 M_⊕ and P_j+1/P_j = 1.4, but we show the results for different semimajor axes for the innermost planet. In the three examples, all of which have $| g{| }_{\max }/\alpha \gt {(| g| /\alpha )}_{\mathrm{crit}}$, the obliquity starts from 0°, undergoes a transient period of chaotic excitation and settles into libration around Cassini state 2 with $| g| =| g{| }_{\max }$ (the fastest secular eigenfrequency).

Zoom In Zoom Out Reset image size

After performing these integrations for the full M_p—a_1,i—P_j+1/P_j grid, we find good agreement with the analytic parameter boundaries identified earlier in Section 3.3. That is, whenever $| g{| }_{\max }/\alpha \gtrsim 1$ (as represented with the solid contours in Figure 5), the obliquity is always captured into Cassini state 2. About half the time, the Cassini state 2 is with the fastest frequency ($| g| =| g{| }_{\max }$), and the other times, it is with the second-fastest frequency. Cassini state 2 with the second-fastest frequency has a higher obliquity than that with $| g{| }_{\max }$. Accordingly, even if the spin vector temporarily settles into libration around Cassini state 2 with the second-fastest frequency and later breaks out of it, it may still encounter Cassini state 2 with $| g{| }_{\max }$ as the obliquity damps back down. Finally, we also find that whenever $| g{| }_{\max }/\alpha \lesssim 1$, the obliquity is always captured into Cassini state 1 with $| g| =| g{| }_{\max }$.

The results of these integrations thus confirm that the analytical simplification of using $| g{| }_{\max }$ as the dominant mode is appropriate. Cassini state 2 with $| g| =| g{| }_{\max }$ is a frequent occurrence whenever $| g{| }_{\max }/\alpha \gt 1$. However, these secular modes may not be readily distinguishable if the system parameters lead to resonance overlap. In this case, the planetary spin vector would be susceptible to chaotic evolution (see Section 6.1). Further work is therefore required to better understand which initial state is most likely for arbitrary starting parameters.

To better understand the limitations imposed by angular momentum conservation, it is instructive to examine the inner planet's inward migration, while adhering to the above framework. Angular momentum is conserved by way of tidal torques upon the planet being transferred to the orbit. In order to account for these torques within a secular model, we present in Appendix an extension of the secular model in Section 3.2 that incorporates angular momentum conservation into the inner planet's orbital decay (see also Chyba et al. 1989). This is accomplished by way of a forced decay of the orbital inclination over a timescale τ_I while the semimajor axis decays.

In this section, we solve the inclination and semimajor axis evolution Equations (A1) and (A4) subject to initial conditions (A12) and a_1,i. We choose an inclination-damping timescale τ_I such that the inner planet becomes a USP (i.e., reaches P₁ ≲ 1 day) within ∼100 Myr. As discussed above, the true time taken to form a USP via obliquity tides can be longer. However, for the purposes of this simulation, the choice of τ_I is arbitrary provided that it exceeds the timescale of the secular oscillations.

As illustrated in Figure 2, the tidal migration rate increases as the semimajor axis decays, leading to a runaway migration. Accordingly, we modulate the timescale of inclination damping by a factor of ${({a}_{1}/{a}_{1,i})}^{5}$, prescribing ${\tau }_{I}=0.1\,\mathrm{Myr}\,{({{\rm{a}}}_{1}/{{\rm{a}}}_{1,i})}^{5}$. We carry out two secular integrations, each lasting 100 Myr for a star with M_⋆ = M_⊙, R_⋆ = R_⊙, and J_2⋆ = 10⁻⁴. Each planet is given the same mass of M_p1 = M_p2 = 5 M_⊕, and they begin at semimajor axes a_1,i = 0.03 au and a_2,i = 0.05 au. From angular momentum conservation (see Equation (29)), the minimum required tilt between the orbital and stellar angular momenta is I_p⋆ ≳ 24°. Therefore, we illustrate the importance of the angular momentum constraint by choosing one case to begin above the critical misalignment, at I_p⋆ = 30°, and the second case to possess insufficient inclination, with I_p⋆ = 15°. These tilts predict, respectively, innermost achievable periods of P_1,f = 0.63 days and P_1,f = 1.5 days.

The results of our secular integrations are displayed in Figure 11. In the top panel, we show the evolution of the inner planet's orbital period during tidal evolution. The horizontal line indicates the minimum period achievable due to angular momentum constraints (Equation (28)). The bottom panels track the inner (green) and outer (blue) planetary orbital inclinations, with their mutual inclinations shown in gray. Left panels correspond to the greater initial stellar obliquity of I_p⋆ = 30° and the right panels show I_p⋆ = 15°.

Zoom In Zoom Out Reset image size

As expected from the discussion above, only when the initial stellar obliquity is sufficiently large is the inner planet able to migrate far enough inward to become a USP, with P₁ < 1 day (left). From the bottom left panel, we see that the evolution is characterized by transient, oscillatory behavior for t ≲ 0.3 Myr, or a few inclination decay timescales. During this transient period, the system's evolution is governed by a mix of two eigenmodes. At this stage, the inner planet's inclination decays faster than the outer planet's and the mutual planet–planet inclination grows.

Physically, the early growth in mutual inclinations arises from the inner planet becoming more dominated by the stellar quadrupole during its inward migration (Becker et al. 2020; Li et al. 2020). Its orbit reorients to the stellar equator, while the outer planet remains inclined. Thus, if the inner planet becomes a USP at this stage, the two planets will exhibit a large mutual inclination, as observed (Dai et al. 2018). A USP formed in this way is predicted to exhibit a low misalignment with respect to the stellar spin axis, lower than its exterior companion. We return to this prediction in Section 6.3.

After the initial transient evolution, the system collapses onto a single eigenmode (Zhang et al. 2013; Pu & Lai 2019), as the outer planet's inclination begins to undergo substantial decay. The inner planet's period falls well below P₁ = 1 day while the mutual inclination between the planets remains high. In this simple example, tidal migration continues until P₁ = 0.6 days, aligning the orbits. However, as discussed above, the planet typically breaks out of the Cassini state before reaching the ultimate tidal end state (see Figures 3 and 6).

When insufficient orbital inclination exists (right panel of Figure 11), the orbits tidally realign before the inner planet migrates below P₁ = 1 day. Once this state is reached, the planetary obliquity falls to zero, and tidal migration stalls. Accordingly, angular momentum places a strong constraint upon USP formation via obliquity tides. However, here we considered constraints due to one exterior planet at 0.05 au. Angular momentum constraints become significantly less restrictive as the outer planet's angular momentum increases (Figure 10) or if the number of exterior planets increases. Moreover, if some fraction of the orbital decay occurs through eccentricity tides (Petrovich et al. 2019; Pu & Lai 2019), this will also weaken the constraint, since in that case orbital circularization would contribute to angular momentum conservation.

The solution for a multiple-planet system's inclination/node evolution generally involves a set of frequencies {g_i} and amplitudes {I_ji}, such that each planet's orbital evolution can be approximated by a superposition of these modes. As discussed several times throughout this work, a planet's Cassini state may be established with any one of these frequency modes, typically that which is closest to the planet's spin-axis precession constant α. However, the process of identifying the dominant frequency mode is complicated (e.g., Peale 1974), and we have not addressed it in this work.

Moreover, when α is close to several orbital eigenfrequencies, or combinations of these frequencies, obliquity chaos may result from secular spin–orbit resonance overlap (Saillenfest et al. 2019). All of the terrestrial planets in the solar system have wide ranges of possible spin states in which their obliquities undergo large-amplitude chaotic variations (Laskar & Robutel 1993). Although Mars is the only planet that is currently undergoing strong (∼60° amplitude) chaotic variations (Ward 1973; Touma & Wisdom 1993), the obliquities of the other terrestrial planets were likely chaotic in the past as well. The obliquities of Mercury and Venus have been stabilized by tides (e.g., Peale 1974; Correia et al. 2003), whereas Earth's obliquity would undergo chaotic variations if not for the Moon's stabilizing influence (Neron de Surgy & Laskar 1997; Li & Batygin 2014).

More work must be done to understand the prevalence of chaos in obliquity dynamics of short-period exoplanets, as well as the stability of Cassini states when multiple close frequency modes are present. Such dynamics could affect our proposed USP production mechanism. For instance, chaotic dynamics could prematurely knock planets out of their high-obliquity Cassini states before the planets have fully migrated. Chaos is not necessarily destructive for the mechanism, though, since high-amplitude obliquity variations—such as those experienced by Mars—would still lead to large-scale, obliquity-driven orbital decay. Future studies of both chaos and tides within short-period planets will help inform these potential scenarios.

If USP production happens early (≲1 Gyr), as we have postulated, multiple system parameters may be changing simultaneously during the orbital migration process. The resulting evolution is more complex than we have depicted. For instance, the star's initially rapid rotation starts slowing early on. The large initial J_2⋆ aids the inward migration, because it yields faster g frequencies, such that the orbital decay can go further before Cassini state 2 breaks (Figure 9). In this work, we did not parameterize stellar spin-down, but instead considered extremes.

In addition, if the planet is born with a primordial H/He envelope, it will undergo early mass loss through thermal mechanisms (e.g., Owen & Wu 2017; Ginzburg et al. 2018). The effect of this mass loss would only be on the level of a few percent in M_p1 (thus not affecting the g or α frequencies much), but the associated change in R_p1 could be a factor of two or more. For fixed a₁, a decrease in R_p1 would shrink α and the resulting Cassini state 2 obliquity. A shrinking R_p1 during migration could make the orbital decay go further, since it would delay the transition from Cassini state 2 to state 1. At the same time, however, the mass loss could change Q and k₂ unpredictably, as could the strong interior heating (which might lead to a partially molten planet). These effects would be interesting to expand upon in a future exploration.

Both obliquity tides and eccentricity tides are important components of the overall tidal dissipation rate (e.g., Leconte et al. 2010). There is no reason why they cannot operate simultaneously within the same system or with one dominating over the other in specific systems. However, in terms of distinguishing the obliquity tides mechanism of USP production from the corresponding eccentricity-based mechanisms (e.g., Petrovich et al. 2019; Pu & Lai 2019), we can highlight several observational predictions of our hypothesis.

First, the obliquity tides mechanism predicts that USPs should frequently have planetary companions with P < 10 days, whether they are co-transiting or not. This can be relaxed when the companions are more massive, so the condition is better stated as: if the USPs initially started with typical separations from their nearest neighbors, they should have had ∣g∣/α ≳ 1. Our theory predicts no clear trends with present-day eccentricities, whereas the eccentricity-based mechanisms may expect remnant eccentricity enhancements for the companion planets in USP systems. This would be interesting to check with observational constraints of eccentricities with radial velocities, transit timing variations (e.g., Hadden & Lithwick 2017), or stability analyses (Tamayo et al. 2020).

While our theory does not require large eccentricities, we do require nonzero inclinations. More specifically, in order for the inner planet to migrate inward while conserving angular momentum, it must transfer some of its angular momentum to other planets or the star. Typically, the star constitutes a sufficiently large sink of angular momentum when its stellar obliquity exceeds  ∼20° in the two-planet case. However, the presence of additional exterior small planets relaxes this constraint, as would the inclusion of a massive distant planet. As the USP migrates inward, it aligns with the stellar spin-axis, and it does so faster than its exterior planets. Consequently, we suggest that USPs will tend to be misaligned with their exterior planetary companions, but aligned with the stellar spin axis. The former of these features is observed (Dai et al. 2018). The latter could be tested through stellar obliquity measurements, with the larger, misaligned USP companions being more promising observational targets.

An important aspect of our theory is that USP migration is stalled when the planet breaks out of Cassini state 2 at high obliquity. Given that migration happens rapidly, relative to the total system lifetime, we predict that the majority of currently observed USPs are not undergoing significant orbital decay. This is unlike the prediction for a current high planetary obliquity of hot Jupiter WASP-12 b (Millholland & Laughlin 2018), which was hypothesized to perhaps be experiencing ongoing tidal decay from obliquity tides. An additional consequence of tidal breaking is that this mechanism preferentially produces USPs closer to P ∼ 1 day, as opposed to even shorter periods. Observationally, Sanchis-Ojeda et al. (2014) find that the USP occurrence rate peaks at 1 day and decreases with the orbital period, in agreement with our model.

As a final observational connection, we recall that the obliquity tides framework predicts that USPs should be more common around smaller-mass stars (Section 3.3). This is consistent with the empirical trend in this direction among USP stellar hosts (Sanchis-Ojeda et al. 2014).